Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I know that Eigenvalues is already quite well implemented in Mathematica, nor am I foolishly trying to improve on it. In order to improve my programming skills, I am trying to write Mathematica-style code to locate eigenvalues of a tridiagonal symmetric matrix using bisection. This is what I came up with.

tridiagbisec[diag_List, subdiag_List, tol_Real] := 
 Module[{diagonals = 
     Split[Transpose[{diag, Join[{0}, subdiag]}], #2[[2]] =!= 0 &]},
   Flatten[
    bisecnonzero[-myNorm[Sequence @@ Transpose[#]], 
       myNorm[Sequence @@ Transpose[#]], tol, #] & /@ diagonals]] /; 
  Length[diag] - Length[subdiag] == 1

tridiagbisec is the main function I will be calling. It builds an array (diagonals) from the arrays containing the matrix elements of both the leading diagonal and the subdiagonal, then Splits it wherever a zero is found so that each block is evaluated separately, using bisecnonzero.

bisecnonzero[\[Lambda]min_Real, \[Lambda]max_Real, tol_Real, diagonals_List] :=
  Module[{\[Lambda]med = (\[Lambda]min + \[Lambda]max)/2, 
   nmin = numeig[diagonals, \[Lambda]min], 
   nmax = numeig[diagonals, \[Lambda]max], nmed},
  nmed = numeig[diagonals, \[Lambda]med];
  {Which[(nmin > nmed) && (\[Lambda]med - \[Lambda]min > tol), 
    bisecnonzero[\[Lambda]min, \[Lambda]med, tol, diagonals],
    (nmin > nmed) && (\[Lambda]med - \[Lambda]min <= tol), 
    ConstantArray[(\[Lambda]med + \[Lambda]min)/2, nmin - nmed],
    True, {}],
   Which[(nmed > nmax) && (\[Lambda]max - \[Lambda]med > tol), 
    bisecnonzero[\[Lambda]med, \[Lambda]max, tol, diagonals],
    (nmed > nmax) && (\[Lambda]max - \[Lambda]med <= tol), 
    ConstantArray[(\[Lambda]max + \[Lambda]med)/2, nmed - nmax],
    True, {}]}
  ]

This, I believe, is where I may have built my program in a suboptimal way. Are two Which the right way to iterate bisection?

numeig[diagonals_, \[Lambda]_] := 
  numeig[diagonals, \[Lambda]] = 
   Unitize[#].UnitStep[#] &@
    Rest@FoldList[
      If[Not@PossibleZeroQ@#1, #2[[
          1]] - \[Lambda] - #2[[2]]^2/#1, +\[Infinity]] &, 1, diagonals];

numeig computes the number of eigenvalues greater than \[Lambda] (cfr. (5) in Barth, Martin and Wilkinson (1967).

myNorm = Max[Abs@#1 + Abs@#2 + RotateLeft@Abs@#2] &;

I know that I am also reimplementing Norm[#,Infinity] &, but for scholastic purposes I think it may be useful, since for a tridiagonal matrix it has a particularly simple form and this way I can avoid building a SparseArray structure at all.

Since I am trying to be performance-conscious, how could the whole algorithm be improved upon with reasonable effort? ...Besides using Eigenvalues, that is! :D

share|improve this question
    
re: putting together the diagonal and subdiagonal arrays, have you seen Riffle[]? I seem to recall writing a bisection routine in Mathematica a while back; let me see if I can find it... –  J. M. Jun 20 '12 at 8:04
    
Also, note that the EISPACK routine possesses a number of improvements from the original Handbook routine; you might want to consider emulating that instead of the Algol routine from the original article. –  J. M. Jun 20 '12 at 8:06
2  
You could use On["Packing"] to see if the alg. unpacks and try to eliminate those parts of the code. –  user21 Jun 20 '12 at 15:32
    
@ruebenko, I hadn't thought about that! Thanks! What's surprising, however, is that numeig unpacks! And it seems it's because FoldList unpacks! This does not look nice. –  Andrea Colonna Jun 20 '12 at 21:19
add comment

1 Answer

up vote 3 down vote accepted

As a starting point, here's a slight update on old code I wrote for implementing bisection:

n = 10;
d = Table[2 k - 1, {k, n}]; e = Table[k, {k, n - 1}]; (* Laguerre tridiagonal matrix *)

prec = 20;
(* emin and emax are bounds from Gerschgorin's theorem *)
emin = N[Min[Total /@ Partition[Riffle[d, -Abs[e]], 3, 2, {2, 1}, {}]], 3 prec/2];
emax = N[Max[Total /@ Partition[Riffle[d, Abs[e]], 3, 2, {2, 1}, {}]], 3 prec/2];

N[Table[
   a = emin; b = emax; h = Abs[b - a];
   While[h > 10^-prec,
          h /= 2; x = a + (b - a)/2; u = d[[1]] - x;
          k = Boole[Negative[u]];
          Do[
              If[u == 0, u = (e[[j]] + 10^-(2 prec)) 10^-(2 prec)];
              u = (d[[j + 1]] - Abs[e[[j]]]^2/u) - x;
              k += Boole[Negative[u]],
             {j, n - 1}];
          If[k < m, a = x, b = x];
         ];
   a + (b - a)/2,
  {m, n}], prec]
{0.1377934705404924308, 0.7294545495031704982, 1.808342901740316048,
 3.401433697854899515, 5.552496140063803633, 8.330152746764496700,
 11.843785837900065565, 16.27925783137810210, 21.996585811980761951,
 29.920697012273891560}

% == Sort[Eigenvalues[N[SparseArray[{Band[{1, 1}] -> d,
                                     Band[{2, 1}] -> e, Band[{1, 2}] -> e}], prec]]]
True

(Compare the eigenvalues with the output of \[FormalX] /. NSolve[LaguerreL[n, \[FormalX]], \[FormalX], prec] as well.)

The implementation is procedural, as you can see; it should be possible to make the code functional, but how to cleanly do so is escaping me at the moment. I'll edit this answer later if I think of something.

share|improve this answer
    
(P.S. I'm looking at RecurrenceTable[] now; maybe it can help a bit here.) –  J. M. Jun 20 '12 at 14:49
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.