# Wrong generalized eigenvalues of large numeric matrices

I've been working on a toy model of a variational method of a Hamiltonian H with an overlap matrix S. The point is that when H and S are small enough, Mathematica gives the correct eigenvalues (≈ -1,-0.25,-0.111... hydrogen-like eigenvalues).

The problem is that when i increase the basis, the correct eigenvalues become zero. For example, if my basis is made of 34 states i have some correct eigenvalues of the exact problem and some others (positive ones) that are wrong, but are expected, however, if i increase the basis to 35 states a few of the correct eigenvalues become 0, and if the basis is, let's say 50 states, all the correct eigenvalues become 0.

I know i don't need such amount of states for this problem but the real problem i want to address is more complex than this one and i'd like to understand what is happening here.

Here i put the code i am using, S[i,j] and H[i,j] are the matrix elements of Sij[n] and Hij[n], where n is the size of the basis:

S[i_, j_] :=8 Sqrt[1.7^(3 i + 3 j - 6)]/(1.7^(i - 1) + 1.7^(j - 1))^3

H[i_, j_] := -8 Sqrt[1.7^(3 i + 3 j - 6)]/(1.7^(i - 1) + 1.7^(j -
1))^3*(0.01*(1.7^(i - 1) + 1.7^(j - 1)) - 0.01^2*1.7^(i + j - 2))

Sij[n_] := Table[S[i, j], {i, 1, n}, {j, 1, n}]

Hij[n_] := Table[H[i, j], {i, 1, n}, {j, 1, n}]

Eig[n_]:=Sort[Eigenvalues[{Hij[n],Sij[n]}]


I have tried not to use directly the Mathematica default Eigenvalues algorithm, for example using the inverted Sij[n] or the Cholesky decomposition, but i had similar results.

• Probably not an issue but: Do you know if the geometric multiplicity (of your zero eigenvalues)(geometric multiplicity= dimension of the space spanned by the eigenvectors) is equal to the algebraic multiplicity? If they're equal then increasing the precision is probably optimal. If they're unequal then some eigenvectors/adjoint eig vecs are orthogonal. The condition can be improved by lifting the problem into a higher dimensional space. Think of two orthogonal vectors in the z=0 plane, (0,1,0) and (1,0,0). Lift them up to say (0,1,1/2) (1,0,1/2). This can be done cleanly preserving spectrum.
– JEP
Jan 29 '15 at 1:21
• The algebraic multiplicity and the geometric multiplicity are equal. Thank you anyway. Jan 29 '15 at 9:52

This appears to be a precision issue. Increase the precision.

S[i_, j_] =
8 Sqrt[1.7^(3 i + 3 j - 6)]/
(1.7^(i - 1) + 1.7^(j - 1))^3 //
Rationalize // Simplify;

H[i_, j_] =
-8 Sqrt[1.7^(3 i + 3 j - 6)]/
(1.7^(i - 1) + 1.7^(j - 1))^3*
(0.01*(1.7^(i - 1) + 1.7^(j - 1)) -
0.01^2*1.7^(i + j - 2)) //
Rationalize // Simplify;

Sij[n_] := Array[S, {n, n}]

Hij[n_] := Array[H, {n, n}]

Eig[n_, workPrec_: 16] :=
Eigenvalues[N[{Hij[n], Sij[n]}, workPrec]] // Sort

Eig // NumberForm[#, 6] & Eig[50, 20] // NumberForm[#, 6] & Eig[50, 30] // NumberForm[#, 6] & 