[eigen] Feature proposal: fast updates to QR factorizations

[Patrick R Conrad <prconrad@xxxxxxx>]

Hi everyone,
I've been using Eigen happily for a while now, and I've started developing a feature that might be a useful contribution to the library. I'd like to assess whether there's any interest and to ask for some guidance.

In my application, we want to quickly re-solve a least squares problem with deleted rows, providing leave-one-out cross validation. If we solve the least squares system with QR, this amounts to a fast update of the factorization. Someone actually asked for this feature at least once:
https://forum.kde.org/viewtopic.php?f=74&t=95534&p=198193&hilit=update+qr#p198193

Apparently it's not all that hard to do. This paper outlines suitable algorithms for updating a QR for adding/deleting rows/columns:
http://eprints.ma.man.ac.uk/1192/01/covered/MIMS_ep2008_111.pdf
The process involves additional householder or givens rotations on a carefully chosen set of entries. Since Eigen already has these operations, it's not a ton of work to apply them.

For reasonably conditioned matrices, this is fairly accurate and can be much faster, especially if you have to do a lot of re-solves. All I need is the relatively specific case of removing one row and resolving, which I believe I have implemented successfully. I've attached a sample implementation that shows how to do this update for HouseholderQR and ColPivHouseholderQR.

A sample output looks like this:
*******************************
>> g++ fast_qr_update_sample.cpp --std=c++11 -O3
>> ./a.out
Standard HouseholderQR: test matrix = 1001 x 1000
Complete refactoring time: 0.65119 seconds.
Fast update to QR to resolve: 0.111785 seconds.
Max coefficient error in solve, compared to resolve with ColPivQR: 1.6378e-12
Max coefficient error reconstructing the LHS with updated QR 4.55191e-15

ColPivHouseholderQR: test matrix = 1001 x 1000
Complete refactoring time: 0.558148 seconds.
Fast update to QR to resolve: 0.110247 seconds.
Max coefficient error in solve, compared to resolve with ColPivQR: 1.68043e-12
*******************************

Note that C++11 is only needed for timing features from chrono, nothing else depends upon it. Both version include re-solving a least squares for a given RHS without forming  the revised Q, and for HouseholderQR there's a method that explicitly forms the new Q and R. These results on random matrices look favorable to me.

This is obviously just a proof of concept, and requires much more work before it would be widely useful. If there's suitable interest, I'm willing to work towards that. Before I can proceed, however, I need some help figuring out how these features should work - a suitable interface, and how it fits into the existing QR system. I look forward to your thoughts.

Regards,
Patrick Conrad

#include <iostream>
#include <chrono>
#include <assert.h>

#include <Eigen/SVD>
#include <Eigen/QR>

using namespace Eigen;
using namespace std;

//Proposed low rank updates to QR methods for Eigen.
//Copyright Patrick R. Conrad 2014

MatrixXd RemoveARow(MatrixXd mat, int rowIndex)
{
  MatrixXd removed(mat.rows() - 1, mat.cols());

  removed << mat.topRows(rowIndex), mat.bottomRows(mat.rows() - rowIndex - 1);
  return removed;
}

//Resolve the least squares problem from a QR and a row to remove. The later methods are just variants on this one.
MatrixXd ReSolveWithRemovedRow(HouseholderQR<MatrixXd>& lhs_qr, MatrixXd const& rhs, int rowToRemove)
{
  MatrixXd r =  lhs_qr.matrixQR().template triangularView<Upper>();

  assert(r.rows() == rhs.rows());
  assert(rowToRemove < r.rows());
  assert(r.rows() > r.cols());

  MatrixXd newRhs = rhs;

  newRhs.applyOnTheLeft(lhs_qr.householderQ().transpose());


  //Get the row of Q we need by applying it to a unit vector
  VectorXd q_to_zero = VectorXd::Unit(rhs.rows(), rowToRemove);
  q_to_zero.applyOnTheLeft(lhs_qr.householderQ().transpose());


  //Work through q_to_zero, also applying the rotation to R and the rhs.
  JacobiRotation<double> rot;
  for (int i = q_to_zero.rows() - 1; i > 0; --i) {
    rot.makeGivens(q_to_zero(i - 1), q_to_zero(i));
    q_to_zero.applyOnTheLeft(i - 1, i, rot.adjoint()); //the actual vector defining the work we have to do
    newRhs.applyOnTheLeft(i - 1, i, rot.adjoint());    //rotate the rhs
    r.applyOnTheLeft(i - 1, i, rot.adjoint());         //work on R
  }


  r.block(1, 0, r.cols(),
          r.cols()).template triangularView<Upper>().solveInPlace(newRhs.block(1, 0, r.cols(), newRhs.cols()));
  return newRhs.block(1, 0, r.cols(), newRhs.cols());
}

//Take a Householder QR, and make a new QR (newq and newr) for the matrix with a row removed. Return Q*R so we can test
// its effectiveness.
MatrixXd ApproximateQRWithRemovedRow(HouseholderQR<MatrixXd>& lhs_qr, int rowToRemove)
{
  MatrixXd q_old =  lhs_qr.householderQ();
  MatrixXd r     =  lhs_qr.matrixQR().template triangularView<Upper>();

  assert(rowToRemove < r.rows());
  assert(r.rows() > r.cols());


  //Permute q to move the row we're deleting to the top
  MatrixXd q(q_old.rows(), q_old.cols());
  if (rowToRemove < q_old.rows() - 1) {
    q << q_old.row(rowToRemove), q_old.block(0, 0, rowToRemove, q_old.cols()), q_old.block(rowToRemove + 1,
                                                                                           0,
                                                                                           q_old.rows() - rowToRemove - 1,
                                                                                           q_old.cols());
  } else { //special case when accessing row = (rowToRemove+1) would throw an error
    q << q_old.row(rowToRemove), q_old.block(0, 0, rowToRemove, q_old.cols());
  }


  VectorXd q_to_zero = q_old.row(rowToRemove);

  //Each iteration, work on
  JacobiRotation<double> rot;                          //q_to_zero(i), q_to_zero(i+1));
  for (int i = q_to_zero.rows() - 1; i > 0; --i) {
    rot.makeGivens(q_to_zero(i - 1), q_to_zero(i));
    q_to_zero.applyOnTheLeft(i - 1, i, rot.adjoint()); //the actual vector defining the work we have to do
    q.applyOnTheRight(i - 1, i, rot);                  //this time, work on Q, not the
    r.applyOnTheLeft(i - 1, i, rot.adjoint());         //work on R
  }


  MatrixXd newq = q.block(1, 1, q.rows() - 1, q.cols() - 1);
  MatrixXd newr = r.block(1, 0, r.rows() - 1, r.cols());

  return newq * newr;
}

void HouseholderQR_delete_example()
{
  //Generate random test inputs
  MatrixXd lhs    = MatrixXd::Random(1001, 1000);
  MatrixXd rhs    = MatrixXd::Random(1001, 100);
  int rowToRemove = 65;

  HouseholderQR<MatrixXd> qr(lhs);


  //Time the fast version
  std::chrono::steady_clock::time_point fastTime = std::chrono::steady_clock::now();
  MatrixXd fastSolve                             = ReSolveWithRemovedRow(qr, rhs, rowToRemove);
  std::chrono::duration<double> fastDuration     = std::chrono::steady_clock::now() - fastTime;

  MatrixXd newLhs = RemoveARow(lhs, rowToRemove);
  MatrixXd newRhs = RemoveARow(rhs, rowToRemove);


  //Time with direct solving
  std::chrono::steady_clock::time_point slowTime = std::chrono::steady_clock::now();
  ColPivHouseholderQR<MatrixXd> fullLhsFactor(newLhs);
  MatrixXd slowSolve = fullLhsFactor.solve(newRhs);

  std::chrono::duration<double> slowDuration = std::chrono::steady_clock::now() - slowTime;

  cout << "Standard HouseholderQR: test matrix = " << lhs.rows() << " x " << lhs.cols() << endl;
  cout << "Complete refactoring time: " << std::chrono::duration_cast<std::chrono::microseconds>(slowDuration).count() /
  1.0e6 << " seconds." << endl;
  cout << "Fast update to QR to resolve: " <<
  std::chrono::duration_cast<std::chrono::microseconds>(fastDuration).count() / 1.0e6 << " seconds." << endl;

  cout << "Max coefficient error in solve, compared to resolve with ColPivQR: " <<
  (fastSolve - slowSolve).array().abs().matrix().maxCoeff() << endl;
  cout << "Max coefficient error reconstructing the LHS with updated QR " <<
  (newLhs - ApproximateQRWithRemovedRow(qr, rowToRemove)).array().abs().matrix().maxCoeff() << endl;
}

//Same as above, except use a ColPivHouseholderQR as the input. The only change in the
//algorithm to to apply the pivoting permutation at the end.
MatrixXd ReSolveWithRemovedRow(ColPivHouseholderQR<MatrixXd>& lhs_qr, MatrixXd const& rhs, int rowToRemove)
{
  MatrixXd r =  lhs_qr.matrixQR().template triangularView<Upper>();

  assert(r.rows() == rhs.rows());
  assert(rowToRemove < r.rows());
  assert(r.rows() > r.cols());
  assert(lhs_qr.rank() == r.cols()); //if not full rank, we might have to deal with zero pivots that become non-zero?

  MatrixXd newRhs = rhs;

  newRhs.applyOnTheLeft(lhs_qr.householderQ().transpose());

  //Get the row of Q we need by applying it to a unit vector
  VectorXd q_to_zero = VectorXd::Unit(rhs.rows(), rowToRemove);
  q_to_zero.applyOnTheLeft(lhs_qr.householderQ().transpose());

  JacobiRotation<double> rot;
  for (int i = q_to_zero.rows() - 1; i > 0; --i) {
    rot.makeGivens(q_to_zero(i - 1), q_to_zero(i));
    q_to_zero.applyOnTheLeft(i - 1, i, rot.adjoint()); //the actual vector defining the work we have to do
    newRhs.applyOnTheLeft(i - 1, i, rot.adjoint());    //rotate the rhs
    r.applyOnTheLeft(i - 1, i, rot.adjoint());         //work on R
  }


  r.block(1, 0, r.cols(),
          r.cols()).template triangularView<Upper>().solveInPlace(newRhs.block(1, 0, r.cols(), newRhs.cols()));


  MatrixXd dst(r.cols(), rhs.cols());

  //Now undo the column pivoting, offset by 1, since we're discarding the first row
  for (int i = 0; i < lhs_qr.cols(); ++i) {
    dst.row(lhs_qr.colsPermutation().indices().coeff(i)) = newRhs.row(i + 1);
  }

  return dst;
}

void ColPivHouseholderQR_delete_example()
{
  //Generate random test inputs

  MatrixXd lhs    = MatrixXd::Random(1001, 1000);
  MatrixXd rhs    = MatrixXd::Random(1001, 100);
  int rowToRemove = 65;

  ColPivHouseholderQR<MatrixXd> qr(lhs);


  //Time the fast version

  std::chrono::steady_clock::time_point fastTime = std::chrono::steady_clock::now();
  MatrixXd fastSolve                             = ReSolveWithRemovedRow(qr, rhs, rowToRemove);
  std::chrono::duration<double> fastDuration     = std::chrono::steady_clock::now() - fastTime;


  MatrixXd newLhs = RemoveARow(lhs, rowToRemove);
  MatrixXd newRhs = RemoveARow(rhs, rowToRemove);

  //Time with direct solving
  std::chrono::steady_clock::time_point slowTime = std::chrono::steady_clock::now();
  ColPivHouseholderQR<MatrixXd> fullLhsFactor(newLhs);
  std::chrono::duration<double> slowDuration = std::chrono::steady_clock::now() - slowTime;
  MatrixXd slowSolve                         = fullLhsFactor.solve(newRhs);

  cout << "ColPivHouseholderQR: test matrix = " << lhs.rows() << " x " << lhs.cols() << endl;
  cout << "Complete refactoring time: " << std::chrono::duration_cast<std::chrono::microseconds>(slowDuration).count() /
  1.0e6 << " seconds." << endl;
  cout << "Fast update to QR to resolve: " <<
  std::chrono::duration_cast<std::chrono::microseconds>(fastDuration).count() / 1.0e6 << " seconds." << endl;
  cout << "Max coefficient error in solve, compared to resolve with ColPivQR: " <<
  (fastSolve - slowSolve).array().abs().matrix().maxCoeff() << endl;
}

int main()
{
  HouseholderQR_delete_example();
  cout << endl;
  ColPivHouseholderQR_delete_example();
}