# Fast trace of inverse of a square matrix

Which would be the most efficient way (in computational time) to compute tr(inv(H)), where H is a (dense) square matrix?

In my particular problem I also have a LU decomposition of H already available, which was used in a previous context to solve a system of linear equations. My current approach is to simply use the LU to compute the inverse and then calculate the trace. Is there any other more efficient way to achieve this, considering I already have this matrix factorized?

I could have first computed the inverse and then made a multiplication by the inverse to solve my previous system of linear equations, but I was trying to avoid multiplication by the inverse to avoid numerical inaccuracies.

Edit: Forgot to mention that under normal conditions H should be symmetric.

-
Is there any structure inherent in your $H$ (e.g. sparse, banded, symmetric, etc.)? –  J. M. Nov 19 '10 at 0:00
The pivoting inherent in LU decomposition complicates a lot of things in what you want to do. I'm going to spitball and suggest the slightly more expensive QR decomposition, but the payoff (hopefully) is that you don't have to do the required backsubstitution on the rows of the orthogonal matrix Q in full, but only up until the component needed for the running sum of the trace. You'll have to do your own testing to be sure, of course, since the speed can be environment-dependent. –  J. M. Nov 19 '10 at 3:31
I forgot to mention that under normal conditions, H will be symmetric. But it will be always dense. –  César Nov 19 '10 at 12:17
Aha. People, people, remember, when dealing with matrix problems, please mention if there's anything special about them. Symmetry is a big thing! –  J. M. Nov 19 '10 at 12:25
Maybe it is also positive definite? This would help (no pivoting) –  Federico Poloni Nov 19 '10 at 12:58

Given that the poster has specified that his matrix is symmetric, I offer a general solution and a special case:

1. Eigendecomposition actually becomes more attractive here: the bulk of the work is in reducing the symmetric matrix to tridiagonal form, and finding the eigenvalues of a tridiagonal matrix is an O(n) process. Assuming that the symmetric matrix is nonsingular, summing the reciprocals of the eigenvalues nets you the trace of the inverse.

2. If the matrix is positive definite as well, first perform a Cholesky decomposition. Then there are methods for generating the diagonal elements of the inverse.

-
For the unsymmetric case, I am not sure, and you might want to investigate my suggestion to use QR. –  J. M. Nov 19 '10 at 12:36
Thanks, I cannot guarantee the matrix will be positive definite, but since I will probably be able to cache the matrix in a tridiagonal form prior to many calculations, this will probably speed up things considerably. –  César Nov 19 '10 at 17:56
@César: it might prove profitable to attempt a Cholesky decomposition first on your matrix to check if it is indeed positive definite, as it is cheaper to perform than the reduction to tridiagonal form. –  J. M. Nov 20 '10 at 10:52

If you have an $n\times n$ matrix $A$ with eigenvalues $\lambda_1,\ldots, \lambda_n$, then the characteristic polynomial of $A$ is $det(\lambda I-A)=\prod_{i=1}^n (\lambda-\lambda_i) = \sum_{i=0}^n a_i \lambda^i$. Then $det(A^{-1})=- a_1/a_0$, since the coefficient $a_{n-i}$ is the sum of all products of $i$ eigenvalues, so $a_1/a_0 = \sum_{i=1}^n \lambda_i^{-1}= trace(A^{-1})$. I'm not sure if it helps to compute the characteristic polynomial if you have an LU decomposition though.

-
Eigendecomposition is often a more expensive endeavor than Gaussian elimination, so I don't really see the practicality of this compared to just computing the inverse and summing the diagonal elements. –  J. M. Nov 19 '10 at 2:50
@ J.M.: I didn't mean to imply that one needs to compute the eigenvalues, just the bottom two coefficients of the characteristic polynomial. $a_0=\pm det(A)$, and $a_1$ is a sum of $n$ cofactor determinants for each diagonal entry. So the LU decomposition probably doesn't help. –  Ian Agol Nov 19 '10 at 6:09
Computing the coefficients of the characteristic polynomial is an even more unstable task than computing eigenvalues. The task of computing the coefficients of the characteristic polynomial is equivalent to finding a similarity transformation that turns an arbitrary matrix into a Frobenius companion matrix, and the matrices required for this transformation can be arbitrarily ill-conditioned. (See Wilkinson's The Algebraic Eigenvalue Problem for more details). –  J. M. Nov 19 '10 at 7:37