From 47c0f368bd1d6d1b279ea95e3dbf2ccabb879b75 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Tue, 4 Feb 2014 02:19:04 -0500
Subject: [PATCH] Implement the QR algorithm for computing eigenvalues.

---
 src/Linear/QR.hs | 52 ++++++++++++++++++++++++++++++++++++++++++------
 1 file changed, 46 insertions(+), 6 deletions(-)

diff --git a/src/Linear/QR.hs b/src/Linear/QR.hs
index 283512a..79ca5f1 100644
--- a/src/Linear/QR.hs
+++ b/src/Linear/QR.hs
@@ -4,13 +4,14 @@
 -- | QR factorization via Givens rotations.
 --
 module Linear.QR (
+  eigenvalues,
   givens_rotator,
   qr )
 where
 
 import qualified Algebra.Ring as Ring ( C )
 import qualified Algebra.Algebraic as Algebraic ( C )
-import Data.Vector.Fixed ( ifoldl )
+import Data.Vector.Fixed ( N1, S, ifoldl )
 import Data.Vector.Fixed.Cont ( Arity )
 import NumericPrelude hiding ( (*) )
 
@@ -19,6 +20,7 @@ import Linear.Matrix (
   (*),
   (!!!),
   construct,
+  diagonal,
   identity_matrix,
   transpose )
 
@@ -133,7 +135,7 @@ givens_rotator i j xi xj =
 --   >>> is_upper_triangular' 1e-10 r
 --   True
 --
-qr :: forall m n a. (Arity m, Arity n, Eq a, Algebraic.C a, Ring.C a, Show a)
+qr :: forall m n a. (Arity m, Arity n, Eq a, Algebraic.C a, Ring.C a)
    => Mat m n a -> (Mat m m a, Mat m n a)
 qr matrix =
   ifoldl col_function initial_qr columns
@@ -152,12 +154,50 @@ qr matrix =
     --   thing as col_dunction does. It updates the QR factorization,
     --   maybe, and returns the current one.
     f col_idx (q,r) idx _ -- ignore the current element
-      | idx <= col_idx = (q,r)
---          trace ("---------------\nidx: " ++ (show idx) ++ ";\ncol_idx: " ++ (show col_idx) ++ "; leaving it alone") (q,r) -- leave it alone.
+      | idx <= col_idx = (q,r) -- leave it alone
       | otherwise = (q*rotator, (transpose rotator)*r)
---          trace ("---------------\nidx: " ++ (show idx) ++ ";\ncol_idx: " ++ (show col_idx) ++ ";\nupdating Q and R;\nq: " ++ (show q) ++ ";\nr " ++ (show r) ++ ";\nnew q: " ++ (show $ q*rotator) ++ ";\nnew r: " ++ (show $ (transpose rotator)*r) ++ ";\ny: " ++ (show y) ++ ";\nr[i,j]: " ++ (show (r !!! (col_idx, col_idx))))
---          (q*rotator, (transpose rotator)*r)
           where
             y = r !!! (idx, col_idx)
             rotator :: Mat m m a
             rotator = givens_rotator col_idx idx (r !!! (col_idx, col_idx)) y
+
+
+
+-- | Determine the eigenvalues of the given @matrix@ using the
+--   iterated QR algorithm (see Golub and Van Loan, \"Matrix
+--   Computations\").
+--
+--   Examples:
+--
+--   >>> import Linear.Matrix ( Col2, Col3, Mat2, Mat3 )
+--   >>> import Linear.Matrix ( frobenius_norm, fromList, identity_matrix )
+--
+--   >>> let m = fromList [[1,1],[-2,4]] :: Mat2 Double
+--   >>> let actual = eigenvalues 1000 m
+--   >>> let expected = fromList [[3],[2]] :: Col2 Double
+--   >>> frobenius_norm (actual - expected) < 1e-12
+--   True
+--
+--   >>> let m = identity_matrix :: Mat2 Double
+--   >>> let actual = eigenvalues 10 m
+--   >>> let expected = fromList [[1],[1]] :: Col2 Double
+--   >>> frobenius_norm (actual - expected) < 1e-12
+--   True
+--
+--   >>> let m = fromList [[0,1,0],[0,0,1],[1,-3,3]] :: Mat3 Double
+--   >>> let actual = eigenvalues 1000 m
+--   >>> let expected = fromList [[1],[1],[1]] :: Col3 Double
+--   >>> frobenius_norm (actual - expected) < 1e-2
+--   True
+--
+eigenvalues :: forall m a. (Arity m, Algebraic.C a, Eq a)
+             => Int
+             -> Mat (S m) (S m) a
+             -> Mat (S m) N1 a
+eigenvalues iterations matrix =
+  diagonal (ut_approximation iterations)
+  where
+    ut_approximation :: Int -> Mat (S m) (S m) a
+    ut_approximation 0 = matrix
+    ut_approximation k = rk*qk where (qk,rk) = qr (ut_approximation (k-1))
+
-- 
2.43.2