Suggestions for speeding up this eigen code (woodbury formula for a very sparse low-rank update)

edit- i am compiling with the following flags: g++ -march=native -O3

I am writing code to compute the inverse of a matrix $A+UU^T$ with inverse $A^{-1}$ already computed. here U,V are both very sparse matrices. A is NxN, U is Nxk, V=U^T=Nxk where k~=N/100, and V,U have only k non-zero entries

Here’s the comparative speed for vanilla matrix inversion of $A+UV$, then with the woodbury formula, and then using the sparse matrix construction for U,V.

Are there any clear routes to speed up the calculation further? one thing to note is that A_inv is a covariance (and so symmetric) matrix, and UU^T only has diagonal entries

vanilla calculation ms: 1451.49
loop (10) woodbury: 66.3399
loop (10) woodbury sparse: 51.3955

Eigen::MatrixXd woodburyInverseF(const Eigen::MatrixXd& A_inv,
                                 const Eigen::MatrixXd& identity_Matrix,
                                 const Eigen::SparseMatrix<double>& U,
                                 const Eigen::SparseMatrix<double>& V,
                                 const Eigen::MatrixXd& C_inv) {
    auto start = std::chrono::high_resolution_clock::now();
    auto start_middle = std::chrono::high_resolution_clock::now();

    Eigen::MatrixXd middle = (identity_Matrix + (V * A_inv) * U).inverse();
    Eigen::MatrixXd vainv = V * A_inv;
    Eigen::MatrixXd ainvu = A_inv * U;
    Eigen::MatrixXd wood_A_inv = A_inv - ainvu *(middle *vainv);
//    Eigen::MatrixXd wood_A_inv = A_inv - (A_inv * U) * ((identity_Matrix + (V * A_inv) * U).inverse()) * (V * A_inv);

    return wood_A_inv;

I have tried forcing order of operations to make sure that the terms middle, vainv, ainvu are computed first to keep all the intermediate steps of the calculation as fast as possible (naiive eigen was almost as bad as the original).

I expect there is further improvement on construction as this is day 1 of me using eigen, so I’m not too familiar with the available speedups / optimisation.

  • a) don’t use mathjax, it’s not enabled on stackoverflow b) please post a minimal reproducible example along with input that mimics matrices you use so we can run it

    – 




Leave a Comment