Added accelerated power method to compute eigenpair fast #533

jingwuOUO · 2020-09-20T15:31:12Z

No description provided.

david-m-rosen

Looks like you have all of the necessary pieces here, but this implementation looks to me like it's waaaaay more complicated than it needs to be.

I think it would be much better to refactor this into two simple standalone functions (i.e. not structs, which in my opinion are way overkill for such simple functions):

powerMethod(), which implements only the basic (optionally-accelerated) power method
computeMinimumEigenPair(), which implements Alg 6 in the paper (using powerMethod())

You should be able to implement a basic version of PowerMethod() in ~30 lines of code, and similarly for computeMinimumEigenPair()

Importantly: everything involving initialization, etc. should be done outside of the powerMethod() function -- this function should take as input the starting vector v0 (perhaps with an optional default random initialization of appropriate dimension, in case the user doesn't have, or doesn't want to supply, an explicit initial estimate).

gtsam/sfm/PowerMethod.h

gtsam/sfm/ShonanAveraging.cpp

dellaert · 2020-09-21T12:48:40Z

@david-m-rosen i will take over from @jingwuOUO for a bit but will take your comments into account.

dellaert · 2020-09-21T19:19:46Z

@jingwuOUO back to you. However, I'd like to speak to you first before you start editing, so let's schedule a meeting.

gtsam/sfm/PowerMethod.h

dellaert · 2020-09-21T19:51:17Z

Given our conversation, please:

remove the sorting
split up the PowerMethod into two classes (not structs):

a base class PowerMethod
a derived class AcceleratedPowerMethod

for PowerMethod, just have the one vector that is iterated
for AcceleratedPowerMethod, also add beta_ and the "previousVector_". If actually needed, that is: I think you can probably do without it.
replace S with an optional initial vector (for both classes)

dellaert

OK, This is starting to look good. Still, a large number of cleanup/documentation needed before this can go into GTSAM. I'd also like you to split the header into two files (and the tests) and move them to the linear sub-directory, as they are actually general and have nothing to do with sfm per se.

gtsam/sfm/tests/testPowerMethod.cpp

gtsam/sfm/PowerMethod.h

gtsam/sfm/ShonanAveraging.cpp

dellaert · 2020-10-13T18:18:09Z

Just a comment: please do not force-push a public branch, ever. It breaks the review process, and a number of other things. If this is about having the latest develop, justmerge develop into your branch.

dellaert

Still a number of things. Don't force-push

gtsam/linear/PowerMethod.h

gtsam/sfm/ShonanAveraging.cpp

david-m-rosen · 2020-10-14T02:34:09Z

@dellaert I'd also be happy to give this a close look whenever you have a version you're happy with, if you'd like a second set of eyes on it [I've not been tracking this closely for the past few weeks, as I got the sense that you and Jing wanted to take the lead on this for a bit :-) ].

dellaert · 2020-10-14T17:32:44Z

@dellaert I'd also be happy to give this a close look whenever you have a version you're happy with, if you'd like a second set of eyes on it [I've not been tracking this closely for the past few weeks, as I got the sense that you and Jing wanted to take the lead on this for a bit :-) ].

Yeah, my latest comments are just nits, so maybe you can take a look at the algorithmic aspects and comments? Would like to land this soon. Note, it's in linear subdirectory now and the intent is it is generic, not specific to Shonan. All Shonan-specific stuff should be in ShonanAveraging...

david-m-rosen

@dellaert @jingwuOUO OK I finished giving this another close look.

Overall the Shonan Averaging code looks much cleaner -- there are a few small changes to the minimum eigenvalue computation that should be done to speed up the method, but after having restructured the code, these amount to changing 1 or 2 lines. This part looks good to me.

The power method and accelerated power method code still needs a bit of work -- there's still weird stuff being done in a few places. The unit tests also need to be fixed -- right now the logic in these tests is actually completely broken (see comments).

gtsam/linear/PowerMethod.h

gtsam/linear/tests/testPowerMethod.cpp

gtsam/sfm/ShonanAveraging.cpp

jingwuOUO · 2020-10-22T19:49:13Z

@dellaert @jingwuOUO OK I finished giving this another close look.

Overall the Shonan Averaging code looks much cleaner -- there are a few small changes to the minimum eigenvalue computation that should be done to speed up the method, but after having restructured the code, these amount to changing 1 or 2 lines. This part looks good to me.

The power method and accelerated power method code still needs a bit of work -- there's still weird stuff being done in a few places. The unit tests also need to be fixed -- right now the logic in these tests is actually completely broken (see comments).

Yes. The test previous is not well-written. There are mistakes that I didn't fix. Now I have corrected the input with right value. For powermethod, the initial vector is a random one. And for accelerated, the initial is the row0 of aim matrix.

Also revised the code according to your comments.

david-m-rosen

OK, still a few things that need to be fixed, but I think this is getting pretty close now :-)

gtsam/linear/AcceleratedPowerMethod.h

david-m-rosen · 2020-10-23T21:19:23Z

gtsam/linear/AcceleratedPowerMethod.h

+
+    // initialize beta_
+    if (!initialBeta) {
+      estimateBeta();


I still don't agree with having estimateBeta() here -- I really don't think it's a good idea to have a heavy-duty operation like estimateBeta() be a part of the constructor.

Notice, for example, that it would be impossible to construct an instance of this class for later use without that code immediately running. In contrast, if you remove that one function, the constructor becomes basically a no-op.

I think it would be much better to have beta_ set directly using initialBeta, and leave estimateBeta() as a completely separate function that can be called by the user, if needed.

Update: I still think estimateBeta() should not be called here ...

@jingwuOUO is there a reason estimateBeta can't be moved out as a static function??

gtsam/linear/AcceleratedPowerMethod.h

david-m-rosen · 2020-10-23T23:14:11Z

gtsam/sfm/ShonanAveraging.cpp

+// (1) compute the maximum eigenpair (\lamda_dom, \vect{v}_dom) of A by power
+// method. if \lamda_dom is less than zero, then return the eigenpair. (2)
+// compute the maximum eigenpair (\theta, \vect{v}) of C = \lamda_dom * I - A by
+// accelerated power method. Then return (\lamda_dom - \theta, \vect{v}).
 static bool PowerMinimumEigenValue(
    const Sparse &A, const Matrix &S, double *minEigenValue,


OK this looks pretty good now, but I think it might be a bit cleaner to use double & instead of double* for minEigenValue and size_t& instead of size_t* for numIterations

(I like to avoid pointers whenever I can -- lots of creative ways to cause horrible segfaults with these :-P)

I changed this part. I'm not 100% sure if there is further consideration to use double* than double & and size_t because this form follows previous code by @dellaert . Could you take a look?

gtsam/linear/PowerMethod.h

david-m-rosen · 2020-10-23T23:22:57Z

gtsam/linear/tests/testAcceleratedPowerMethod.cpp

+  A.coeffRef(0, 0) = 2;
+  A.coeffRef(0, 0) = 1;
+  Vector initial = Vector6::Zero();
+  const Vector6 x1 = (Vector(6) << 1.0, 0.0, 0.0, 0.0, 0.0, 0.0).finished();


It might be better here to make initial a hard-coded vector instead of a random sample -- this would make it easier to reproduce failing tests at different times, different machines, etc.

gtsam/linear/tests/testPowerMethod.cpp

gtsam/linear/PowerMethod.h

david-m-rosen

A few final edits (all minor at this point), and I think this will be ready to go.

david-m-rosen · 2020-10-26T20:04:25Z

gtsam/linear/AcceleratedPowerMethod.h

-  // Update the ritzVector with beta and previous two ritzVector, and return
-  // x_{k+1} = A * x_k - \beta * x_{k-1} / || A * x_k - \beta * x_{k-1} ||
-  Vector powerIteration(const Vector &x1, const Vector &x0,
+  // Run accelerated power iteration on the ritzVector with beta and previous


You should change this description so that the docstring matches the function arguments -- use x_1, x_0, etc.

Also added parenthesis in numerator.

I still see a mismatch between doc and argument names. Also, should be a doxygen comment, /** */

gtsam/linear/AcceleratedPowerMethod.h

david-m-rosen · 2020-10-26T20:09:55Z

gtsam/linear/AcceleratedPowerMethod.h

-      Vector x0 = Vector::Random(this->dim_);
-      x0.normalize();
-      Vector x00 = Vector::Zero(this->dim_);
+    size_t T = 10;


You might consider making the parameter T an argument of this function (maybe with the default value 10)

gtsam/linear/AcceleratedPowerMethod.h

david-m-rosen · 2020-10-26T20:44:36Z

gtsam/linear/tests/testAcceleratedPowerMethod.cpp

+}
+
+/* ************************************************************************* */
+TEST(AcceleratedPowerMethod, useFactorGraph) {


This test still looks excessively complicated to me to be a good unit test ...

david-m-rosen · 2020-10-26T20:47:42Z

gtsam/linear/tests/testAcceleratedPowerMethod.cpp

+  // Check that we get zero eigenvalue and "constant" eigenvector
+  EXPECT_DOUBLES_EQUAL(0.0, solver.eigenvalues()[0].real(), 1e-9);
+  auto v0 = solver.eigenvectors().col(0);
+  for (size_t j = 0; j < 3; j++) EXPECT_DOUBLES_EQUAL(-0.5, v0[j].real(), 1e-9);


Again, you should never enforce getting a specific eigenvector

david-m-rosen · 2020-10-26T20:49:03Z

gtsam/linear/tests/testPowerMethod.cpp

+  // Check that we get zero eigenvalue and "constant" eigenvector
+  EXPECT_DOUBLES_EQUAL(0.0, solver.eigenvalues()[0].real(), 1e-9);
+  auto v0 = solver.eigenvectors().col(0);
+  for (size_t j = 0; j < 3; j++) EXPECT_DOUBLES_EQUAL(-0.5, v0[j].real(), 1e-9);


Same as before: you should never enforce getting a specific eigenvector

david-m-rosen · 2020-10-26T20:50:17Z

gtsam/linear/tests/testPowerMethod.cpp

+  pf.compute(50, 1e-4);
+  EXPECT_DOUBLES_EQUAL(ev1, pf.eigenvalue(), 1e-8);
+  auto actual2 = pf.eigenvector();
+  EXPECT(assert_equal(ev2, actual2, 3e-5));


Again, you should never enforce getting a specific eigenvector

About why using this unittest, the idea comes from that we want to use the graph.hessian() later to extract the A matrix. There are some tests using the EigenSolver from Eigen to see whether the eigenpair result is always "constant" and also to compare the results of power/accelerated method at initial (but now more clear about the correctness of power method).
But I'm still a little confused here, for the specific matrix A in this unittest, there is no two eigenvectors for one eigenvalue, and the eigenvector is normalized, the only uncertainty here is the sign of the eigenvector, whether +/-. So I guess it's fine to keep this. I'm not sure. I'm ok with deleting the checking specific eigenvector part.

The problem with:

EXPECT(assert_equal(ev2, actual2, 3e-5));

is that you are trying to enforce the condition that the (unit) eigenvector actual2 returned by the power method is exactly equal to the eigenvector ev2 returned by Eigen. But there are always at least 2 unit eigenvectors associated with any given eigenvalue, and in general there will be infinitely many (if that eigenvalue happens to have multiplicity > 1).

So the question here is: why would we ever expect the power method to return that specific eigenvector?

The answer is: we wouldn't, because eigenvectors are never unique. And if they're not unique, we shouldn't be trying to enforce uniqueness in a unit test. This doesn't make mathematical sense.

The correct test of "eigenvector-ness" is the eigenvector equation A*x = lambda*x. Specifically, we should be checking something like a normalized Ritz residual:

eps := ||A*x - lambda*x|| / ||x||

Note that this error measure does not depend upon which specific eigenvector we get, since any eigenvector x for eigenvalue lambda would have eps = 0.

david-m-rosen · 2020-10-26T20:57:42Z

gtsam/sfm/ShonanAveraging.cpp

+// ref : Part III. C, Rosen, D. and Carlone, L., 2017, September. Computational
+// enhancements for certifiably correct SLAM. In Proceedings of the
+// International Conference on Intelligent Robots and Systems.
+Vector perturb(const Vector &initialVector) {


Ah ok got it -- I missed this earlier.

But this construction looks a bit weird to me: you have perturb() as a function in ShonanAveraging.cpp, but not in ShonanAveraging.h -- this means that this function is only "visible" to functions in ShonanAveraging.cpp that occur lower down in this file.

I could easily see a case in which you might want to use this functionality outside of ShonanAveraging.cpp -- for example, it looks like you do this kind of perturbation in some of the unit tests?

I would suggest exposing this function externally by including its declaration in ShonanAveraging.h

jingwuOUO · 2020-11-02T01:50:37Z

Updated according to last review by David.

Something I'm not sure about. I see previous SparseMinimumEigenValue is not included in Shonan, so I did the same with PowerMinimumEigenValue, not sure why. Do we need to add PowerMinimumEigenValue func into Shonan so that I can also add purturb func in Shonan. @dellaert

dellaert

Still some comments by David that have not been addressed. Also, please read for 5 minutes about doxygen and then please apply that to entire PR.

gtsam/linear/AcceleratedPowerMethod.h

dellaert · 2020-11-09T13:39:12Z

gtsam/linear/AcceleratedPowerMethod.h

+
+    // initialize beta_
+    if (!initialBeta) {
+      estimateBeta();


@jingwuOUO is there a reason estimateBeta can't be moved out as a static function??

dellaert · 2020-11-09T13:40:16Z

gtsam/linear/AcceleratedPowerMethod.h

-  // Update the ritzVector with beta and previous two ritzVector, and return
-  // x_{k+1} = A * x_k - \beta * x_{k-1} / || A * x_k - \beta * x_{k-1} ||
-  Vector powerIteration(const Vector &x1, const Vector &x0,
+  // Run accelerated power iteration on the ritzVector with beta and previous


I still see a mismatch between doc and argument names. Also, should be a doxygen comment, /** */

gtsam/linear/AcceleratedPowerMethod.h

gtsam/linear/PowerMethod.h

dellaert · 2020-11-09T13:56:56Z

gtsam/linear/tests/testAcceleratedPowerMethod.cpp

+}
+
+/* ************************************************************************* */
+TEST(AcceleratedPowerMethod, useFactorGraph) {


It's a test for factor graph operator, at least that part should be tested somewhere. We should have tests for a sparse and dense matrix, as well. Simplicity would be good.

gtsam/linear/tests/testPowerMethod.cpp

jingwuOUO · 2020-11-16T04:00:56Z

So we need to add unittest for sparse and dense matrix, can I create some sparse and dense matrix myself? And then add it in these two test file.

dellaert · 2020-11-30T15:51:45Z

So we need to add unittest for sparse and dense matrix, can I create some sparse and dense matrix myself? And then add it in these two test file.

Addressed in meeting. Please re-request review.

dellaert

I will approve but trust you to address the final few comments I have.

gtsam/linear/tests/testPowerMethod.cpp

dellaert · 2020-12-01T04:14:52Z

gtsam/linear/PowerMethod.h

+ * \brief Compute maximum Eigenpair with power method
+ *
+ * References :
+ * 1) Rosen, D. and Carlone, L., 2017, September. Computational


Note sure this reference is appropriate for power method. More for accelerated. Power method itself is as old as the street.

Added Matrix Computation as reference and noted pages, which is also cited by Yunlun_arxiv paper

dellaert · 2020-12-01T04:16:02Z

gtsam/linear/AcceleratedPowerMethod.h

+  }
+
+  /// Tuning the momentum beta using the Best Heavy Ball algorithm in Ref(3)
+  double estimateBeta() const {


Explain again why this is still a method? I thought David had asked to make this a free function or at least a static method.

@dellaert @jingwuOUO I think it's fine to have this as a method (in fact, probably necessary, since we need access to the matrix-vector multiplication operator that's passed in during construction). What I don't want is this function being called inside the constructor (as currently happens in line 71), since this is a heavy-duty operation that is not actually required in order to construct an instance of this object -- permitting a default value of beta = 0.0 is perfectly fine, and I think arguably the "right" default behavior, since that will exactly correspond to the usual (simple) power iteration

dellaert · 2020-12-01T04:16:47Z

gtsam/sfm/ShonanAveraging.cpp

+// ref : Part III. C, Rosen, D. and Carlone, L., 2017, September. Computational
+// enhancements for certifiably correct SLAM. In Proceedings of the
+// International Conference on Intelligent Robots and Systems.
+Vector perturb(const Vector &initialVector) {


Actually, if it's a local function this can just be made static in this file.

david-m-rosen · 2020-12-13T22:33:12Z

Hi @jingwuOUO, the only changes I see that remain to be done here are:

Adding the number of iterations T used in estimateBeta() as an (optional) argument (with default value e.g. 10)
Not calling estimateBeta() from within the constructor for the accelerated power method, since this is a heavy-duty operation. It would be better to allow the user to specify beta directly, as an optional argument with a default value of 0 (corresponding to a standard power iteration). The user can always call estimateBeta() themselves if they're not sure what a good value of beta is.

dellaert · 2021-01-17T22:55:40Z

@jingwuOUO pls fix last two issues so we can finally merge?

david-m-rosen

Alright, looks good to me 😎

jingwuOUO added the feature New proposed feature label Sep 20, 2020

david-m-rosen requested changes Sep 20, 2020

View reviewed changes

dellaert reviewed Sep 21, 2020

View reviewed changes

gtsam/sfm/PowerMethod.h Outdated Show resolved Hide resolved

gtsam/sfm/PowerMethod.h Outdated Show resolved Hide resolved

gtsam/sfm/PowerMethod.h Outdated Show resolved Hide resolved

varunagrawal assigned jingwuOUO Sep 23, 2020

jingwuOUO force-pushed the jing/powermethod branch from 8f3fbd2 to 86255f6 Compare September 30, 2020 15:28

jingwuOUO requested a review from dellaert September 30, 2020 15:29

dellaert requested changes Oct 9, 2020

View reviewed changes

jingwuOUO added 6 commits October 11, 2020 02:15

Add acc power method in Alg6

2afbaaa

Completed and tested power method with unittest

e02633c

Added PowerMinimumEigenValue method and unittest.

fc112a4

feat: Add Best Heavy Ball alg to set beta

d6e2546

Refactor power and accelerated method

758c4b0

Reformat code.

b7f29a0

jingwuOUO force-pushed the jing/powermethod branch from 415bdef to b7f29a0 Compare October 13, 2020 06:34

Small fix

fbebd3e

dellaert requested changes Oct 13, 2020

View reviewed changes

david-m-rosen requested changes Oct 16, 2020

View reviewed changes

jingwuOUO added 2 commits October 22, 2020 14:19

Revise according to Frank and David's comments

5f50e74

Correct unittest input error

f86ad95

jingwuOUO added 2 commits October 22, 2020 15:51

Added more detailed documentation

32bf81e

Fixed not solve errors and detailed documentation

56300ca

jingwuOUO requested review from dellaert and david-m-rosen October 22, 2020 22:06

david-m-rosen requested changes Oct 23, 2020

View reviewed changes

jingwuOUO added 2 commits October 26, 2020 15:40

Revised as David's second review

cf813c5

Refined unittest

4ee4014

jingwuOUO requested a review from david-m-rosen October 26, 2020 19:41

david-m-rosen requested changes Oct 26, 2020

View reviewed changes

jingwuOUO added 2 commits November 1, 2020 20:45

Delete forcing compare eigenvector in unittest

f604a97

Revised documentation

70cecb3

Fixed forcing comparing eigenvector.

abfc98e

dellaert requested changes Nov 9, 2020

View reviewed changes

Fixed doxygen

6d9f95d

jingwuOUO added 2 commits November 30, 2020 18:52

Added dense matrix test case in power/acc

844cbea

Merge branch 'develop' into jing/powermethod

420272e

jingwuOUO requested a review from dellaert December 1, 2020 03:36

dellaert approved these changes Dec 1, 2020

View reviewed changes

jingwuOUO added 3 commits December 11, 2020 01:01

Added the example graph in powerMethodExample.h

afb6ebb

Refined reference documentation

55ad1f1

Make purturb static

525dbb6

Set estimateBeta() as optional

9c781b6

jingwuOUO requested a review from david-m-rosen January 19, 2021 15:08

david-m-rosen approved these changes Jan 21, 2021

View reviewed changes

jingwuOUO merged commit d267368 into develop Jan 21, 2021

dellaert deleted the jing/powermethod branch August 22, 2021 15:18

Added accelerated power method to compute eigenpair fast #533

Added accelerated power method to compute eigenpair fast #533

Conversation

jingwuOUO commented Sep 20, 2020

david-m-rosen left a comment

Choose a reason for hiding this comment

dellaert commented Sep 21, 2020

dellaert commented Sep 21, 2020

dellaert commented Sep 21, 2020

dellaert left a comment

Choose a reason for hiding this comment

dellaert commented Oct 13, 2020

dellaert left a comment

Choose a reason for hiding this comment

david-m-rosen commented Oct 14, 2020

dellaert commented Oct 14, 2020

david-m-rosen left a comment

Choose a reason for hiding this comment

jingwuOUO commented Oct 22, 2020

david-m-rosen left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

david-m-rosen left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jingwuOUO commented Nov 2, 2020

dellaert left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jingwuOUO commented Nov 16, 2020

dellaert commented Nov 30, 2020

dellaert left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

david-m-rosen commented Dec 13, 2020 • edited Loading

dellaert commented Jan 17, 2021

david-m-rosen left a comment

Choose a reason for hiding this comment

david-m-rosen commented Dec 13, 2020 •

edited

Loading