ArnoldiOp.h, negative vnorm

[hverhelst]

Hi,

I have some problems with using Spectra on the following matrices:

A = 54.821043,0.000000,-20.557891,10.278946
0.000000,109.642086,20.557891,30.836837
-20.557891,20.557891,356.336781,-171.315760
10.278946,30.836837,-171.315760,359.077833

B = -12.041839,-5.880219,18.558189,-23.373410
-5.880499,-31.779469,-18.670997,-29.203559
17.349251,-18.077255,-342.592370,170.787904
-21.685584,-27.818669,170.786817,-343.479375

I am using the ShiftInvert method for my computations and something seems to go wrong in the Arnoldi iterations. In particular, the TridiagEigen fails because of failing Arnoldi iterations.

It goes wrong in line 149 of Lanczos.h among other places, caused by line 83 of ArnoldiOp.h where the sqrt gets a negative value.

I looked a bit in the theory and from Algorithm 2 on page 51 of the ARPACK User's Guide and from page 44 of the book, it seems that the norm ||v||=sqrt(v.transpose()*v) is 'just' defined as the Euclidean 2 norm. However, in many places in Arnoldi.h, the ArnoldiOp::norm() is used, which is the norm sqrt(v.transpose()*A*v) . I am wondering what the reason is to use this norm, because it seems to cause trouble if v.transpose()*A*v<0. I tried to change it locally in some places, but my iterations now do not seem to converge, so maybe I am wrong with the changes.

Thanks for your help!

[yixuan]

Hi @hverhelst, page 44 and 51 (assuming they are pages of the book, not PDF) of ARPACK User's Guide are talking about regular eigenvalue problems, and Section 4.5 in page 59 is on the generalized eigenvalue problem.

And in page 60, there is one paragraph:

If $A$ is symmetric then one can maintain symmetry in the Arnoldi/Lanczos process by taking the inner product to be $\langle x,y \rangle=x^H My$.

This is essentially the reason why the norm $\Vert x \Vert_B=\sqrt{x'Bx}$ is used.

If you use SymGEigsSolver, then B must be positive definite. If A is positive definite but B is not, then you need to use SymGEigsShiftSolver.

In fact, at least one of A and B need to be positive definite for current solvers.

[yixuan]

After taking a look at your matrices, it turns out that A is positive definite and B is negative definite. You can simply work on -B, which is positive definite.