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: