Implementing Lanczos algorithm in Mathematica

[fearsmathematic2015 (Feras Moh)]

I have to create a method or function based on Lanczo’s procedure. This method has to calculate eigenvalue and eigenvectors. The picture shows the algorithm that I am trying to use in Mathematica.
The idea of the algorithm can be explained as follows:
First an arbitrary array (q1) is selected and then a
new array (r0) calculated by this array. The (beta) is for eigenvectors and (k) is for iteration. The array (qk) ist calculated. With (Alpha) you calculate the eigenvalues. And with (Beta) the norm of the array (rk) is calculated so that (qk) becomes smaller and smaller and then (q[k+1]) is approximated to eigenvectors.

• I am not allowed to use the functions: “Eigenvalues, Eigensystem, Eigenvector”.

• I know my code has a lot of bugs in it, but maybe I can be helped.

• The problem is that I don’t get correct values of the eigenvectors and eigenvalues in
  comparison with the mathematicas above mentioned functions.

Also my code:

   A={{1,2,3},{2,4,5},{3,5,6}} 
   n=Length[A];
   q=RandomReal[{-1,1},n] ;
   r= q;
   BetaVe=1;
   While[BetaVe!=0.0000001,
   q= r/BetaVe ;
   AlphaVa=(q.A).q;
   r=((A-(AlphaVa *IdentityMatrix[n])).q)-(BetaVe*q);
   BetaVe = Norm[r,2];
   Print[{AlphaVa}]
   Print[q]
   ]

Update:

A={{1,2,3},{2,4,5},{3,5,6}} 
n=Length[A];
k=0.;
q[k_]:=q[0]=0.;
q[k+1]={0.55,0.63,0.5482700065};
r[k_]:=r[0]= q[k+1];
AlphaVa[k_] :=AlphaVa[0]=0.; 
BetaVe[k_]:=BetaVe[0]= 1.;
While[BetaVe[k]>0.000001,
q[k+1]= r[k]/BetaVe [k];
k=k+1;
AlphaVa[k]=(q[k].A).q[k];
r[k]=((A-(AlphaVa[k] *IdentityMatrix[n])).q[k])-(BetaVe[k-1]*q[k-1]);
BetaVe[k] =Norm[r[k]];

Print[{AlphaVa[k]}]
Print[q[k]]
];

The results of my Code are now as follows with Print[r[k]]

Mathematica’s function “Eigensystem”:

W={{1,2,3},{2,4,5},{3,5,6}} 
Eigensystem[W]

The results of Mathematica functions are as follows:

I think that the eigenvalues can be calculated or approximated better with the updated algorithm, but the eigenvectors cannot. could someone possibly tell me why it is?

Update:

A={{1,2,3},{2,4,5},{3,5,6}} 
n=Length[A];
k=0;
q=Normalize[{-1,1,1}];
qq={q};
alpha=0;
alphaList={alpha};
beta= 1;
While[beta>0.99,
alpha=q.(A.q);
q=Normalize[((A-(alpha *IdentityMatrix[n])).q)-(beta*qq[[-1]])];
r=((A-(alpha *IdentityMatrix[n])).q)-(beta*qq[[-1]]);
beta=Norm[r];
AppendTo[alphaList,alpha];
AppendTo[qq,q];
];
alphaList  // N
qq // N

The result :

Out[10]=
{0.,3.66667,5.40606}
Out[11]=
{{-0.57735,0.57735,0.57735},{0.905204,0.243709,0.348155},{-0.920022,0.250721,0.301161}}

But I have bad convergence again.

[trounev (Alexander Trounev)]

As it well known, Lanczos algorithm is unstable in general case. For stabilization we can use the Gram-Schmidt process to orthogonalize intermediate set of vectors. So we can use two step algorithm. First, we precompute with the Lancozs algorithm as follows

A = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}};
n = Length[A];
q0 = RandomInteger[{-1, 1}, n];
q[1] = q0/Norm[q0, 2]; r[0] = q[1]; b[0] = 1.; 
q[0] = ConstantArray[0, n]; k = 0;
Do[If[b[k] != 0, q[k + 1] = r[k]/b[k], Break[]];
 k = k + 1; a[k] = (q[k] . A) . q[k];
 r[k] = ((A - (a[k]*IdentityMatrix[n])) . q[k]) - (b[k - 1]*q[k - 1]);
 b[k] = Norm[r[k], 2];, {k, 0, 2 n}] 

Second, we use orthogonalization and put condition for break

Do[If[b[k] != 0., q[k + 1] = r[k]/b[k], Break[]];
 k = k + 1; a[k] = (q[k] . A) . q[k];
 r[k] = ((A - (a[k]*IdentityMatrix[n])) . q[k]) - (b[k - 1]*q[k - 1]);
  ort = Orthogonalize[Table[r[k - n + i], {i, n}], 
   Method -> "GramSchmidt"]; Do[r[k - n + i] == ort[[i]];, {i, n}];
 b[k] = Norm[r[k], 2]; id = k; 
 If[Sqrt[Sum[(a[k - 2 n + i] - a[k - n + i])^2, {i, n}]] < 10^-5, 
  Break[], Continue[]];, {k, n + 1, 100 n^2}]

In the next example we take (by random choice) q0={-1, 1, -1}. In this case number of iterations is id=23, and the result of the last 2n iterations is

Table[{a[k], q[k]}, {k, id - 2 n + 1, id}] // TableForm

We can compare result for eigenvalues with standard function Eigenvalues[] as follows

 {Sort[Table[a[k], {k, id - n + 1, id}]], 
 val = Eigenvalues[A // N] // Sort}

Out[]= {{-0.515729, 0.170912, 11.3448}, {-0.515729, 0.170915, 
  11.3448}}

Finally we plot result on every iteration to check convergence

Show[Plot[val, {x, 0, id}, Frame -> True, Axes -> False], 
 ListPlot[Table[a[k], {k, 1, id}], PlotRange -> All, 
  PlotStyle -> Black]]

[rohitnamjoshi]

Crossposted here.

[trounev (Alexander Trounev)]

For numerical result we need to convert matrix in numerical form, for example for A = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}} we use

 Eigensystem[A // N]

Out[]= {{11.3448, -0.515729, 
  0.170915}, {{-0.327985, -0.591009, -0.736976}, {-0.736976, \
-0.327985, 0.591009}, {0.591009, -0.736976, 0.327985}}}

[fearsmathematic2015 (Feras Moh)]

Thank you very much for your comment. I know the Lanczo’s Algorithmus is not stable. Therefore I have to discuss these procedures and compare them with Mathematica’s approach.