Simultaneous Iteration:
Real Sym

if we think of the subspace as a list of basis vectors then the following code implements one step of Simultaneous iteration on the list of vectors.

SimIter[A_,Q_]:=Orthogonalize[Table[A.Q[[i]],{i,1,Length[Q]}]];

{n,k}={23,5};
A=RandomReal[{-1,1},{n,n}];
ASym=A.Transpose[A];
Q=Orthogonalize[RandomReal[{-1,1},{k,n}]];
Do[
Q=SimIter[ASym,Q],
{34}]

{q1,q2}=Q[[1;;2]]
Eigenvectors[ASym,2]

To compute the eigenvectors is easy.

ASym.q1/q1
ASym.q2/q2

Eigenvalues[ASym,2]

Simultaneous Iteration:
Non-Sym

For a non-sym matrix it is less obvious that it has worked.

SimIter[A_,Q_]:=Orthogonalize[Table[A.Q[[i]],{i,1,Length[Q]}]];

{n,k}={23,5};
A=RandomReal[{-1,1},{n,n}];
ANonSym=A;
Q=Orthogonalize[RandomReal[{-1,1},{k,n}]];
Do[
Q=SimIter[ANonSym,Q],
{240}]

{q1,q2}=Q[[1;;2]]
{qq1,qq2}=Eigenvectors[ANonSym,2]

And it is less easy to compute the eigenvalues if they are complex.

ANonSym.q1/q1
ANonSym.q2/q2

Eigenvalues[ANonSym,2]

The trick is to use Theorem 6.1.4, build the little matrix BHat, an compute the evals.

BHat={
{(ANonSym.q1).q1, (ANonSym.q1).q2},
{(ANonSym.q2).q1, (ANonSym.q2).q2}
};
MatrixForm[BHat]
Eigenvalues[BHat]