# Determinant of Large Sparse Matrix

I've looked through previous posts and although there appear to be somewhat similar questions to the ones I'm asking here, they aren't quite the same. Nevertheless, if there is a relevant question(s) that I've missed I apologize for the duplication.

I'm interested in computing the determinant of a large sparse matrix ($80\times80$) to be specific. In a somewhat predictable fashion, if the non-zero entries are kept as formal variables, Mathematica runs and I have yet to see it finish. Some of the formal variables can be swapped out for floating point numbers (not all of them however—there are 16 variables which have to be kept as variables), and after the swap Mathematica runs and finishes in about 1.5 mins, but the output is unwieldy (to say the least), moreover, there are many many 0. expressions appearing in the output (which I believe are contributing substantially to the unwieldy nature of the output).

I have reason to believe that many (if not all) of the 0. expressions should indeed be simply 0's, this is because the formal variable entries in the matrix which can be swapped out for real numbers come from a using NIntegrate, and there are really only a handful (8) distinct Integrals which get numerically integrated while there are ~600 non-zero entries.

I have a few things I'd like to know and would really appreciate any input anyone has:

• What is Mathematica's implementation of the determinant? (it appears to involve Gaussian elimination (?) as the output is a gigantic quotient, but perhaps its something else...)

• I would love to keep the speed of the computation of the determinant which uses the real numbers, but I'd also like for mathematica to instead output integer zeroes and not real zeroes when a cancellation happens, is there a way to do this? (I've naively just tried Chop[Det[A]], however it seems as though somewhere in the denominator Mathematica chops a little to much and outputs warnings regarding Infinities and then stops the computation)

• Is there a way to simply speed up Mathematica (for example: in the way Compile together with CompilationTarget -> "C" can speed up numerical computations substantially), while carrying the 16 formal variables?

Again, any help/new perspective would be much appreciated.

Edit: This is the code I use to populate the matrix. The following lines basically define all the parameters and variables associated with the matrix in question.

s = 1/4; p = 4; q = 4;
k = {
{1, 0},
{0, -1},
{-1, 0},
{0, 1}
};
\[CapitalDelta] = {-1, 1};
\[CapitalLambda] = Array[Subscript[\[Lambda], ##] &, {4, 4}];
Array[Subscript[\[ScriptI], ##] &, {1, 2, 2}];
\[ScriptCapitalI] =
{
{
{Subscript[\[ScriptI], 1, 1, 1], Subscript[\[ScriptI], 1, 1, 2]},
{Subscript[\[ScriptI], 1, 2, 1], Subscript[\[ScriptI], 1, 2, 2]}
},
{
{Subscript[\[ScriptI], 1, 2, 2], Subscript[\[ScriptI], 1, 2, 1]},
{Subscript[\[ScriptI], 1, 1, 1], Subscript[\[ScriptI], 1, 1,2]}
},
{
{Subscript[\[ScriptI], 1, 1, 1], Subscript[\[ScriptI], 1, 1, 2]},
{Subscript[\[ScriptI], 1, 2, 2], Subscript[\[ScriptI], 1, 2,1]}
},
{
{Subscript[\[ScriptI], 1, 2, 1], Subscript[\[ScriptI], 1, 2, 2]},
{Subscript[\[ScriptI], 1, 1, 1], Subscript[\[ScriptI], 1, 1, 2]}
},
{
{s, s}, {s, s}
}
};
Array[Subscript[\[ScriptJ], ##] &, {1, 2, 2}];
\[ScriptCapitalJ] =
{
{
{Subscript[\[ScriptJ], 1, 1, 1], Subscript[\[ScriptJ], 1, 1, 2]},
{Subscript[\[ScriptJ], 1, 2, 1], Subscript[\[ScriptJ], 1, 2, 2]}
},
{
{-Subscript[\[ScriptJ], 1, 2, 2], -Subscript[\[ScriptJ], 1, 2, 1]},
{Subscript[\[ScriptJ], 1, 1, 1], Subscript[\[ScriptJ], 1, 1, 2]}
},
{
{Subscript[\[ScriptJ], 1, 1, 1], Subscript[\[ScriptJ], 1, 1, 2]},
{-Subscript[\[ScriptJ], 1, 2, 2], -Subscript[\[ScriptJ], 1, 2, 1]}
},
{
{Subscript[\[ScriptJ], 1, 2, 1], Subscript[\[ScriptJ], 1, 2, 2]},
{Subscript[\[ScriptJ], 1, 1, 1], Subscript[\[ScriptJ], 1, 1, 2]}
},
{
{0, 0}, {0, 0}
}
};
A = Table[0, {i, 1, 5 p q}, {j, 1, 5 p q}];


The following lines replace the symbolic vairables with real numbers (or real numbers multiplied by a variable).

(*  Fill in the matrices \[ScriptCapitalI] and \[ScriptCapitalJ]  *)

Do[
Do[
Subscript[\[ScriptI], i, 1, j] =
Round[NIntegrate[
Log[(x - s k[[i, 1]])^2 + ((\[CapitalDelta][[j]] s)/2 +
s k[[i, 2]])^2], {x, -s/2 , s/2}], 10^-6];
Subscript[\[ScriptI], i, 2, j] =
Round[NIntegrate[
Log[(x - s k[[i, 2]])^2 + ((\[CapitalDelta][[j]] s)/2 +
s k[[i, 1]])^2], {x, -s/2, s/2}], 10^-6];
Subscript[\[ScriptJ], i, 1, j] =
Round[NIntegrate[-(((\[CapitalDelta][[j]] s)/2 +
s k[[i, 2]])/((x -
s k[[i, 1]])^2 + ((\[CapitalDelta][[j]] s)/2 +
s k[[i, 2]])^2)), {x, -s/2, s/2}], 10^-6];
Subscript[\[ScriptJ], i, 2, j] =
Round[NIntegrate[-(((\[CapitalDelta][[j]] s)/2 +
s k[[i, 1]])/(((\[CapitalDelta][[j]] s)/2 +
s k[[i, 1]])^2 + (x - s k[[i, 2]])^2)), {x, -s/2, s/2}],
10^-6],
{j, 1, 2}],
{i, 1, 4}];


Then we fill in the 80x80 matrix, A.

(*  Fill in the first fifth of A  *)
Do[
Do[
Do[
i = 2 p (l2dash - 1) + 2 (l1dash - 1) + d + 1;
Do[
Do[
j = 5 p (l2  - 1) + 5 (l1 - 1);
Which[
(*---*){l1, l2, d} == {l1dash, l2dash, 0},
A[[i, j + 1 ;; j + 5]] = \[ScriptCapitalI][[;; , 1, 1]],

(*---*){l1, l2, d} == {l1dash, l2dash + 1, 0},
A[[i, j + 1 ;; j + 5]] = -\[ScriptCapitalI][[;; , 1, 2]],

(*---*){l1, l2, d} == {l1dash, l2dash, 1},
A[[i, j + 1 ;; j + 5]] =
Subscript[\[Lambda], l1, l2] \[ScriptCapitalJ][[;; , 1, 1]],

(*---*){l1, l2, d} == {l1dash, l2dash + 1, 1},
A[[i,
j + 1 ;; j + 5]] = -Subscript[\[Lambda], l1,
l2] \[ScriptCapitalJ][[;; , 1, 2]],

(*---*)True,
A[[i, j + 1 ;; j + 5]] = {0, 0, 0, 0, 0}
]
(*  End of main block  *)
, {l1, 1, p}]
, {l2, 1, q}]
, {d, 0, 1}]
, {l1dash, 1, p}]
, {l2dash, 1, q - 1}];

(*  Fill in second fifth of A  *)
Do[
Do[
Do[
i = 2 p (q - 1) + 2 (l2dash - 1) + 2 q (l1dash - 1) + d + 1;
Do[
Do[
j = 5 p (l2  - 1) + 5 (l1 - 1);
Which[
(*---*){l1, l2, d} == {l1dash, l2dash, 0},
A[[i, j + 1 ;; j + 5]] = \[ScriptCapitalI][[;; , 2, 1]],

(*---*){l1, l2, d} == {l1dash + 1, l2dash , 0},
A[[i, j + 1 ;; j + 5]] = -\[ScriptCapitalI][[;; , 2, 2]],

(*---*){l1, l2, d} == {l1dash, l2dash, 1},
A[[i, j + 1 ;; j + 5]] =
Subscript[\[Lambda], l1, l2] \[ScriptCapitalJ][[;; , 2, 1]],

(*---*){l1, l2, d} == {l1dash + 1, l2dash, 1},
A[[i,
j + 1 ;; j + 5]] = -Subscript[\[Lambda], l1,
l2] \[ScriptCapitalJ][[;; , 2, 2]],

(*---*)True,
A[[i, j + 1 ;; j + 5]] = {0, 0, 0, 0, 0}
]
(*  End of main block  *)
, {l1, 1, p}]
, {l2, 1, q}]
, {d, 0, 1}]
, {l1dash, 1, p - 1}]
, {l2dash, 1, q}];

(* Fill in third fifth of A *)
Do[
Do[
row = 4 p q - 2 (p + q) + p (j - 1) + i;
col = 5 p (j - 1) + 5 (i - 1);
Do[
Which[
(*---*)(k >= col + 1) && (k <= col + 4),
A[[row, k]] = 1,
(*---*)True,
A[[row, k]] = 0
]
, {k, 1, 5 p q }]
, {j, 1, q}]
, {i, 1, p}];

(*  Fill in fourth fifth of A  *)
Do[
Do[
i = 2 (p (q - 1) + q (p - 1)) + p q + 2 (l1dash - 1) + d;
Do[
Do[
j = 5 p (l2  - 1) + 5 (l1 - 1);
Which[
(*---*){l1, l2, d} == {l1dash, 1, 1},
A[[i, j + 1 ;; j + 5]] = \[ScriptCapitalJ][[;; , 1, 2]],

(*---*){l1, l2, d} == {l1dash, q, 2},
A[[i, j + 1 ;; j + 5]] = \[ScriptCapitalJ][[;; , 1, 1]],

(*---*)True,
A[[i, j + 1 ;; j + 5]] = {0, 0, 0, 0, 0}
]
(*  End of main block  *)
, {l1, 1, p}]
, {l2, 1, q}]
, {d, 1, 2}]
, {l1dash, 1, p}];

(*  Fill in last fifth of A  *)
Do[
Do[
i = 2 (p q + q (p - 1)) + p q + 2 (l2dash - 1) + d;
Do[
Do[
(*  Enter code here ...  *)
j = 5 p (l2  - 1) + 5 (l1 - 1);
Which[
(*---*){l1, l2, d} == {p, l2dash, 1},
A[[i, j + 1 ;; j + 5]] = \[ScriptCapitalJ][[;; , 2, 1]],

(*---*){l1, l2, d} == {1, l2dash, 2},
A[[i, j + 1 ;; j + 5]] = \[ScriptCapitalJ][[;; , 2, 2]],

(*---*)True,
A[[i, j + 1 ;; j + 5]] = {0, 0, 0, 0, 0}
]
(*  End of main block  *)
, {l1, 1, p}]
, {l2, 1, q}]
, {d, 1, 2}]
, {l2dash, 1, q}];

A[[5 p q, 5 p (q - 1) + 1 ;; 5 p (q - 1) + 5]] = {0, 0, 0, 0, 1};
A = SparseArray[A];

• You could try applying Rationalize to the matrix to replace the floating point numbers with exact values. – mikado Sep 12 '16 at 22:22
• I tried doing something which I think should be effectively the same (using Round[ ,10^(-5)]), however, it seems like mathematica treats these exact numbers in a similar way to how it treats symbolic variables ... and the program doesn't finish – Andrew Ross Sep 12 '16 at 23:05
• It might be more useful if you could actually show us the matrix. – Igor Rivin Sep 13 '16 at 11:00
• Do the same variables appear in multiple places? That is to say, is there some low bound on degree of the result with respect to each variable? – Daniel Lichtblau Sep 13 '16 at 15:19
• @IgorRivin I'll add the code to create the matrix. – Andrew Ross Sep 13 '16 at 15:23

this ultimately doesn't work sucessfully, but might lead somewhere. Simplify the matrix as much as possible by identifying repeated elements:

AA = SparseArray[(Normal@A) /.
Select[Flatten[
MapIndexed[ {# -> a[First@#2], -# -> -a[First@#2]} &,
Union[Flatten[Normal@A]]], 1],
! (#[[1]] === 0 || #[[1]] === 1 || #[[1]] ===
1/4 || #[[1]] === -1/4) &]] ;


at least now we can see there are ~50 unique elements and we can visualize what it looks like.

MatrixForm[AA]


I let Det run on this quite a while without a result.

• This is interesting. It think I might fiddle around with something like this tomorrow. – Andrew Ross Sep 15 '16 at 3:08