Tridiagonal symmetric matrix eigenvalue using bisection

Question

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

Answer 1

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.