# 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 Jan 31 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 Jan 31 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. Jan 31 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. Jan 31 at 1:25
• And, it looks like dlaed9() will do all the work to update the eigenvalues and Q. Jan 31 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);
v = mxGetDoubles(prhs);
c = mxGetDoubles(prhs);
#else
d = mxGetPr(prhs);
v = mxGetPr(prhs);
c = mxGetPr(prhs);
#endif
n = mxGetM(prhs);

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

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

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

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

/* create output matrix C */
plhs = mxCreateDoubleMatrix(n, n, mxREAL);
plhs = mxCreateDoubleMatrix(n, 1, mxREAL);;
#if MX_HAS_INTERLEAVED_COMPLEX
Q = mxGetDoubles(plhs);
l = mxGetDoubles(plhs);
#else
Q = mxGetPr(plhs);
l = mxGetPr(plhs);
#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) != n || mxGetN(prhs) != 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: 