# 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. 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 incomparison 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.
4 Replies
Sort By:
Posted 1 month ago
 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]]
Posted 1 month ago
 Crossposted here.