4
$\begingroup$

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?

$\endgroup$
5
  • $\begingroup$ Does this answer your question? Cheap recalculation of eigenvalues and eigenvectors for a low-rank update of the matrix $\endgroup$ Jan 31, 2023 at 0:55
  • $\begingroup$ @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 $\endgroup$
    – eepperly16
    Jan 31, 2023 at 1:03
  • 1
    $\begingroup$ 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. $\endgroup$ Jan 31, 2023 at 1:19
  • $\begingroup$ 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. $\endgroup$ Jan 31, 2023 at 1:25
  • $\begingroup$ And, it looks like dlaed9() will do all the work to update the eigenvalues and Q. $\endgroup$ Jan 31, 2023 at 1:29

1 Answer 1

6
$\begingroup$

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:

enter image description here

$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$ Sep 15, 2023 at 8:54
  • $\begingroup$ @MarcelFerrari I fully agree. Looking forward to the beta $\endgroup$
    – eepperly16
    Sep 16, 2023 at 15:23
  • 1
    $\begingroup$ 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. $\endgroup$ Sep 20, 2023 at 15:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.