Message Boards Message Boards

GROUPS:

Implementing Lanczos algorithm in Mathematica

Posted 1 month ago
458 Views
|
4 Replies
|
8 Total Likes
|

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. Lanczos Algorithmus

  • 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]]

Result my Code

Mathematica's function "Eigensystem":

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

The results of Mathematica functions are as follows:

Resualt of function

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.

4 Replies

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]]

Figure 1

Posted 1 month ago

Crossposted here.

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}}}
Posted 26 days ago

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.

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract