# inertia count sparse matrix with dense low-rank perturbation

I would like to determine the number of negative eigenvalues (inertia count) of the $(N \times N)$ symmetric real matrix $K - \sigma M$, with $K$ a positive-definite sparse matrix and $M$ a positive-semidefinite matrix such that $M = Z Z^T$ with $Z$ a dense $(N \times r)$ full column-rank matrix, $r \ll N$, $\sigma > 0$.

The matrix $K - \sigma M$ is a low-rank perturbation of $K$ and cannot be assembled into memory because $N$ is very large. If $M$ was sparse, the LDL factorization $K - \sigma M = L D L^T$ would be possible and would give the inertia count through the inertia count of $D$.

• Even if $N$ is huge, can you factorize $K$? – user7440 Aug 20 '18 at 3:02
• Yes, I can. It is not possible to factorize $K - \sigma M$ because it $M$ is dense and cannot be assembled. I edited and replaced 'huge' with 'very large'. – Olivier Aug 21 '18 at 4:58
• If you can factor K and store/apply Z, you should be able to use the Woodbury identity to build the solution operator (inverse) for K-sigma*M. This isn't quite the finish line, but could be a start because it lets you probe for near-zero eigenvalues using inverse iteration. – rchilton1980 Aug 21 '18 at 12:53
• Yes, thanks to Woodbury matrix identity I am able to efficiently solve for $\textbf{x}$ the linear system ${( K - \sigma M )}^{-1} \textbf{f}$ for a given real vector $\textbf{f}$. What I would like to do is to solve the generalized eigenvalue problem $K \boldsymbol{\varphi} = \lambda M \boldsymbol{\varphi}$ for many eigenpairs by spectrum slicing. To do so, I need to divide the sprectum in subintervals for which I know the number of eigenvalues present. An accurate inertia count at several $\sigma$'s would be perfect; a rough approximation of the spectral density would be better than nothing. – Olivier Aug 21 '18 at 18:02

Upon further examination, I do think the Woodbury identity can be used to solve this problem. With it we can write:

$\left( \mathbf K - \sigma \mathbf Z \mathbf Z^T \right)^{-1} = \mathbf K^{-1} - \mathbf K^{-1} \mathbf Z \left(-\frac{1}{\sigma}\mathbf I + \mathbf Z^{T} \mathbf K^{-1} \mathbf Z\right)^{-1}\mathbf Z^T \mathbf K^{-1}$

Since $\mathbf K$ is symmetric positive definite, you can write $\mathbf K^{-1}$ in terms of the Cholesky triangle $\mathbf L$:

$\mathbf K^{-1} = \left(\mathbf L \mathbf L^T \right)^{-1} = \mathbf L^{-T} \mathbf L^{-1}$

Insert this representation, then simplify by introducing $\mathbf Y = \mathbf L^{-1} \mathbf Z$:

$\left( \mathbf K - \sigma \mathbf Z \mathbf Z^T \right)^{-1} = \mathbf L^{-T} \mathbf L^{-1} - \mathbf L^{-T} \mathbf Y \left(-\frac{1}{\sigma}\mathbf I + \mathbf Y^{T} \mathbf Y\right)^{-1}\mathbf Y^T \mathbf L^{-1}$

Now regroup:

$\left( \mathbf K - \sigma \mathbf Z \mathbf Z^T \right)^{-1} = \mathbf L^{-T} \left[ \mathbf I - \mathbf Y \left(-\frac{1}{\sigma}\mathbf I + \mathbf Y^{T} \mathbf Y\right)^{-1}\mathbf Y^T \right] \mathbf L^{-1}$

This is progress. It's not quite an $\mathbf L \mathbf D \mathbf L^T$ decomposition, but it might be good enough in the sense that it's a symmetric real congruence transform of the perturbed $\mathbf K$. So the square-bracketed part, despite being non-diagonal, still embeds the inertia. To proceed, introduce the matrix $\mathbf S = \left(-\frac{1}{\sigma}\mathbf I + \mathbf Y^{T} \mathbf Y\right)$, which we should be able to tabulate explicitly because it's so small (only $r \times r$). We can also readily compute/apply its inverse (using dense $\mathbf L \mathbf D \mathbf L^T$):

$\left( \mathbf K - \sigma \mathbf Z \mathbf Z^T \right)^{-1} = \mathbf L^{-T} \left[ \mathbf I - \mathbf Y \mathbf S^{-1}\mathbf Y^T \right] \mathbf L^{-1}$

Continue teasing at the bracketed part to get it to display its eigenvalues. Compute the $\mathbf Q \mathbf R$ decomposition of the tall/skinny matrix $\mathbf Y$:

$\left( \mathbf K - \sigma \mathbf Z \mathbf Z^T \right)^{-1} = \mathbf L^{-T} \left[ \mathbf I - \left(\mathbf Q \mathbf R\right)\mathbf S^{-1} \left( \mathbf Q \mathbf R \right)^T \right] \mathbf L^{-1}$

Regroup, then compute the eigendecomposition of the (symmetric) matrix $\mathbf R \mathbf S^{-1} \mathbf R^T = \mathbf V \mathbf \Lambda \mathbf V^T$:

$\left( \mathbf K - \sigma \mathbf Z \mathbf Z^T \right)^{-1} = \mathbf L^{-T} \left( \mathbf I - \left(\mathbf Q \mathbf V\right) \mathbf \Lambda \left( \mathbf Q \mathbf V \right)^T \right) \mathbf L^{-1}$

The inner data structure, $\left(\mathbf Q \mathbf V\right) \mathbf \Lambda \left( \mathbf Q \mathbf V \right)^T$, is like a "thin" eigendecomposition. For sake of illustration, fatten it back up by including the implicit zero eigenpairs. Group the handful of non-null eigenvectors into the orthogonal matrix $\mathbf U = \mathbf Q \mathbf V$ and denote the null eigenvectors as $\mathbf N = \text{null}(\mathbf U^T)$:

$\left( \mathbf K - \sigma \mathbf Z \mathbf Z^T \right)^{-1} = \mathbf L^{-T} \left( \mathbf I - \left[\begin{array}{cc} \mathbf U & \mathbf N \end{array}\right] \left[\begin{array}{cc} \mathbf \Lambda & \mathbf 0 \\ \mathbf 0 & \mathbf 0 \end{array}\right] \left[\begin{array}{cc} \mathbf U & \mathbf N \end{array}\right]^T \right) \mathbf L^{-1}$

Of course $\mathbf I$ can be simultaneously diagonalized with any other matrix, so fold it right into this structure:

$\left( \mathbf K - \sigma \mathbf Z \mathbf Z^T \right)^{-1} = \mathbf L^{-T} \left( \left[\begin{array}{cc} \mathbf U & \mathbf N \end{array}\right] \left[\begin{array}{cc} \mathbf I - \mathbf \Lambda & \mathbf 0 \\ \mathbf 0 & \mathbf I \end{array}\right] \left[\begin{array}{cc} \mathbf U & \mathbf N \end{array}\right]^T \right) \mathbf L^{-1}$

This displays eigenvalues/inertia. Since the perturbation was only low rank, most of them will stay positive as indicated the lower right $\mathbf I$ block. However, from any $\mathbf \Lambda_i$ that is greater than one, $\left( \mathbf K - \sigma \mathbf Z \mathbf Z^T \right)^{-1}$ will pick up a negative eigenvalue from the upper left $\mathbf I - \mathbf \Lambda$ block. Since the eigenvalues of $\mathbf K - \sigma \mathbf Z \mathbf Z^T$ are just the reciprocals, the signs/inertia are the same.

To summarize, here are the essential calculations:

• Factor $\mathbf K = \mathbf L \mathbf L^T$, using a sparse direct package.
• Form $\mathbf Y = \mathbf L^{-1} \mathbf Z$, basically a "fat" backsolve across multiple RHS's. But note that this is only "half" of the normal sparse forward/backward solution. A high quality sparse solver package will expose this as an independent calculation (MKL PARDISO does, and so does my own package, MyraMath).
• Initialize the small matrix $\mathbf S = \mathbf Y^T \mathbf Y$. This is basically a "fat" dot product, can be done with BLAS [dsyrk].
• Downdate $\mathbf S = \mathbf S -\frac{1}{\sigma} \mathbf I$.
• Form the $\mathbf L \mathbf D \mathbf L^T$ decomposition of $\mathbf S$. Can be done with LAPACK [dsytrf].
• Form the $\mathbf Q \mathbf R$ decomposition of $\mathbf Y$. Can be done with LAPACK [dgeqrf].
• Form the small matrix $\mathbf T = \mathbf R \mathbf S^{-1} \mathbf R^{T}$. You can use [dystrs] to apply $\mathbf S^{-1}$ to $\mathbf R$ into a temporary, then use [dgemm] to form $\mathbf R \cdot \left( \mathbf S^{-1} \mathbf R \right)$
• Compute the eigendecomposition of $\mathbf T = \mathbf V \mathbf \Lambda \mathbf V^T$, using LAPACK's [dsyevd].
• The number of negative eigenvalues of $\mathbf K - \sigma \mathbf Z \mathbf Z^T$ is the number of $\mathbf \Lambda_i$'s above 1. All the others eigenvalues are positive.

... I think? It seems reasonably efficient to me, as long as you are able to explicitly store $\mathbf Z$ (and therefore $\mathbf Y$). As an aside, note that $\mathbf N$ is never needed/computed, it's only introduced for explanatory purposes.

• Thank you very much for this great answer. And yes, it is reasonably efficient of a procedure. So, although $L^{-T} [U N]$ is not lower triangular, the small-dimension matrix in the middle contains the inertia? I am going to implement and see for myself. In addition, I have the feeling that the small eigendecomposition might be unnecessary, because an LDL factorization would be enough instead? Let me go through the algebra. – Olivier Aug 21 '18 at 18:33
• Yes, my claim about equal inertia is all drawing from: en.wikipedia.org/wiki/Matrix_congruence – rchilton1980 Aug 21 '18 at 18:39
• Highly encourage you to try it out, hope it works. For small problem sizes you can probably implement this in a few lines of Matlab, and perhaps even ignore the sparsity of $\mathbf K$. Not efficient, but it could be a fast way to find a random counterexample, in the event there is an error in the procedure. – rchilton1980 Aug 21 '18 at 18:42
• Yes, that's what I was planning on doing. I will let you know. – Olivier Aug 21 '18 at 18:45
• I have lost the thread here a bit, would you mind writing your own answer that summarizes the best/streamlined procedure? Then feel free to accept it.. you've certainly made the good faith effort to become the expert on your problem. Overall I think future visitors would be better served if they didn't have to dig so deeply into the commentary. – rchilton1980 Aug 24 '18 at 12:29

Let $K = L L^T$ be the Cholesky factorization of $K$.

Let $Y = L^{-1} Z$ , computed by forward substitution.

Let $\mathcal{I} (A)$ denote the number of negative eigenvalues of a given matrix $A$ (negative index of inertia).

We seek to compute $\mathcal{I}_\sigma = \mathcal{I} (K - \sigma M) = \mathcal{I} (L L^T - \sigma Z Z^T)$ .

Since the inertia is invariant under a change of basis (see https://en.wikipedia.org/wiki/Matrix_congruence and https://en.wikipedia.org/wiki/Sylvester%27s_law_of_inertia), inertia count $\mathcal{I}_\sigma$ can be written as

$\mathcal{I}_\sigma = \mathcal{I} \left( L^{-1} \left( L L^T - \sigma Z Z^T \right) L^{-T} \right) = \mathcal{I} \left( I_N - \sigma Y Y^T \right) = \mathcal{I} \left( \frac{1}{\sigma} I_N - Y Y^T \right)$ .

Now, it should be noted that $Y Y^T$ and $Y^T Y$ have the same nonzero eigenvalues (see for instance https://en.wikipedia.org/wiki/Singular-value_decomposition). The eigendecomposition

$\left( Y Y^T \right) \Phi_Y = \Phi_Y \begin{bmatrix} \Lambda_Y & 0\\ 0 & 0 \end{bmatrix}$ such that $Y Y^T = \Phi_Y \begin{bmatrix} \Lambda_Y & 0\\ 0 & 0 \end{bmatrix} \Phi_Y^T$ and $\Phi_Y \Phi_Y^T = I_N$ yields

$\mathcal{I}_\sigma = \mathcal{I} \left( \Phi_Y \left( \frac{1}{\sigma} I_N - \begin{bmatrix} \Lambda_Y & 0\\ 0 & 0 \end{bmatrix} \right) \Phi_Y^T \right)$ and consequently,

$\mathcal{I}_\sigma = \mathcal{I} \left( \frac{1}{\sigma} I_N - \begin{bmatrix} \Lambda_Y & 0\\ 0 & 0 \end{bmatrix} \right)$ , which is given by $\mathcal{I}_\sigma = \mathcal{I} \left( \frac{1}{\sigma} I_r - \Lambda_Y \right)$ .

Matrix $\Lambda_Y$ is the diagonal matrix of the eigenvalues of $Y^T Y$. Once $\Lambda_Y$ is calculated (independent of $\sigma$), then $\mathcal{I}_\sigma$ can be computed at no cost for any $\sigma > 0$.

Alternatively, we now show how it is possible to compute all the $r$ finite eigenvalues and associated eigenvectors of the generalized eigenvalue problem (GEP) $K \Phi = M \Phi \Lambda$ , with the same cost.

Recalling that $M = Z Z^T$, it can be deduced that the eigenvectors belong to the image (column space) of $(N \times r)$ real matrix $S = K^{-1} Z$ (see for instance Large-scale generalized eigenvalue problem with low rank LHS matrix).

Applying Rayleigh-Ritz (see https://en.wikipedia.org/wiki/Rayleigh%E2%80%93Ritz_method) to above GEP (Galerkin projection onto $S$), a $(r \times r)$ GEP is obtained,

$(S^T K S) X = (S^T M S) X \Lambda$.

That is to say, an eigenvector $\boldsymbol{\varphi}$ of the $(N \times N)$ GEP is written as $\boldsymbol{\varphi} = S \textbf{x}$ with $\textbf{x}$ an eigenvector of the $(r \times r)$ GEP (that is, $\Phi = S X$). This procedure yields exact eigenpairs since $S$ spans the sought eigenspace.

Now, it can be observed that $S^T K S = Y^T Y$ and $S^T M S = {(Y^T Y)}^2$. Introducing the notation $A = Y^T Y$, the $(r \times r)$ GEP writes $A X = A^2 X \Lambda$, which is equivalent to the standard eigenvalue problem (SEP) $A X = X \Lambda^{-1}$ . It can be recognized that $\Lambda = \Lambda_Y^{-1}$.

As a conclusion, the proposed approach for inertia count is as costly as directly computing all the eigenvalues of the associated GEP. The cost (dense eigenvalue problem with dimension $r$) rapidly increases with $r$.

It should be noted that the cost for computing $\mathcal{I}_\sigma$ for one given $\sigma > 0$ can be less than that of computing all the eigenvalues, as one also has $\mathcal{I}_\sigma = \mathcal{I} \left( \frac{1}{\sigma} I_r - Y^T Y \right)$, which can be obtained by LDL factorization of $\frac{1}{\sigma} I_r - Y^T Y$.

My next question would be how to determine the inertia count of $K_\sigma - \sigma Z Z^T$ with $K_\sigma = K - \sigma M_0$ , in which $M_0$ is a sparse $(N \times N)$ symmetric positive-semidefinite matrix. The difference with the above is that $K_\sigma$ is not positive definite. I am presently using the above developments to answer this question.

Edit: See below for the answer to this question.

We seek to compute $\mathcal{I}_\sigma = \mathcal{I} \left( K - \sigma M \right)$ with $M = M_0 + Z Z^T$.

We have $\mathcal{I}_\sigma = \mathcal{I} \left( K_\sigma - \sigma Z Z^T \right)$ with $K_\sigma = K - \sigma M_0$ a full-rank sparse symmetric matrix.

The LDL factorization of $K_\sigma$ is performed, which yields $K_\sigma = L_\sigma D_\sigma L_\sigma^T$.

Thus, we obtain $\mathcal{I}_\sigma = \mathcal{I} \left( L_\sigma^{-1} \left( K_\sigma - \sigma Z Z^T \right) L_\sigma^{-T} \right) = \mathcal{I} \left( D_\sigma - \sigma Y_\sigma Y_\sigma^T \right)$ with $Y_\sigma = L_\sigma^{-1} Z$. Which gives $\mathcal{I}_\sigma = \left( \frac{D_\sigma}{\sigma} - Y_\sigma Y_\sigma^T \right)$.

Introducing the $(N \times N)$ GEP $\left( \frac{D_\sigma}{\sigma} \right) \Phi = \left( Y_\sigma Y_\sigma^T \right) \Phi \Lambda$ such that $\Phi^T \left( \frac{D_\sigma}{\sigma} \right) \Phi = \Lambda$ and $\Phi^T \left( Y_\sigma Y_\sigma^T \right) \Phi = I_N$ (unit mass eigenvector normalization), we obtain $\mathcal{I}_\sigma = \mathcal{I} \left( \Phi^T \left( \frac{D_\sigma}{\sigma} - Y_\sigma Y_\sigma^T \right) \Phi \right) = \mathcal{I} \left( \Lambda - I_N \right)$. This GEP cannot be solved.

Similarly to before, the $r$ eigenvectors associated with the $r$ finite eigenvalues belong to the column space of $S_\sigma = {\left( \frac{D\sigma}{\sigma} \right)}^{-1} Y_\sigma$. The $(r \times r)$ GEP $\left( S_\sigma^T \left( \frac{D_\sigma}{\sigma} \right) S_\sigma \right) X = \left( S_\sigma^T \left( Y_\sigma Y_\sigma^T \right) S_\sigma \right) X \Lambda_f$ allows to obtain the finite eigenvalues $\Lambda_f$ of the $(N \times N)$ GEP. Let $A_\sigma = S_\sigma^T \left( \frac{D_\sigma}{\sigma} \right) S_\sigma$. One has $A_\sigma = Y_\sigma^T {\left( \frac{D_\sigma}{\sigma} \right)}^{-1} Y_\sigma$ and $S_\sigma^T \left( Y_\sigma Y_\sigma^T \right) S_\sigma = A_\sigma^2$. Thus, the $r$ finite eigenvalues $\Lambda_f$ can be deduced from the $(r \times r)$ SEP $A_\sigma X = X \Lambda_f^{-1}$.

Let $\Lambda_i$ denote the matrix of the $N - r$ infinite eigenvalues. We have $\mathcal{I}_\sigma = \mathcal{I} \left( \Lambda - I_N \right) = \mathcal{I} \left( \Lambda_f - I_r \right) + \mathcal{I} \left( \Lambda_i - I_{N - r} \right)$. Let $\mathcal{I}_f = \mathcal{I} \left( \Lambda_f - I_r \right)$ and $\mathcal{I}_i = \mathcal{I} \left( \Lambda_i - I_{N - r} \right)$. We have $\mathcal{I}_f = \mathcal{I} \left( X^T A_\sigma^{-1} X - X^T X \right)$ because $X^{-1} = X^T$, and consequently, $\mathcal{I}_f = \mathcal{I} \left( A_\sigma^{-1} - I_r \right)$. Since for a real scalar $a$ we have $a^{-1} \in ] - \infty , 1 [ \,\,\, \Rightarrow \,\, a \in ] - \infty , 0 [ \, \cup \,] 1 , + \infty [\,$, it can be seen that $\mathcal{I}_f = \mathcal{I} \left( A_\sigma^{-1} - I_r \right) = \mathcal{I} \left( A_\sigma \right) + \mathcal{I} \left( I_r - A_\sigma \right)$.

On the other hand, we have $\mathcal{I}_i = \mathcal{I} \left( \Lambda_i - I_{N - r} \right) = \mathcal{I} \left( \Lambda_i \right)$ ($\mathcal{I}_i$ is given by the number of negative infinite eigenvalues). Since $\mathcal{I} \left( \Lambda \right) = \mathcal{I} \left( \Lambda_f \right) + \mathcal{I} \left( \Lambda_i \right)$ and $\mathcal{I} \left( \Lambda \right) = \mathcal{I} \left( \Phi^T \left( \frac{D_\sigma}{\sigma} \right) \Phi \right) = \mathcal{I} \left( \frac{D_\sigma}{\sigma} \right) = \mathcal{I} \left( D_\sigma \right)$, we obtain $\mathcal{I}_i = \mathcal{I} \left( D_\sigma \right) - \mathcal{I} \left( \Lambda_f \right)$. In addition, we have $\mathcal{I} \left( \Lambda_f \right) = \mathcal{I} \left( \Lambda_f^{-1} \right) = \mathcal{I} \left( A_\sigma \right)$. Therefore, we obtain $\mathcal{I}_i = \mathcal{I} \left( D_\sigma \right) - \mathcal{I} \left( A_\sigma \right)$.

As a conclusion, it can be deduced that $\mathcal{I}_\sigma = \mathcal{I} \left( D_\sigma \right) + \mathcal{I} \left( I_r - A_\sigma \right)$.

Summary of the operations:

• LDL factorization of $(N \times N)$ sparse matrix $K_\sigma = L_\sigma D_\sigma L_\sigma^T$

• Solve $Y_\sigma = L_\sigma^{-1} Z$ by forward substitutions

• Assemble $A_\sigma = Y_\sigma^T {\left( \frac{D_\sigma}{\sigma} \right)}^{-1} Y_\sigma$

• LDL factorization of $(r \times r)$ dense matrix $I_r - A_\sigma$