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 Split
s it wherever a zero is found so that each block is evaluated separately, using bisecnonzero
.
bisecnonzero[λmin_Real, λmax_Real, tol_Real, diagonals_List] :=
Module[{λmed = (λmin + λmax)/2,
nmin = numeig[diagonals, λmin],
nmax = numeig[diagonals, λmax], nmed},
nmed = numeig[diagonals, λmed];
{Which[(nmin > nmed) && (λmed - λmin > tol),
bisecnonzero[λmin, λmed, tol, diagonals],
(nmin > nmed) && (λmed - λmin <= tol),
ConstantArray[(λmed + λmin)/2, nmin - nmed],
True, {}],
Which[(nmed > nmax) && (λmax - λmed > tol),
bisecnonzero[λmed, λmax, tol, diagonals],
(nmed > nmax) && (λmax - λmed <= tol),
ConstantArray[(λmax + λ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_, λ_] :=
numeig[diagonals, λ] =
Unitize[#].UnitStep[#] & @
Rest @ FoldList[
If[Not @ PossibleZeroQ @ #1, #2[[
1]] - λ - #2[[2]]^2/#1, +∞] &, 1, diagonals];
numeig
computes the number of eigenvalues greater than λ
(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[#, ∞] &
, 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
Riffle[]
? I seem to recall writing a bisection routine in Mathematica a while back; let me see if I can find it... $\endgroup$ – J. M.'s ennui♦ Jun 20 '12 at 8:04On["Packing"]
to see if the algorithm unpacks and try to eliminate those parts of the code. $\endgroup$ – user21 Jun 20 '12 at 15:32numeig
unpacks! And it seems it's becauseFoldList
unpacks! This does not look nice. $\endgroup$ – Editortoise-Composerpent Jun 20 '12 at 21:19