# Efficiently compute this Table of NullSpace

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 ;;]]}] • I don't think NullSpace[] respects sparsity in the exact case, which is why it's slow for a sufficiently large size. – J. M.'s torpor Jun 15 '16 at 23:11
• @J.M. What do you mean by exact case (here I use N)? Do you see any workaround? Thank you. – anderstood Jun 16 '16 at 0:10
• Ah, sorry, I missed the N[]. In any case: did you try looking at what MatrixPower[A1,p].MatrixPower[A2,q] looks like for small n and slightly larger p and q? I suspect some severe fill-in is happening during the powering. – J. M.'s torpor Jun 16 '16 at 0:14
• @J.M. No, it does not fill-in for the following reason. $A_1$ and $A_2$ are both orthogonal with entries in $\{-\frac{1}{2},0,\frac{1}{2}\}$. Hence, they both have exactly 4 non-zero values on each row/column (each row/column is unitary). When we do the product of such matrices, the result is also such a matrix, because each element is the sum of an even number of products of halves... I am not sure if that's understandable, but I recommend MatrixPloting $A_1^pA_2^q$, it's quite interesting (the diagonal band kinds of reflects... hard to explain with words...). – anderstood Jun 16 '16 at 0:24
• @J.M. unfortunately mathematica can currently (V 10.4.1) not not exploit sparse structure for NullSpace computation. – user21 Jun 16 '16 at 6:43

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...

• You could use id = IdentityMatrix[2 n, SparseArray] - but that does not make a big difference here. – user21 Jun 17 '16 at 22:07