# Eigenvalues of diagonal plus rank-one

I need to compute an eigendecomposition of an $$n\times n$$ matrix

$$D + c vv^\top = Q\Lambda Q^\top \tag{1}$$

in MATLAB, where $$D$$ is a real diagonal matrix, $$c$$ is a scalar, and $$v$$ is a real vector. As is detailed in section 5.3.3 of Demmel's Applied Numerical Linear Algebra, this can be achieved in $$\mathcal{O}(n^2)$$ operations. As discussed in this book, the numerical implementation of this algorithm is quite subtle, and I would like to leverage a high-quality implementation such as in LAPACK. Unfortunately, even if I am willing to call LAPACK from MATLAB with a MEX file, I am not whether the LAPACK subroutines are written in such a way as to isolate the $$\mathcal{O}(n^2)$$ computation of (1) that I need. For instance, dlaed1 appears to run in $$\mathcal{O}(n^3)$$ operation because it is integrated with steps in the divide and conquer eigensolver. This brings me to my question:

Is there a high-quality piece of software for computing (1) in $$\mathcal{O}(n^2)$$ operations? Can this be called in MATLAB?

• Does this answer your question? Cheap recalculation of eigenvalues and eigenvectors for a low-rank update of the matrix Commented Jan 31, 2023 at 0:55
• @BrianBorchers This does not answer my question. The linked article is in the case where $D$ is a dense matrix whose eigendecomposition is known; in my case, $D$ is just a diagonal matrix. This makes the $\mathcal{O}(n^2)$ algorithm suggested in the linked question lead to a speedup for my problem, whereas in that problem the speedup is destroyed by needing to compute the product of two $n\times n$ eigenmatrices. In addition, my question concerns software implementation of the $\mathcal{O}(n^2)$ algorithm, not discussed in the linked question Commented Jan 31, 2023 at 1:03
• The point is that your question should be closed and someone should answer the original question by pointing to Demmel's Algorithm 5.3 and discussing how it might be implemented. Looking at it, I see no particular reason to need LAPACK in any of the steps of that algorithm. Commented Jan 31, 2023 at 1:19
• Actually, dlaed4() could be used to solve the secular equation to get the eigenvalues of the updated matrix. The remaining two steps of Algorithm 5.3 could then be done with your own code. Commented Jan 31, 2023 at 1:25
• And, it looks like dlaed9() will do all the work to update the eigenvalues and Q. Commented Jan 31, 2023 at 1:29

Thanks to Brian Borchers for suggesting dlaed9; this mostly solves my problem.

Unfortunately, dlaed9 does not include deflation if $$D$$ has repeated entries or $$v$$ has zero entries, in which case dlaed9 will either produce an error message or populate the output eigenvalues/eigenvectors.

To help anyone who may follow down the same path I did, here is a MEX file which calls dlaed9. I have not made an effort to implement deflation, but I have written checks for errors so this issue should at least not fail silently.

MATLAB file update_eig.m

function [Q,D] = update_eig(d,v,c)
%% Negate, sort, and normalize
negative = c < 0;
if negative; d = -d; c = -c; end
[d,p] = sort(d,'ascend');
v = v(p);
c = c * norm(v)^2;
v = v / norm(v);

%% Compute update
[Q,D] = update_eig_mex(d,v,c);

%% Undo negation and sorting
Q(p,:) = Q;
if negative; D = -D; end
D = diag(D);
end


MEX file update_eig_mex.c

// Compute eigendecomposition of diag(d) + c*v*v', for a unit vector v

#include "mex.h"
#include "lapack.h"
#include "math.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
#if MX_HAS_INTERLEAVED_COMPLEX
mxDouble *d, *v, *c, *Q, *l;
#else
double *d, *v, *c, *Q, *l;
#endif
ptrdiff_t n;      /* matrix dimensions */

#if MX_HAS_INTERLEAVED_COMPLEX
d = mxGetDoubles(prhs[0]);
v = mxGetDoubles(prhs[1]);
c = mxGetDoubles(prhs[2]);
#else
d = mxGetPr(prhs[0]);
v = mxGetPr(prhs[1]);
c = mxGetPr(prhs[2]);
#endif
n = mxGetM(prhs[0]);

for (int i = 0; i < 3; ++i) {
if( !mxIsDouble(prhs[0]) || mxIsComplex(prhs[0])) {
mexErrMsgIdAndTxt( "MATLAB:update_eig:fieldNotRealMatrix",
"Input arguments must be a real matrix.");
}
}

if (mxGetN(prhs[0]) != 1) {
mexErrMsgIdAndTxt("MATLAB:update_eig:dims",
"First input must be column vector");
}

if (mxGetM(prhs[1]) != n || mxGetN(prhs[1]) != 1) {
mexErrMsgIdAndTxt("MATLAB:update_eig:dims",
"Second input must be column vector, same size as first");
}

if (mxGetM(prhs[2]) != 1 || mxGetN(prhs[1]) != 1) {
mexErrMsgIdAndTxt("MATLAB:update_eig:dims",
"Third input must be scalar");
}

/* create output matrix C */
plhs[0] = mxCreateDoubleMatrix(n, n, mxREAL);
plhs[1] = mxCreateDoubleMatrix(n, 1, mxREAL);;
#if MX_HAS_INTERLEAVED_COMPLEX
Q = mxGetDoubles(plhs[0]);
l = mxGetDoubles(plhs[1]);
#else
Q = mxGetPr(plhs[0]);
l = mxGetPr(plhs[1]);
#endif

double *work = malloc(sizeof(double) * n*n);

ptrdiff_t info;
ptrdiff_t one = 1;

/* Pass arguments to Fortran by reference */
dlaed9(&n, &one, &n, &n, l, work, &n, c, d, v, Q, &n, &info);

if (info != 0) {
if (mxGetM(prhs[1]) != n || mxGetN(prhs[1]) != 1) {
mexWarnMsgIdAndTxt("MATLAB:update_eig:dslaed9_error",
"dlaed9 exited with error %d", info);
}
}

bool found_nans = false;
for (int i = 0; i < n; ++i){
if (isnan(l[i])) {
mexWarnMsgIdAndTxt("MATLAB:update_eig:nans",
"d vector has NaNs, problem not properly deflated");
break;
}
for (int j = 0; j < n; ++j) {
if (isnan(Q[i + n*j])) {
mexWarnMsgIdAndTxt("MATLAB:update_eig:nans",
"Q matrix has NaNs, problem not properly deflated");
found_nans = true;
break;
}
}
if (found_nans) break;
}

free(work);
}


Confirming $$\mathcal{O}(n^2)$$ runtime:

• dlaed9 has its limitations as it is meant for a very specific purpose. I am currently working on a C++ based suite of eigensolvers which, among other special matrices, supports the DPR1 case. Once a beta version is ready I will add a comment with a link to that. Commented Sep 15, 2023 at 8:54
• @MarcelFerrari I fully agree. Looking forward to the beta Commented Sep 16, 2023 at 15:23
• The code is still super messy and early in development, but it runs pretty fast and generates correct results now github.com/MarcelFerrari/Freccia/blob/main/src/DPR1/… I've done 10kx10k DPR1 matrix in ~800 ms using MKL backend and running in parallel. Commented Sep 20, 2023 at 15:52