LOBPCG Algorithm 1 from Duersch et al, 2018 appears to simply be the original Algorithm 5.1 from A. V. Knyazev, Toward the Optimal Preconditioned Eigensolver: Locally Optimal Block Preconditioned Conjugate Gradient Method. SIAM Journal on Scientific Computing 23 (2001), no. 2, pp. 517-541. http://dx.doi.org/10.1137/S1064827500366124 It is probably appropriate the cite the original paper...

Also, this implementation may become unstable for various reasons.
More stable, but still simple, versions are described at
I. Lashuk, M. E. Argentati, E. Ovchinnikov and A. V. Knyazev, Preconditioned Eigensolver LOBPCG in hypre and PETSc. In Lecture Notes in Computational Science and Engineering, v. 55 (Proceedings of the 16th International Conference on Domain Decomposition Methods. 2005). Editors: O.~B. Widlund and D.~E. Keyes, pages 635-642. Springer, Berlin, 2007. ISBN: 3-540-34468-3. (http://math.ucdenver.edu/~aknyazev/research/papers/DD16/knyazev_contrib_revised/dd16.pdf)
or simply follow the latest code
https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py

@lobpcg Thanks for your comment!
Yes, I am aware of the possible instability issues with the basic version of LOBPCG algorithm.
I implemented it for reference and plan to implement also the more robust versions of the LOBPCG.
IIUC, the stability of lobpcg in SciPy is achieved by "soft-locking" the eigenpairs that have converged.
Using soft-locking for batches of eigenproblems is problematic as different batches may have different convergences scenarios and so the soft-locking feature might be made available only for dense and sparse matrix cases.

With the current sparse test matrices, the basic lobpcg works ok. While there are several (large) sparse matrices available for which the basic lobpcg might have problems, using these in pytorch test-site might be unfeasible due to their size. So, I wonder if you could provide any hints for how to generate relatively small sparse matrices where the basic lobpcg is guaranteed to fail?

It is not a matter of being sparse or dense, or even large/small. The instability mostly appears when the vectors in the RR become linearly dependent for any reason, unless one runs expensive explicit full orthogonalization prior to RR.

Easy ways to break RR:

1. run long enough with unreasonably small tolerance, especially in single precision
2. choose the block size large enough, e.g., n/6, where n is the matrix size
3. run in a loop with random initial guesses - some will fail
4. run with random initial guesses, but using a poor random number generator with a short cycle, relative to n, which thus creates a special block structure in the initial guess

And there are few trickier ones, with special choices of the initial guesses.

It is difficult, but still possible to break the latest code
https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py
so you can use it as a base for comparison when your code fails.

@pytorchbot rebase this please

💊 CircleCI build failures summary and remediations

As of commit 76c3cfb (more details on the Dr. CI page):

Commit 76c3cfb was recently pushed. Waiting for builds...

The PR is ready for review.

This PR implementations of LOBPCG do not follow
https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py
which is stable, while still avoiding any expensive QR-based orthogonalization.

Not a blocker, but may want to keep this in mind if there would be a need in the future to improve performance for large block sizes.

Thanks for your implementations!

I'm working on finding some good reviewers for this PR.

@vincentqb has graciously accepted to review this PR. @nikitaved, if you also want to bill Quansight time to help review this PR, that would also be fine by me.

@lobpcg How does this diff look to you? Happy to assign you as a reviewer too.

@ezyang I support this PR, but unfamiliar with pytorch. And this PR code is rather long, so I am unable to formally review, but thanks for asking.

I've looked briefly at a couple of core functions, and not noticed anything troubling, just had a comment #29488 (comment)

@pearu it may also help reviewers quite a bit if we split this PR into one per piece of functionality, if this is easily possible (I haven't read it yet so I don't know if there's some dependence)

Indeed, this PR implements two new features, low-rank SVD and LOBPCG, which are independent methods. However, the expected expertise of reviewers is the same and both features depend on some new utility methods needed for testing (see tests/common_utils.py), and some minor fixes in torch/_six.py and torch/serialization.py.

That's cool. In that case the utility of splitting would mostly be being able to land one part earlier than the rest, but if we're not in a rush let's just keep it together.

Thanks @vincentqb for the review of the SVD/PCA part. I have addressed all your questions and concerns. There is an open question of whether we should keep the complex input support or not because the pytorch complex support is incomplete currently.

@vincentqb Update from a JIT developer: apparently we're NOT supposed to put torch.jit.script annotations in core library #33418 (comment)

So if deleting that annotation suffices, that looks good too.

Let's remove the annotation for now. We can open a separate pull request to make this jitable, and discuss why decorating doesn't work.

@vincentqb I believe my last commit fixes the issue: the decorator is removed and lobpcg is torchscript jit'able.

Phabricator is refusing to update the original internal diff. We'll need to open a new PR to land the change. I'll close this one in the meanwhile.

@pearu -- can you push the change to a new PR, and link to this one? I'll approve it and merge it from there.

@vincentqb done creating a new PR.
Thanks!

@pearu - Amazing work, thanks a lot!
Used it for LBO decompositions, and compared it to Reuter's 2009.
Thank you!

@pearu, @ezyang, @vincentqb What do you think about some advertisement for this method? I could write, for example, a tutorial on Spectral Clustering with different Graph-Laplacians for MNIST/Fashion-MNIST and use LOBPCG/SVD/whatever for it. Could be useful for benchamarking as well...

@nikitaved @lobpcg -- Are you interested in writing about this? :)

@vincentqb , absolutely!

Great! I'm adding @wookim3 to the loop.

this is really cool, thanks for working on this!!

I do not think I can write it, never running pytorch myself, but I would be glad to review/edit. Please feel free to ping me.

@nikitaved - Thank you! To get started, can you share your email address? We can create a Google doc and collaborate together on there.

cc @vincentqb

@wookim3, sure, [email protected] it is. Thank you! But just so it is clear, it is a tutorial we are talking about, right?