11
$\begingroup$

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.

$\endgroup$
  • 6
    $\begingroup$ How slow is too slow? How long does it take on what hardware? Are the trillions of these $Q$s independent so that you can compute many of these determinants in parallel? If so, how big a machine can you run on? What led to this problem? Are you sure you need to compute the determinants? $\endgroup$ – Bill Barth Jun 18 '13 at 11:20
  • 3
    $\begingroup$ How often (for what fraction of cases) are the determinants the same/different? If they're different most of the time, there may be a cheaper test to determine that they may be different and you'd verify that they are the same only if the first test fails. The other way around if they're the same most of the time. $\endgroup$ – Wolfgang Bangerth Jun 18 '13 at 12:09
  • 1
    $\begingroup$ As already asked: could you provide some detail on where $\mathbf Q$ came from? Maybe there's a better approach than blindly computing determinants. $\endgroup$ – J. M. Jun 18 '13 at 14:33
  • 4
    $\begingroup$ The notion that this condition has to be tested "for a trillion of matrices" suggests 1) that $Q$ is known apriori to have some special structure (otherwise the expectation that the condition holds at random is slight), and 2) that a better approach might be to characterize all matrices $Q$ with this property (with an efficiently checkable formulation). $\endgroup$ – hardmath Jun 18 '13 at 22:06
  • 1
    $\begingroup$ @hardmath Yes, $Q$ is an integer matrix having diagonal entries ranging from $1$ to $12$ and $-1$ as off diagonal elements $\endgroup$ – Jernej Jun 20 '13 at 13:57
8
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I second the recommendation of implementing $LDL^T$, aka root-free Cholesky, since Armadillo does not seem to have a way of taking advantage of positive-definiteness/diagonal dominance. $\endgroup$ – hardmath Jun 20 '13 at 13:07
5
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Hm.. $Q$ is in fact $D-A$ where $A$ is the adjacency matrix of $G$ so I think that this result may not be correct. In particular it would imply that the number of spanning trees of a graph $G$ is determined by its degree sequence which does not hold. $\endgroup$ – Jernej Jun 21 '13 at 8:08
  • $\begingroup$ The off-diagonal entries of $Q$ will then generally contain 0 as well as -1. The $LDL^T$ decomposition suggested by Victor takes advantage of symmetry and reduces the leading term in operation count from $\frac{2}{3}n^3$ to $\frac{1}{3}n^3$. There is an exact integer approach, but that is probably not needed for your modest size matrix and entries. If I understand the construction, $12I-Q-J$ is positive definite for the same reason $Q$ is. $\endgroup$ – hardmath Jun 21 '13 at 11:57
  • $\begingroup$ @Jernej: If you believe something I stated is incorrect, I've created a chat room based on this Question where the discussion can be threaded without unnecessary commenting here. $\endgroup$ – hardmath Jun 22 '13 at 13:36
1
$\begingroup$

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}$$

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.