Hello, I'm using Eigen for some Least Squares minimization. I can easily use the x = A.jacobiSVD(ComputeThinU | ComputeThinV).solve(b) form and it works perfectly. However, I'm getting a constant stream of new measurements, and since what I'm doing is an approximation, the old measurements become less valid as time goes on. I find it works well to use the last 60 or so, giving me a 60x6 sized A matrix.

My question is, could I save time by keeping the JacobiSVD around, and somehow replace the oldest measurement without recomputing the whole thing?

I am also looking at a problem that would be just adding new measurements, not replacing older ones.

Another idea I had was thinking doing a weighted least squares using the previous result with a weight of # of measurements and the new measurement with a weight of one. I don't know if that would work though.
IE: an identity for the first six rows of A, then the previous x as the first six rows of b.
But that wouldn't remove old measurements, I don't think. Unless I could put that with a weight of negative one. Just a thought.

I'm sure this has come up before, but my Google-fu is weak today.

Thanks

All right, I tried the Weighted Least Squares thing and it didn't work. That's always assuming I did it right, but I think so.

But! And here is the good news. I figured out a good set of search terms that found what I was looking for. It's called Recursive Least Squares. There are updating and downdating equations. In the interest of helping anyone who finds this in their google search, I'll put some of the useful references here.

http://home.utad.pt/~psal/Mestrado/ficheiros/lsq.pdf
https://www.cs.iastate.edu/~cs577/handouts/recursive-least-squares.pdf

I'll be reading these to figure out how to do it, but if someone already knows and has Eigen code that does it, I'd be happy to see it.

If your problem is well conditioned you can use LDLT::rankUpdate through the normal equation, something like (pseudo code):

Code: Select all

MatrixX6d A = ...;
VectorXd b = ...;
Matrix6d AtA = A^T * A;
Vector6d Atb =  A^T * b;

LDLT<Matrix6d> ldlt;
ldlt.compute(A);
Vector6d x = ldlt.solve(Atb);

while(...) {
    Vector6d a = "new measurement";
   double y = "new measurement";
    ldlt.rankUpdate(a);
   Atb += y*a;
    x = ldlt.solve(Atb);
   // you can also remove measurements:
   ldlt.rankUpdate(a,-1);
   Atb -= y*a;
}


However, this strategy is efficient only if you are adding a very few novel measurements, otherwise it is better to rebuild the normal equation from zero at once.

Actually, since your problem is very small (only 6 unknows), if you add many measurements at once, then better update the AtA matrix like the Atb vector and recompute the LDLT from the new AtA:

Code: Select all

AtA += v^T * v; // add
AtA -= w^T * w; // remove
ldlt.compute(AtA);


Here v and w can be matrices.