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 ;;]]}]
NullSpace[]
respects sparsity in the exact case, which is why it's slow for a sufficiently large size. $\endgroup$ – J. M.'s ennui♦ Jun 15 '16 at 23:11N
)? Do you see any workaround? Thank you. $\endgroup$ – anderstood Jun 16 '16 at 0:10N[]
. In any case: did you try looking at whatMatrixPower[A1,p].MatrixPower[A2,q]
looks like for smalln
and slightly largerp
andq
? I suspect some severe fill-in is happening during the powering. $\endgroup$ – J. M.'s ennui♦ Jun 16 '16 at 0:14MatrixPlot
ing $A_1^pA_2^q$, it's quite interesting (the diagonal band kinds of reflects... hard to explain with words...). $\endgroup$ – anderstood Jun 16 '16 at 0:24NullSpace
computation. $\endgroup$ – user21 Jun 16 '16 at 6:43