Efficient eigendecomposition of DPR1 matrices

Question

I'm finding that the following bit is the bottleneck in my code

d = 2000;
h = 1./Range[d];
evecs = DiagonalMatrix[h] + {h}\[Transpose] . {h} // 
    Eigenvectors; // Timing  (* 3.80948 *)

Matrix in question is known as DPR1 (diagonal + rank1) matrix for which efficient algorithms exist. In my case it's also symmetric + positive definite.

There seems to be a large literature on DPR1 matrices and Matlab / Julia code. I'm looking for something comparable in Mathematica.

Answer 1

Okay, apparently, the root finding for the characteristic polynomial is quite easy. Since the derivative of the characteristic polynomial is easy to compute, too, I try to use Newton's method where possible.

This is the function; I designed it to work for real symmetric, positive-definite matrices. I don't know whether it would work for indefinite ones, too.

This has grown in a more complex project than I thought. Thus for maintenance reasons, I decided to create the following public github repository:

https://github.com/HenrikSchumacher/DPR1Eigensystem

And here is a usage example:

Needs["DPR1Eigensystem`"];

n = 10000;
diag = RandomReal[{1, n}, n];
z = RandomReal[{-1, 1}, n];

{\[Lambda]True, UTrue} = Eigensystem[DiagonalMatrix[diag] + KroneckerProduct[z, z]]; // AbsoluteTiming // First

{\[Lambda], U} =  DPR1Eigensystem[diag, z, MaxIterations -> 20,  Tolerance -> 10^-14]; // AbsoluteTiming // First

Max[Abs[\[Lambda]True - \[Lambda]]]/Max[Abs[\[Lambda]True]]
Max[DPR1TestEigensystem[diag, z, \[Lambda], U]]

51.0864

0.239854

1.26254*10^-14

8.65564*10^-9

More than two hundred times faster than the naive method and apparently quite accurate. I am a bit disappointed by the low accuracy of the eigenvectors, though. I guess one can get a bit better by refining them with a few Raleigh iterations...

It even works decently for matrices so large that Eigensystem fails:

n = 30000;
diag = RandomReal[{1, n}, n];
z = RandomReal[{-1, 1}, n];

{\[Lambda], U} =  DPR1Eigensystem[diag, z, MaxIterations -> 20, Tolerance -> 10^-14]; // AbsoluteTiming // First
Max[DPR1TestEigensystem[diag, z, \[Lambda], U]]

2.03043

1.27448*10^-8

Disclaimer

All this is done under the assumption that no two elements of diag coincide. The Cauchy interlacing theorem implies that no eigenvalue has a multiplicity. Not sure how to deal with multiplicities...

Answer 2

The most important facts about DPR1 matrices of size $n \times n$ are:

1. Their matrix-vector multiplication can be computed in $O(n)$ time instead of $O(n^2)$ as for ordinary dense matrices.

2. By virtue of the Sherman-Morrison formular, also the matrix-vector multiplication with their (shifted) inverse can be computed in $O(n)$ time instead of $O(n^3)$ as for ordinary dense matrices.

The routines DPR1Apply and DPR1ApplyShiftedInverse from the code section at the bottom of this post facilitate this.

Here is a brief usage example that demonstates that the inverse of the shifted matrix can be accurately evaluated several orders of magnitude faster than the naive way:

n = 6000;

diag = Sort@RandomReal[{1., 10.}, n];
z = RandomReal[{-1, 1}, n];
A = DiagonalMatrix[diag] + KroneckerProduct[z, z];

(* We need to detect whether the shifted diagonal diag-\[Mu] has a \
zero entry. *)
(* We use a NearestFunction to speed this up. In practice one should \
reuse it as often as possible.*)
nf = Nearest[diag -> "Index"];

(*some shift*)
\[Mu] = RandomReal[{1, 10.}];
(*some vector*)
u = RandomReal[{-1, 1}, n];

vTrue = LinearSolve[A - \[Mu] IdentityMatrix[n], u]; // 
  AbsoluteTiming // First
v = DPR1ApplyShiftedInverse[diag, z, \[Mu], u, nf]; // 
  AbsoluteTiming // First
Norm[vTrue - v]/Norm[vTrue]

0.6447

0.000151

1.47582*10^-12

Hence every iterative eigenvalue solver that uses only the action of the matrix and its (shifted) inverse can be speed up by incorporating fast routines for their actions. For computing only few eigenvalues of a positive-definite DPR1 matrix, one could use, for example, the Arnoldi's method.

I don't have a editable implementation of Arnoldi's method at hand. Thus I demonstrate this with the easy to implement Raleigh iterations instead. It is a cubically convergent variant of the inverse power iteration (which converges only linearly). With an initial guess \[Mu]0 for the eigenvalue to be approximated, the algorithm can be implemented like in the function DPR1RaleighIterations in the code section at the bottom. It is certainly not the most efficient way to do so, though.

For symmetric, positive-definite input, it is very effective in determining the extreme eigenvalues of A.

Here the preparations:

n = 10000;

diag = Sort@RandomReal[{1., 10.}, n];
z = RandomReal[{-1, 1}, n];
A = DiagonalMatrix[diag] + KroneckerProduct[z, z];

nf = Nearest[diag -> "Index"];

For the smallest eigenvalue, we can simply take the smallest element on the diagonal:

\[Mu]SmallestGuess = Min[diag];
\[Mu]True = 
    Eigenvalues[A, {-1}, 
      Method -> {"Arnoldi", "Shift" -> \[Mu]SmallestGuess}][[1]]; // 
  AbsoluteTiming // First
\[Mu] = DPR1RaleighIterations[diag, z, \[Mu]SmallestGuess, nf]; // 
  AbsoluteTiming // First
"Relative Error" -> Abs[\[Mu]True - \[Mu]]/Abs[\[Mu]True]

3.9683

0.016788

"Relative Error" -> 5.77044*10^-15

This is over 200 times faster than Arnoldi's method and very, very accurate. It could be made faster by removing some redundant vector divisions in the code. Compiling the the Raleigh iterations and replacing the NearestFunction by reordering diag and employing a binary search might help, too.

Also the largest eigenvalue can be computed quickly and reliably:

\[Mu]LargestGuess = Max[diag] + z . z;
\[Mu]True = 
    Eigenvalues[A, {1}, 
      Method -> {"Arnoldi", "Shift" -> \[Mu]LargestGuess}][[1]]; // 
  AbsoluteTiming // First
\[Mu] = DPR1RaleighIterations[diag, z, \[Mu]LargestGuess, nf]; // 
  AbsoluteTiming // First
"Relative Error" -> Abs[\[Mu]True - \[Mu]]/Abs[\[Mu]True]

3.45909

0.01716

"Relative Error" -> 2.15118*10^-15

Note that I use some cheap upper bound for the spectral radius as ininital guess. The largest diagonal element is typically a good initial guess for the second largest eigenvalue, not for the largest.

Intermediate eigenvalues are an issue. There is no guarantee that Raleigh iterations using RankedMin[diag, k] as initial guess will converge to the $k$-th smallest eigenvalue. But by Cauchy's interlace theorem, it should be almost always converge towards RankedMin[diag, k-1] or RankedMin[diag, k]. The following experiments seem to back this up:

Precompute the eigenvalues once:

\[CapitalLambda]True = Eigenvalues[A];

Run the following as often as you want:

k = RandomInteger[{1, n}];
\[Mu]IntermediateGuess = RankedMin[diag, k];
\[Mu] = DPR1RaleighIterations[diag, z, \[Mu]IntermediateGuess, nf];

"Relative Error k-th" -> 
 Abs[RankedMin[\[CapitalLambda]True, k] - \[Mu]] / Abs[RankedMin[\[CapitalLambda]True, k]]

"Relative Error (k-1)-st" -> Min[{
   If[k > 0, 
    Abs[RankedMin[\[CapitalLambda]True, k - 1] - \[Mu]] / Abs[RankedMin[\[CapitalLambda]True, k - 1]], Nothing]
   }]

Always, either "Relative Error k-th" or "Relative Error (k-1)-st" will be zero.

To use this to predict the full spectrum without omitting any eigenvalues is a bit more involved. But since one can say quite precisely where the eigenvalues have to lie (again by Cauchy's interlacing theorem), one can use interval bisection to finally find all eigenvalues.

Finally, the eigenvectors ccan be computed from the eigenvalues with the function DPR1EigenvectorFromEigenvalue below.

Usage example (assuming that \[Mu] is an eigenvalue):

v = DPR1EigenvectorFromEigenvalue[diag, z, \[Mu]];
Norm[DPR1Apply[diag, z, v] - \[Mu] v]/Norm[v]

2.94404*10^-9

Not super-duper accurate, but maybe it suffices for your application.

Literature

I was able to derive this quite quickly from the well-written introduction of this paper by Stor et. al.. Instead of an interative solver, the authors try to find the smallest root of the characteristic polynomial for the shifted matrix. They can do it by interval bisection, and it is probably way faster than Raleigh iterations. (The point is that the characteristic polynomial can be computed in $O(n)$ time without involving any determinants.) But I have to wrap my head around that on another day.

Complete code

I got fed up with copy-pasting all the updates to the code here. Thus I simply created the following public github repository:

https://github.com/HenrikSchumacher/DPR1Eigensystem

Answer 3

First attempt. We can use Eigenvalues[] instead of Eigenvectors[] to compute the eigenvectors in a case of DPR1 as follows

dpr10[d_] := 
  Module[{h, mat, dm, id, l, w, x}, h = 1./Range[d]; 
   id = IdentityMatrix[d]; dm = DiagonalMatrix[h];
   mat = (dm + {h}\[Transpose] . {h}); 
   l = Eigenvalues[mat, Method -> "Direct"]; 
   w = Table[x = h/(h - l[[i]]); Normalize[x], {i, d}]; w];

This function we compare to

vec[d_] := Module[{h, mat, e}, h = 1./Range[d];
  mat = (DiagonalMatrix[h] + {h}\[Transpose] . {h}); 
  e = Eigenvectors[mat]; e]

Note, that vec is actually code to be improved. We have for d=10^3 practically same results

dpr10[10^3]; // AbsoluteTiming

Out[]= {0.0751313, Null}

vec[10^3]; // AbsoluteTiming

Out[]= {0.0960007, Null}

For d=10^4 function dpr10 looks better

dpr10[10^4]; // AbsoluteTiming
{14.0283, Null}

vec[10^4]; // AbsoluteTiming
{76.4033, Null}

But for d=2 10^4 function vec10 looks much better

dpr10[2 10^4];//AbsoluteTiming
 {83.6977,Null}

 vec[2 10^4];//AbsoluteTiming
 {588.388,Null}