# Calculating determinant while solving $Ax=b$ using CG

I am solving $Ax=b$ for a huge sparse positive definite matrix $A$ using the conjugate gradient (CG) method. It is possible to compute the determinant of $A$ using the information produced during the solve?

• Why would you wish to compute the determinant? Such a result will surely be either an underflow or an overflow for a huge matrix anyway. I'd be more charitable had you asked to compute the condition number, but don't waste your time on the determinant!
– user840
Dec 20 '12 at 11:52
• You probably already know that, but the Ritz values during the conjugate gradient process converge to the eigenvalues of the matrix, and you can derive simple estimates for the determinant from this. Dec 20 '12 at 21:09

Computing the determinant of a sparse matrix is typically as expensive as a direct solve, and I am skeptical that CG would be of much help in computing it. It would be possible to run CG for $n$ iterations (where $A$ is $n \times n$) in order to generate information for the entire spectrum of $A$, and to then compute the determinant as the product of the eigenvalues, but this would be both slow and numerically unstable.

It would be a better idea to compute the sparse-direct Cholesky factorization of your matrix, say $A = L L^H$, where $L$ is lower-triangular. Then $$\text{det}(A) = \text{det}(L) \text{det}(L^H) = |\text{det}(L)|^2,$$ where $\text{det}(L)$ is simply the product of the diagonal entries of the lower-triangular matrix $L$ since the eigenvalues of a triangular matrix lie along its diagonal.

In the case of a general non-singular matrix, a pivoted LU decomposition should be used, say $PA=LU$, where $P$ is a permutation matrix, so that $$\text{det}(A) = \text{det}(P^{-1}) \cdot \text{det}(L) \cdot \text{det}(U).$$ Since $P$ is a permutation matrix, $\text{det}(P)=\pm 1$, and, by construction, $L$ will typically have a diagonal of all ones, which implies that $\text{det}(L)=1$. You can thus compute $\text{det}(A)$ as $\pm \text{det}(U)$ and again recognize that the determinant of a triangular matrix is simply the product of its diagonal entries. Thus, the cost of computing the determinant is essentially just that of a factorization.

• This would be one of the possibilities (although I would use cholesky factorization) if the matrix is smallish, however $A$ has size ~$10^6 x 10^6$ and hence it is not possible to do a decomposition Dec 19 '12 at 17:48
• @ManuelSchmidt Sparse matrices of that size resulting from finite-element type discretizations can usually be easily factored with (for instance) multifrontal methods. I agree that Cholesky factorization should be used if your matrix is HPD (and the generalization of my above argument is obvious). Dec 19 '12 at 18:45
• Thanks for your fast answer & reply. Unfortunately the matrix has no spezial structure (which would allow a easy factorization). Dec 19 '12 at 19:33
• I am curious as to why you need to compute the determinant of the matrix. Are the highest and lowest eigenvalues not sufficient? Dec 19 '12 at 19:52
• It is part of a complex probability distribution function and not only a normalization constant. I know distributions can be factored (and thats what we are doing at the moment) however we have tons of data to model and each of the factors become huge. Dec 20 '12 at 6:09

Others have already noted this but I think it's still worth pointing out that the determinant is not a useful quantity almost always when your matrix is large. The problem is that large matrices are most often approximations to things that are even larger dimensional (statistical samples of large populations, finite dimensional approximations to infinite dimensional things like partial differential equations, etc). Let us say your matrix is $A$ and that it, in some way, approximates some process $B$ where $\textrm{dim}A \ll \textrm{dim} B$ and possibly $\textrm{dim}B=\infty$.

Now, the operator $B$ will also have a spectrum, and the spectrum of $A$ will typically be related in some sense to the spectrum of $B$ -- e.g., the eigenvalues of $A$ approximate (a subset of) those of $B$. But this does not mean that $\textrm{det}A$ approximates $\textrm{det}B$. In fact, if, for example, the eigenvalues of $A$ approximate some of the eigenvalues of $B$, then $$\textrm{det}A = \prod_{j=1\ldots \textrm{dim}A} \lambda_i(A) \approx \prod_{j=1\ldots \textrm{dim}A} \lambda_i(B) \prod_{j=\textrm{dim}A+1\ldots \textrm{dim}B} \lambda_i(B)$$ but the second factor in this product is typically either the product of the small or of the large eigenvalues of $B$ that aren't resolved by $A$. In other words, this factor is typically either very small (possibly zero, if $\textrm{dim}B=\infty$) or very large (possibly infinite). In other words, in all of these cases, we typically have that $\textrm{det}A$ does not approximate $\textrm{det}B$ and is not usually a useful thing to compute.

Without getting (again) into why and how determinants are evil, let's assume that your operator is either not easily factorizable or simply not available as a matrix at all and that you really need to estimate its determinant.

One way to obtain an estimate of the determinant is to use an implementation of CG that is directly based on the Lanczos process. See for instance Algorithm 6.17 (D-Lanczos) in Saad's book "Iterative Methods for Sparse Linear Systems" (SIAM), page 189. There, the operator $A$ is semi-explicitly factored and you can use the determinant of the factors as an approximation of the determinant of $A$. Let me emphasize that I've never tried to estimate a determinant and I've no idea whether or not what I'm suggesting is a good idea---problably not given all the trouble ghost eigenvalues give us. However, that's what comes to mind when reading your question.

You can probably reverse-engineer how this estimate of the determinant comes about in the standard implementation of CG by following closely Section 6.7.3 of the book.

Determinant of matrix A can be computed as $$det(A) = \prod_{i = 1}^n \alpha_k^{-1},$$ where $\alpha_k = \cfrac{r_k^Tr_k}{p_k^TAp_k}$ - are steps computed during CG iteration. But this works only if $r_k \ne 0 \;\,\forall k=1, \ldots, n$. The proof is the following. Let $R$ be a matrix consisting of vectors $r_k$ and let $P$ be a matrix consisting of $p_k$. $$p_k = -r_k + \sum_{i=1}^{k-1}\gamma_i r_i.$$ So, by the property of determinant $det(P) = (-1)^n det(R)$. Now notice that vectors $r_k$ are orthgonal and vectors $p_k$ are conjugate with respect to $A$. Therefore $$\prod_{k = 1}^n \alpha_k = \prod_{k = 1}^n \frac{r_k^T r_k}{p_k^TAp_k} = \frac{det(R^TR)}{det(P^TAP)} = \frac{det(R^TR)}{det(A)det(P^TP)}= (det(A))^{-1}.$$