# IterativeSolvers and Positive Definite Matrices

**URL:** <https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292>\
**Category:** Numerics\
**Tags:** question\
**Created:** [November 17, 2020, 11:41am UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292 "2020-11-17T11:41:16Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jonas](https://avatars.discourse-cdn.com/v4/letter/j/aeb1de/32.png) [@Jonas](https://discourse.julialang.org/u/Jonas)\
**Post date:** [November 17, 2020, 11:41am UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/1 "2020-11-17T11:41:16Z")

</div>

Hello dear Julia Users,

In order to calculate eigenvalues and vectors I am currently using the package IterativeSolvers, but I often get the error that my input matrices are not positive definite. My code is like this:

```julia
result = lobpcg(K, G, false, nte1 ; maxiter = 200)

```

where K is the stiffness matrix of a mesh and G the mass matrix. So from a mathematical point of view the matrices have to be positive definite. Somehow the function lobpcg() disagrees and throws an error. Before IterativeSolvers, I used the function eigs() from the package Arpack and I never had any issue of that kind. My matrix implementations are correct, since I compared my results of eigs() to commercial solvers. Does anybody know, how I can solve this?

Best regards,  
Jonas

---

<div class="post-metadata">

**Author:** ![mfh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mfh/32/13787_2.png) [@mfh](https://discourse.julialang.org/u/mfh)\
**Post date:** [November 17, 2020, 6:54pm UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/2 "2020-11-17T18:54:04Z")

</div>

Hi Jonas

From what you describe it sounds a lot like this is related to some similar issues I had with IterativeSolver’s lobpcg some time ago in using it for our [DFTK](https://dftk.org) code. See [this issue](https://github.com/JuliaMath/IterativeSolvers.jl/issues/246) I opened up in the repo which describes the symptoms I observed.

Since in our applications we sometimes need rather low tolerances (1e-10 and below), we ended up coding up our own lobpcg with specifically numerical stability in mind. It’s unfortunately not yet properly packaged for use outside of DFTK at the moment, but the implementation in the file [src/eigen/lobpcg\_hyper\_impl.jl](https://github.com/JuliaMolSim/DFTK.jl/blob/master/src/eigen/lobpcg_hyper_impl.jl) of our repo should be fairly self-contained (modulo the `@timing` annotations which you can just remove). Perhaps it helps you to solve your problem. I should add, though, that we currently mainly use it for standard eigenproblems, so even though the code has worked for generalised eigenproblems at some point, that’s currently not tested in our CI, so that might well be failing.

Hope that helps  
Michael

---

<div class="post-metadata">

**Author:** ![Jonas](https://avatars.discourse-cdn.com/v4/letter/j/aeb1de/32.png) [@Jonas](https://discourse.julialang.org/u/Jonas)\
**Post date:** [November 18, 2020, 9:00am UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/3 "2020-11-18T09:00:29Z")

</div>

Hello Michael,

thank you very much for your answer, I will definitely look into your code 😃  
The strange thing is that even if I change nothing at all in the code (only the mesh might change slightly), the code sometimes is able to run without any problem and in the next run I get the error. This is the first time that annoyed, frequent running of the code does indeed benefit me haha.

---

<div class="post-metadata">

**Author:** ![mfh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mfh/32/13787_2.png) [@mfh](https://discourse.julialang.org/u/mfh)\
**Post date:** [November 18, 2020, 10:18am UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/4 "2020-11-18T10:18:55Z")

</div>

Yes I observed the same in the issue I reported. Sometimes you seem to hit magic numbers in the size of the matrix or something like that and it works and sometimes it does not. That’s the tricky part with numerical instabilities … if you have some randomness (e.g. in the guesses) sometimes you get lucky and sometimes not.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [November 18, 2020, 10:34am UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/5 "2020-11-18T10:34:02Z")

</div>

Try wrapping the matrix in `Symmetric`. A non-symmetric matrix in Julia is considered to be not positive definite. The matrix could be symmetric up to some small error which is what the `Symmetric` wrapper is for to tell Julia it is actually symmetric.

---

<div class="post-metadata">

**Author:** ![Jonas](https://avatars.discourse-cdn.com/v4/letter/j/aeb1de/32.png) [@Jonas](https://discourse.julialang.org/u/Jonas)\
**Post date:** [November 18, 2020, 11:42am UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/6 "2020-11-18T11:42:22Z")

</div>

Thank you for the idea, unfortunately this did not work either. I tried it with `Hermitian`, too, again with no positive result.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 18, 2020, 1:23pm UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/7 "2020-11-18T13:23:11Z")

</div>

> [@Jonas](#):
>
> So from a mathematical point of view the matrices have to be positive definite. Somehow the function lobpcg() disagrees and throws an error.

I haven’t looked at the `lobpcg()` code, but are you sure it’s complaining about your `K` and `G` matrices?

The LOBPCG algorithm builds up its own subspace basis B as it goes along, and my recollection is that if the implementor is not careful to occasionally re-orthogonalize this basis then it can become numerically ill-conditioned as the algorithm progresses. Maybe the `lobpcg()` code is missing this re-orthogonalization check and is running into a problem where B^T B (which it has to invert at some point) is nearly singular.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [November 18, 2020, 1:34pm UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/8 "2020-11-18T13:34:49Z")

</div>

That’s precisely the issue discussed in [https://github.com/JuliaMath/IterativeSolvers.jl/issues/246](https://github.com/JuliaMath/IterativeSolvers.jl/issues/246).

---

<div class="post-metadata">

**Author:** ![Jonas](https://avatars.discourse-cdn.com/v4/letter/j/aeb1de/32.png) [@Jonas](https://discourse.julialang.org/u/Jonas)\
**Post date:** [November 19, 2020, 10:00am UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/9 "2020-11-19T10:00:30Z")

</div>

Thank you for your effort.  
Best ragards.

---

<div class="post-metadata">

**Author:** ![Jonas](https://avatars.discourse-cdn.com/v4/letter/j/aeb1de/32.png) [@Jonas](https://discourse.julialang.org/u/Jonas)\
**Post date:** [November 19, 2020, 10:02am UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/10 "2020-11-19T10:02:00Z")

</div>

I think you are right. I will look into the discussed issues that mohamed82008 has mentioned.  
Best regards and thank you 😃

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [November 19, 2020, 10:49am UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/11 "2020-11-19T10:49:54Z")

</div>

You can also try this:  
[https://jutho.github.io/KrylovKit.jl/stable/man/eig/#Generalized-eigenvalue-problems-1](https://jutho.github.io/KrylovKit.jl/stable/man/eig/#Generalized-eigenvalue-problems-1)

Not the same as lobpcg but should be more general. Though it has been a while that I looked at this and it is not extensively tested.

---

<div class="post-metadata">

**Author:** ![Jonas](https://avatars.discourse-cdn.com/v4/letter/j/aeb1de/32.png) [@Jonas](https://discourse.julialang.org/u/Jonas)\
**Post date:** [November 20, 2020, 9:19am UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/12 "2020-11-20T09:19:46Z")

</div>

I will look into it, thanks 😃

---

<div class="post-metadata">

**Author:** ![LysandrosAn](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lysandrosan/32/12528_2.png) [@LysandrosAn](https://discourse.julialang.org/u/LysandrosAn)\
**Post date:** [November 22, 2020, 1:44pm UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/13 "2020-11-22T13:44:26Z")

</div>

Since it is a stiffness matrix K, the deformation vector u due to a force vector F should be:

u=inv(K)\*F.

Try to invert the stiffness matrix and multiply it with a force vector where only one element will be unit, the other being zero. The response vector u might give you a hint, on what part of the matrix needs to be corrected. It sounds as if the system is not constrained correctly and can freely move. Good luck

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 22, 2020, 3:54pm UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/14 "2020-11-22T15:54:43Z")

</div>

You worrie me when you say “it is not extensively tested”

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [November 22, 2020, 8:40pm UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/15 "2020-11-22T20:40:35Z")

</div>

By which I meant: there are tests to make sure that the correct result is obtained for simple problems, but I am not using the generalized eigenvalue problem often in my own applications. Furthermore, I am not specialized in numerical mathematics so maybe I am unaware of shortcomings or possible improvements.

---

<div class="post-metadata">

**Author:** ![Jonas](https://avatars.discourse-cdn.com/v4/letter/j/aeb1de/32.png) [@Jonas](https://discourse.julialang.org/u/Jonas)\
**Post date:** [November 24, 2020, 1:03pm UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/16 "2020-11-24T13:03:41Z")

</div>

Hello and thanks for your answer. The matrices in my problem seem to be correct. I am solving a generalized eigenvalue problem for determining the cutoff frequency of cylindrical waveguides. And in comparison with commerical tools, I get the same results. So as it was mentioned, the problem should lie in the used algorithm itself 😃

Best regards,  
Jonas

---

<div class="post-metadata">

**Author:** ![lobpcg](https://avatars.discourse-cdn.com/v4/letter/l/e19b73/32.png) [@lobpcg](https://discourse.julialang.org/u/lobpcg)\
**Post date:** [November 26, 2020, 10:37pm UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/17 "2020-11-26T22:37:24Z")

</div>

@Jonas I guess that you can try to use a Julia wrapper to [scipy.sparse.linalg.lobpcg — SciPy v1.9.1 Manual](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.lobpcg.html)  
which includes many extra tricks to improve stability compared to the current native Julia code.

---

<div class="post-metadata">

**Author:** ![Jonas](https://avatars.discourse-cdn.com/v4/letter/j/aeb1de/32.png) [@Jonas](https://discourse.julialang.org/u/Jonas)\
**Post date:** [November 27, 2020, 8:58am UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/18 "2020-11-27T08:58:44Z")

</div>

Thank you for that suggestion.

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [November 27, 2020, 9:36am UTC](https://discourse.julialang.org/t/iterativesolvers-and-positive-definite-matrices/50292/19 "2020-11-27T09:36:24Z")

</div>

Out of curiosity, have you tried KrylovKit.jl 's `geneigsolve` ?