Let V = C^m and G = span{G1,…,Gn} ⊆ End(V). For a subspace I ⊆ V and a linear subspace S ⊆ G, “I is invariant under every S ∈ S” means S(I) ⊆ I for all S ∈ S.

Key observation (reduction to a rank condition)

• Fix k with 1 ≤ k ≤ m−1 and a k-dimensional subspace I ⊆ V.
• Let i_I: I → V be the inclusion and q_I: V → V/I the quotient map.
• Define the linear map Φ_I: G → Hom(I, V/I) by Φ_I(A) = q_I ∘ A ∘ i_I.

Then:

• I is invariant under A ∈ G iff Φ_I(A) = 0.
• The largest subspace of G whose every element leaves I invariant is S(I) = ker Φ_I.
• Therefore, the admissible pairs (S, I) are exactly those with S ⊆ ker Φ_I. In particular, to get all “maximal” pairs it suffices to list all I with ker Φ_I ≠ {0} and take S = ker Φ_I (any subspace of ker Φ_I also works).

This reduces the problem to finding all I with ker Φ_I ≠ 0 (equivalently rank Φ_I < dim G = n). For a fixed k, this set is a determinantal locus in the Grassmannian Gr(k, V).

Concrete computational recipe

• Fix k. Work on an affine chart of Gr(k, m). Represent I by an m×k matrix U of full column rank (a basis), and choose an (m−k)×m matrix W of full row rank with WU = 0 (its rows span a complement to I). For instance, in the standard chart one may take U = [I_k; Y] and W = [−Y^T, I_{m−k}], with Y an (m−k)×k matrix of coordinates.
• For each i = 1,…,n, form the (m−k)×k matrix B_i(U) = W G_i U. Note that I is invariant under A = ∑ α_i G_i iff there exists α = (α_i) ≠ 0 such that ∑ α_i B_i(U) = 0.
• Stack the B_i(U) as columns of a linear map B(U): C^n → Hom(C^k, C^{m−k}) by α ↦ ∑ α_i B_i(U). Then ker Φ_I ≠ 0 iff ker B(U) ≠ 0, i.e. rank B(U) < n.
• Algebraically: the locus of I with ker Φ_I ≠ 0 is defined on the chart by the vanishing of all n×n minors of the matrix whose columns are the vectorizations vec(B_i(U)). These minors are polynomial in the chart coordinates (entries of Y).
• Solve these determinantal equations (for each k) to obtain the admissible I. For each such I (given by U), compute S(I) = {∑ α_i G_i : ∑ α_i B_i(U) = 0}, i.e. S(I) = ker B(U) lifted back to G.

What this gives you

• A complete characterization: all pairs (S, I) with S ⊆ ker Φ_I, and the “maximal” S for each I is ker Φ_I.
• A practical method: for each k, find the determinantal locus {I ∈ Gr(k, m): rank Φ_I < n} by vanishing of minors; for each I in it, compute ker Φ_I and thereby all S.

Additional remarks

• Dimension control: If you want S to have dimension at least d, require rank Φ_I ≤ n − d. On charts this means all (n − d + 1)×(n − d + 1) minors vanish.
• Generic behavior: For fixed k, the linear constraints AU ⊆ I amount to k(m−k) scalar conditions. Consequently, when n ≤ k(m−k), generically the rank of Φ_I equals n for most I, so ker Φ_I = 0 except possibly on a proper algebraic subset of the Grassmannian. When n > k(m−k), generically dim ker Φ_I ≥ n − k(m−k) for every I.
• Complexity: Solving the determinantal equations on Grassmannians can be computationally heavy in general (there can be infinitely many I). The formulation above shows the problem is algebraic; it also explains why there is no simple “finite list” in general.
• Relation to your observation: Your “fix I and compute the largest S” is exactly S(I) = ker Φ_I. The simultaneous viewpoint packages all I together as the determinantal loci on the Grassmannians.

In short: All admissible pairs are obtained by taking I in the determinantal locus {I: rank Φ_I < n} and S any subspace of ker Φ_I. The locus and the kernels are computable from the given G_i via the recipe above.