Testing if two 12x12 matrices have the same determinant

Question

I am given a $12 \times 12$ matrix $Q$ that is symmetric, invertible, positive definite and dense. I need to test if $$\det(Q) = \det(12I-Q-J) \; \; (1)$$ where $J$ is the all ones matrix.

I am currently doing this with the armadillo library but it turns out to be too slow. The thing is that I need to do this for a trillion of matrices and it turns out that computing the two determinants is the bottleneck of my program. Hence I have two questions

1. Is there any trick I could use to compute the determinant faster given that I known their size? Is perhaps a messy expansion for $12 \times12$ matrices that could work in this case?

2. Is there any other efficient way to test the equality $(1)$

Edit. To answer the comments. I need to compute all connected non self-complementary graphs $G$ of order $13$ such that $G$ and $\overline{G}$ have the same number of spanning trees. The motivation for this can be found in this mathoverflow post. As for the machine I am running this on a 8 core 3.4GHh machine in parallel.

Edit. I was able to reduce the expected running time by 50% by making a C program for specifically computing the determinant of a $12 \times 12$ matrix. Suggestions are still welcome.

Answer 1

Since you're already using C++ and your matrices are symmetric positive definite, I would perform an unpivoted $LDL^T$ factorization of $Q$ and also of $12I-Q-J$. Here I'm assuming that $12I-Q-J$ is also positive definite, otherwise the $LDL^T$ will require pivoting for numerical stability (it's also possible that even though it's not positive definite, pivoting is not needed, but you have to try it).

This is faster than an LU factorization, and also faster than Cholesky because square roots are avoided. The determinant is simply the product of the elements of the diagonal $D$ matrix. The code to perform an LDL factorization is so simple you can write it in less than 50 lines of C. The Wikipedia page on it describes the algorithm, and I have some simple templated code to do Cholesky here. You can vastly simplify that and modify it to avoid the square root to implement the $LDL^T$ factorization.

Since you can also control the storage format, you can further optimize the routine to store only half the matrix and pack it in a linear array to maximize memory locality. I would also write simple custom dot product and rank-1 update routines since the problem sizes are so small, you should let the compiler inline the routines to reduce call overhead. Since it's a fixed size loop, the compiler should be able to automatically inline and unroll things when appropriate.

I would avoid trying to play tricks to take advantage of the fact that $12I-Q-J$ contains $Q$ inside the expression. It is likely that for such small problem sizes, these tricks end up being slower than simply performing two separate determinant computations. Of course, the only way to verify these claims is to try it.

Answer 2

Without some information about the construction of these $12\times 12$ positive definite real symmetric matrices, the suggestions to be made are of necessity fairly limited.

I downloaded the Armadillo package from Sourceforge and took a look at the documentation. Try to improve performance of separately computing $\det(Q)$ and $\det(12I - Q - J)$, where $J$ is the rank one matrix of all ones, by setting e.g det(Q,slow=false). The documentation notes that this is the default for matrices up to size $4\times 4$, so by omission I assume the slow=true option is a default for the $12\times 12$ case.

What slow=true presumably does is partial or full pivoting in getting a row echelon form, from which the determinant is easily found. However you know in advance the matrix $Q$ is positive definite, so pivoting is unnecessary for stability (at least presumptively for the bulk of your computations. It's unclear if the Armadillo package throws an exception if the pivots become unduly small, but this should be a feature of a reasonable numerical linear algebra package. EDIT: I found the Armadillo code that implements det in header file include\armadillo_bits\auxlib_meat.hpp, using C++ templates for substantial functionality. The setting slow=false doesn't appear to affect how a $12\times 12$ determinant will be done because the computation gets "thrown over a wall" to LAPACK (or ATLAS) at that point with no indication that pivoting is not required; see det_lapack and its invocations in that file.

The other point would be to follow their recommendation of building the Armadillo package linking to high speed replacements for BLAS and LAPACK, if you are indeed using those; see Sec. 5 of the Armadillo README.TXT file for details. [The use of a dedicated 64-bit version of BLAS or LAPACK is also recommended for speed on current 64-bit machines.]

Row reduction to echelon form is essentially Gaussian elimination, and has arithmetic complexity $\frac{2}{3} n^3 + O(n^2)$. For both matrices this then amounts to twice that work, or $\frac{4}{3} n^3 + O(n^2)$. These operations may well be the "bottleneck" in your processing, but there's little hope that without special structure in $Q$ (or some known relationships among the trillion test cases allowing amortization) the work could be reduced to $O(n^2)$.

For comparison, expansion by cofactors of a general $n\times n$ matrix involves $n!$ multiplication operations (and roughly as many additions/subtractions), so for $n=12$ the comparison ($12! = 479001600$ vs. $\frac{2}{3} n^3 = 1152$) clearly favors elimination over cofactors.

Another approach requiring $\frac{4}{3} n^3 + O(n^2)$ work would be reducing $Q$ to tridiagonal form with Householder transformations, which also puts $12I - Q$ into tridiagonal form. Computing $\det(Q)$ and $\det(12I - Q - J)$ can thereafter be done in $O(n)$ operations. [The effect of the rank one update $-J$ in the second determinant can be expressed as a scalar factor given by solving one tridiagonal system.]

Implementing such an independent computation might be worthwhile as a check on the results of successful (or failed) calls to Armadillo's det function.

Special Case: As suggested by a Comment of Jernej, suppose that $Q = D - J$ where $J$ as before is the (rank 1) matrix of all ones and $D=\text{diag}(d_1,\ldots,d_n)$ is a nonsingular (positive) diagonal matrix. Indeed for the proposed application in graph theory these would be integer matrices. Then an explicit formula for $\det(Q)$ is:

$$ \det(Q) = \left(\prod_{i=1}^n d_i \right)\left(1 - \sum_{i=1}^n d_i^{-1} \right) $$

A sketch of its proof affords an opportunity to illustrate wider applicability, i.e. whenever $D$ has a known determinant and the system $Dv = (1\ldots 1)^T$ is quickly solved. Begin by factoring out:

$$ \det(D - J) = \det(D) \cdot \det(I - D^{-1}J) $$

Now $D^{-1}J$ is again rank 1, namely $(d_1^{-1}\ldots d_n^{-1})^T(1\ldots 1)$. Note that the second determinant is simply:

$$ f(1) = \det(I- D^{-1}J) $$

where $f(x)$ is the characteristic polynomial of $D^{-1}J$. As a rank 1 matrix, $f(x)$ must have (at least) $n-1$ factors of $x$ to account for its nullspace. The "missing" eigenvalue is $\sum d_i^{-1}$, as may be seen from the computation:

$$ D^{-1}J\; (d_1^{-1}\ldots d_n^{-1})^T = \left(\sum d_i^{-1} \right) (d_1^{-1}\ldots d_n^{-1})^T $$

It follows that the characteristic polynomial $f(x) = x^{n-1} (x - \sum d_i^{-1})$, and $f(1)$ is as shown above for $\det(I - D^{-1}J)$, $1 - \sum d_i^{-1}$.

Also note that if $Q = D-J$, then $12I - Q - J = 12I - D + J - J = 12I - D$, a diagonal matrix whose determinant is simply the product of its diagonal entries.

Answer 3

If you have a structured way of enumerating the graphs you wish to calculate the determinants of, perhaps you could find low-rank updates which transfer you from one graph to another.

If so, then you could use the matrix determinant lemma to cheaply calculate the determinant of the subsequent graph to be enumerated using your knowledge of the current graph's determinant.

That is, for a matrix $A$ and vectors $u,v$: $$\det(A+uv^T)=(1+v^T A^{-1} u)\det(A)$$ This can be generalized if U and V are $n \times m$ matrices and $A$ is $n \times n$: $$\det(A+UV^T)=\det(I_m + V^TA^{-1}U)\det(A)$$

To efficiently calculate the inverse, you can use the Sherman-Morrison formula to obtain the inverse of the subsequent matrix from the current one: $$(A+uv^T)^{-1}=A^{-1}-\frac{A^{-1} uv^T A^{-1}}{1+v^T A^{-1} u}$$