lobpcg has spurious PosDefExceptions

[jpfairbanks]

I'm trying to use the lobpcg for the eigenvalues of adjacency matrices in LightGraphs.jl but I ran into some spurious PosDefExceptions in lobpcg.

Here is a MWE:

julia> A = randn(5,5)
       println("Eigenvalues: $(eigen(A'A).values))")
       println(cholesky(A'A))
       esol = lobpcg((A'A), true, 3)

Eigenvalues: [0.0902153, 1.42337, 5.49986, 13.4877, 17.4287])
Cholesky{Float64,Array{Float64,2}}([2.6168 -0.761333 -0.847569 -0.535911 1.36087; -1.99225 3.46792 0.3024 -0.222808 1.10018; -2.21792 1.69398 1.10087 -0.858487 -1.06607; -1.40237 -0.364672 -0.558235 0.691611 2.68809; 3.56111 2.77927 -1.99433 1.79989 1.86484], 'U', 0)
ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.

A'A is positive definite as show by the positive eigenvalues and successful call to cholesky(A'A)

[haampie]

Does this maybe have to do with the small size of the matrix @mohamed82008? Or some galerkin projection of the matrix that is nearly singular? Or missing a Symmetric / Hermitian wrapper?

[mohamed82008]

Very peculiarly, this only happens with matrices of size 5, not 6 not 4. Basic investigation shows that the residuals matrix is not full rank. This is probably a shortcoming of the LOBPCG algorithm with Cholesky based orthogonalization since there is no theoretical guarantee that all the basis will be independent. Ideally, I should try Cholesky, if it doesn't work then we should fallback on the Householder QR orthogonalization. I brought this up before in a Slack discussion with @haampie on inner and outer othogonalization of LOBPCG. This was supposed to go on a blog post for which a draft was started, but that never saw the light.

As a workaround you can do:

x0 = rand(5,1); esol = lobpcg((A'A), true, x0, 3);

This will decouple the block size, 1, from the number of vectors you want to find, so it will not be dealing with 3 residual vectors simultaneously. With a block size of 1, all the basis are guaranteed to be independent unless you are very very close to convergence, which you won't have to worry about if the tolerance is well over machine precision.

[mohamed82008]

Actually something else is also happening here. Even with using QR to orthogonalize if n < 3*bs, where bs is the block size, the Rayleigh-Ritz matrices will still be singular. It is not recommended to use this method with n < 3*bs.

To understand why this is let me say a thing about how LOBPCG works. LOBPCG tries to find the target eigenvectors in a reduced subspace of dimension bs in the first iteration, 2bs in the second and 3bs in the third iteration onwards. If the basis of this subspace are not independent then the algorithm is in trouble. This usually leads to failed inner or outer orthogonalization of the block basis. In some cases, I think it may also fail silently by not converging. So with bs = 3, if you are lucky enough and the algorithm converges in the first iteration, then you get your solution. If it goes into the second iteration, then there is no way the basis can be independent since these are 6 vectors of length 5. This will lead to a failure like the one you witnessed.

So when will the algorithm converge in the first iteration? If you are lucky enough, the initial random guess of 3 vectors will span the same subspace as the 3 eigenvectors you are looking for. If this happens, the algorithm converges in a single iteration. Of course the chances of this happening are pretty low. To "remedy" this, one thing to do is to never use a bs such that 3bs > n. LOBPCG was designed for huge matrices so this is a corner case that was not accounted for in the code. Using the case where n < bs also needs to be protected against, and if the random initial vectors are not independent. All these are corner cases that can lead to the failure of the algorithm. If bs << n, then these events are all very unlikely or not relevant.

[jpfairbanks]

Thanks for probing this for me. That makes sense, I don't need to solve 5x5 systems with LOBPCG I was just making an MWE.

Maybe checking for 3bs <= n and raising an exception telling the user to use a different solver would be a good way to go to prevent people from experiencing these random failures. The PosDefException definitely made me double check the definition of positive definite, which turned out to not be the problem at all.

[mohamed82008]

That sounds reasonable. Alternatively, we could resort to eigen for these small problems.

[lobpcg]

From https://www.mathworks.com/matlabcentral/fileexchange/48-lobpcg-m

% Checking if the problem is big enough for the algorithm, % i.e. n-sizeY > 5blockSize % Theoretically, the algorithm should be able to run if % n-sizeY > 3blockSize, % but the extreme cases might be unstable, so we use 5 instead of 3 here.

[mohamed82008]

@lobpcg I removed this check when porting the implementation because technically this branch would not be the "lobpcg" algorithm, maybe worth adding again. Thanks for the algorithm btw :)

[lobpcg]

Thank you for your kind words and for implementing lobpcg in Julia - the only one currently available, to my knowledge. Feel free to ask if you get any questions on lobpcg.