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

  • 1
    $\begingroup$ Even if $N$ is huge, can you factorize $K$? $\endgroup$ – user7440 Aug 20 '18 at 3:02
  • $\begingroup$ 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'. $\endgroup$ – Olivier Aug 21 '18 at 4:58
  • 1
    $\begingroup$ 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. $\endgroup$ – rchilton1980 Aug 21 '18 at 12:53
  • $\begingroup$ 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. $\endgroup$ – 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.

  • 1
    $\begingroup$ 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. $\endgroup$ – Olivier Aug 21 '18 at 18:33
  • $\begingroup$ Yes, my claim about equal inertia is all drawing from: en.wikipedia.org/wiki/Matrix_congruence $\endgroup$ – rchilton1980 Aug 21 '18 at 18:39
  • $\begingroup$ 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. $\endgroup$ – rchilton1980 Aug 21 '18 at 18:42
  • $\begingroup$ Yes, that's what I was planning on doing. I will let you know. $\endgroup$ – Olivier Aug 21 '18 at 18:45
  • 1
    $\begingroup$ 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. $\endgroup$ – rchilton1980 Aug 24 '18 at 12:29

This answer is based on rchilton1980 answer and the comments that have followed.

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$


Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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