Implements Alan's fast histogrammer for SymTridiagonal matrices

Author: jiahao

src/FastHistogram.jl

@@ -0,0 +1,57 @@
+#Uses the method of Sturm sequences to compute the histogram of eigenvalues of a matrix
+#
+# Reference
+#   Cy P. Chan, "Sturm Sequences and the Eigenvalue Distribution of
+#   the Beta-Hermite Random Matrix Ensemble", M.Sc. thesis (MIT), 2007
+#
+# For better numerical stability, this algorithm is modified so that only the sign
+# of the determinant is stored.
+#
+#For general matrices, this is slower than hist(eigvals(M), bins) and is NOT recommended
+function hist_eig{GridPoint <: Number}(M::AbstractMatrix, bins::Vector{GridPoint})
+    n = size(M)[1]
+    NumBins = length(bins)
+    histogram = zeros(NumBins)
+    SturmSequence = zeros(n)
+    for BinId=1:NumBins
+        K = M - bins[BinId]*eye(n)
+        #Progression of determinants of lower-right submatrices
+        #SturmSequence = [det(K[end-j+1:end,end-j+1:end]) for j=1:n]
+        #SturmRatioSequence = [SturmSequence[j+1]/SturmSequence[j] for j=1:n]
+
+        SturmSequenceSigns = [int(sign(det(K[end-j+1:end,end-j+1:end]))) for j=1:n]
+        SturmSequenceSigns = [1; SturmSequenceSigns]
+        SturmRatioSequence = [SturmSequenceSigns[j+1]/SturmSequenceSigns[j] for j=1:n]
+        histogram[BinId] = sum([r < 0 for r in SturmRatioSequence])
+    end
+    histogram = int(diff([histogram; n]))
+    return histogram
+end
+
+
+
+#Uses the method of Sturm sequences to compute the histogram of eigenvalues of
+#a symmetric tridiagonal mmatrix
+#
+# Reference
+#   Cy P. Chan, "Sturm Sequences and the Eigenvalue Distribution of
+#   the Beta-Hermite Random Matrix Ensemble", M.Sc. thesis (MIT), 2007
+#
+function hist_eig{GridPoint <: Number}(M::SymTridiagonal, bins::Vector{GridPoint})
+    n = size(M)[1]
+    NumBins = length(bins)
+    histogram = zeros(NumBins)
+    r = zeros(n)
+    for BinId=1:NumBins
+        a, b = M.dv - bins[BinId], M.ev
+        #Formula (2.3)
+        r[1] = a[1]
+        for i=2:n
+            r[i] = a[i] - b[i-1]^2/r[i-1]
+        end
+        histogram[BinId] = sum([SturmRatio < 0 for SturmRatio in r])
+    end
+    histogram = int(diff([histogram; n]))
+    return histogram
+end
+