Eigenvalues of diagonal plus rank-one

Question

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?

Answer 1

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: