Efficiently compute this Table of NullSpace

Question

I have two $(2n,2n)$ matrices, $A_1$ and $A_2$, and I would like to compute $$\ker(A_1^p A_2^q -I)$$ for $p,q\leq 2n$.

Both matrices are orthogonal and have exactly four non-zeros values on each line and columns (sparse matrices). This is what a typical $A_1$ looks like, for $n=5$:

$$A_1=\dfrac{1}{2}\ \left( \begin{array}{cccccccccc} -1 & -1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & -1 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 \\ -1 & 1 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & -1 & 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & -1 & 1 & 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 & -1 & 1 & 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & -1 \\ 0 & 0 & 0 & 0 & 0 & 0 & -1 & 1 & 1 & -1 \\ \end{array} \right)$$

Here is the code to produce $A_1$ and $A_2$:

n = 5;
mat = 1/2*{{1, -r,0,0,1,r},{-1/r, 1, 0, 0, 1/r,1}} ;
A1 = ConstantArray[0,{2n,2n}];
Table[A1[[i;;i+1,i-2;;i-2+5]] = mat,{i,3,(2n-2),2}];
A1[[;;2,;;6]] = 1/2*{{-1, -r, 1, r,0,0},{1/r, 1,1/r,1,0,0}};
A1[[-2;;,-4;;]] = 1/2*{{1,-r,1,-r},{-r,1,r,-1}};
A2 = A1 + SparseArray[{{1,1}->1,{2,2}->-1,{1,2}->r,{2,1}->-1/r},{2n,2n}];
r = 1;
A1 = SparseArray[A1//N]//Chop;
A2 = SparseArray[A2//N]//Chop;

The problem I face is that I am intersted in these kernels for matrices corresponding to $n\approx 500$, and doing:

Table[NullSpace[MatrixPower[A1,p].MatrixPower[A2,q] -
        IdentityMatrix[2n]], {p,1,2n}, {q,1,2 n}]

takes a long time... I tried precomputing all the matrix powers to avoid computing them several times, but it was even worse, probably because of the memory calls (despiste the SparseArray structure).

Any idea how I could speed up this computation?

Note that I also tried to simplify the problem on the math side (see this math.SE post) without success.

Note that $A_1^p A_2^q$ remains a sparse matrix (see comments), as it can be illustrated:

tab = Table[MatrixPower[A1, i].MatrixPower[A2, 12] 
  // MatrixPlot[#, ImageSize -> Small] &, {i, 1, 50, 5}]
GraphicsGrid[{tab[[;; 5]], tab[[6 ;;]]}]

Answer 1

The most I've managed to improve the speed is about a factor of 2, but I thought I would share my attempt anyhow.

First, let's just compute the identity matrix once, instead of once each iteration step:

id = SparseArray[IdentityMatrix[2n]];

Second, since we are going to need all the powers of both A1 and A2, we can gain some speed by simply bumping the power on each iteration in the Table. You can see that this approach makes sense by comparing the timing of MatrixPower[A1,10] to MatrixPower[A1,100].

Third, we can use QRDecomposition since it seems to be faster than all the other decompositions (LU, SVD, RowReduce), at least on my machine. The null space should be spanned by the components of Q corresponding to tiny values on the diagonal of R.

The code:

Block[{a1p,a2p,tab,q,r},
a1p = A1; a2q = A2;
 nullspacebases = Table[
    tab = Table[
      {q, r} = QRDecomposition[a1p.a2q - id];
      a2q = a2q.A2;
      Pick[q, Chop[Diagonal[r]], 0],
      {q, 2 n}
      ];
    a2q = A2;
    a1p = a1p.A1;
    tab,
    {p, 2 n}
    ];
 ]

I've only been able to test this for n = 80, so there is still some patience needed :p Note that for n = 80 the Table in the OP which uses NullSpace gives SingularValueDecomposition::cfail and NullSpace::cfail messages, while QRDecomposition gives none, so I'm not sure what to trust here...