[Math] Calculating the inertia of a real symmetric (or tridiagonal) matrix

eigenvalues-eigenvectorslinear algebramatricesnumerical linear algebra

I'm trying to find a quick method for evaluating the inertia of a real symmetric matrix, though I don't need to evaluate eigenvalues directly. The inertia of a matrix is a triple of the number of positive eigenvalues, negative eigenvalues and eigenvalues equal to zero.

Thus far I have implemented a method of using Householder matrices to reduce a real symmetric matrix to tridiagonal form (whilst preserving inertia). A textbook by Gilbert Stewart I found through a simple Google search suggests that this is in the right direction. Does anyone have any ideas for the last step?

Best Answer

Yes, you're in the right direction. Recall that a symmetric matrix's inertia is always preserved under congruence transformations $\mathbf A\mapsto\mathbf X\mathbf A\mathbf X^\top$, and your first step of reducing to tridiagonal form via Householder transformations (with $\mathbf X$ being orthogonal) is exactly a congruence transformation.

Now, from your original matrix $\mathbf A$ you have $\mathbf A=\mathbf Q\mathbf T\mathbf Q^\top$, with $\mathbf T$ tridiagonal. From here, there are a number of methods you can use for determining the inertia of a symmetric tridiagonal matrix. In most cases, you can try computing the $\mathbf L\mathbf D\mathbf L^\top$ decomposition of your tridiagonal matrix, where $\mathbf L$ unit lower bidiagonal, and $\mathbf D$ is diagonal, from which the inertia of your tridiagonal matrix is given by the number of positive, zero, and negative elements of the diagonal matrix $\mathbf D$.

Of course, that naïve strategy won't work if the leading principal submatrix of $\mathbf T$ is singular. There is, however, an algorithm by James Bunch (see this for more details, and this for a MATLAB realization (but written for clarity at the expense of storage efficiency)), where $\mathbf T$ is decomposed as $\mathbf L\mathbf B\mathbf L^\top$, where this time $\mathbf B$ is block diagonal, with either scalar or $2\times 2$ blocks. The inertia is then evaluated from $\mathbf B$. See the links for more details.