# Tridiagonal symmetric matrix eigenvalue using bisection

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

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

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.

-
(P.S. I'm looking at RecurrenceTable[] now; maybe it can help a bit here.) – J. M. Jun 20 '12 at 14:49