X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FLinear%2FQR.hs;h=043b17256e6af341dbe1fe3a34f18dff649a20a7;hb=2a2db25a6667b2b078390e9ddfadad8c367839ee;hp=79ca5f1a5d394cde0ee1ae469df8f5dba30a5282;hpb=47c0f368bd1d6d1b279ea95e3dbf2ccabb879b75;p=numerical-analysis.git

diff --git a/src/Linear/QR.hs b/src/Linear/QR.hs
index 79ca5f1..043b172 100644
--- a/src/Linear/QR.hs
+++ b/src/Linear/QR.hs
@@ -5,23 +5,27 @@
 --
 module Linear.QR (
   eigenvalues,
+  eigenvectors_symmetric,
   givens_rotator,
   qr )
 where
 
 import qualified Algebra.Ring as Ring ( C )
 import qualified Algebra.Algebraic as Algebraic ( C )
-import Data.Vector.Fixed ( N1, S, ifoldl )
+import Control.Arrow ( first )
+import Data.Vector.Fixed ( S, ifoldl )
 import Data.Vector.Fixed.Cont ( Arity )
 import NumericPrelude hiding ( (*) )
 
 import Linear.Matrix (
+  Col,
   Mat(..),
   (*),
   (!!!),
   construct,
   diagonal,
   identity_matrix,
+  symmetric,
   transpose )
 
 
@@ -151,7 +155,7 @@ qr matrix =
       ifoldl (f col_idx) (q,r) col
 
     -- | Process the entries in a column, doing basically the same
-    --   thing as col_dunction does. It updates the QR factorization,
+    --   thing as col_function 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) -- leave it alone
@@ -167,6 +171,9 @@ qr matrix =
 --   iterated QR algorithm (see Golub and Van Loan, \"Matrix
 --   Computations\").
 --
+--   Warning: this may not converge if there are repeated eigenvalues
+--   (in magnitude).
+--
 --   Examples:
 --
 --   >>> import Linear.Matrix ( Col2, Col3, Mat2, Mat3 )
@@ -193,11 +200,89 @@ qr matrix =
 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))
+             -> Col (S m) a
+eigenvalues iterations matrix
+  | iterations <  0 = error "negative iterations requested"
+  | iterations == 0 = diagonal matrix
+  | otherwise =
+      diagonal (ut_approximation (iterations - 1))
+      where
+        ut_approximation :: Int -> Mat (S m) (S m) a
+        ut_approximation 0 = matrix
+        ut_approximation k = ut_next
+          where
+            ut_prev = ut_approximation (k-1)
+            (qk,rk) = qr ut_prev
+            ut_next = rk*qk
+
+
+
+-- | Compute the eigenvalues and eigenvectors of a symmetric matrix
+--   using an iterative QR algorithm. This is similar to what we do in
+--   'eigenvalues' except we also return the product of all \"Q\"
+--   matrices that we have generated. This turns out to me the matrix
+--   of eigenvectors when the original matrix is symmetric. For
+--   references see Goluv and Van Loan, \"Matrix Computations\", or
+--   \"Calculation of Gauss Quadrature Rules\" by Golub and Welsch.
+--
+--   Warning: this may not converge if there are repeated eigenvalues
+--   (in magnitude).
+--
+--   Examples:
+--
+--   >>> import Linear.Matrix ( Col2, Col3, Mat2, Mat3 )
+--   >>> import Linear.Matrix ( column', frobenius_norm, fromList )
+--   >>> import Linear.Matrix ( identity_matrix, vec3d )
+--   >>> import Normed ( Normed(..) )
+--
+--   >>> let m = identity_matrix :: Mat3 Double
+--   >>> let (vals, vecs) = eigenvectors_symmetric 100 m
+--   >>> let expected_vals = fromList [[1],[1],[1]] :: Col3 Double
+--   >>> let expected_vecs = m
+--   >>> vals == expected_vals
+--   True
+--   >>> vecs == expected_vecs
+--   True
+--
+--   >>> let m = fromList [[3,2,4],[2,0,2],[4,2,3]] :: Mat3 Double
+--   >>> let (vals, vecs) = eigenvectors_symmetric 1000 m
+--   >>> let expected_vals = fromList [[8],[-1],[-1]] :: Col3 Double
+--   >>> let v0' = vec3d (2, 1, 2) :: Col3 Double
+--   >>> let v0  = (1 / (norm v0') :: Double) *> v0'
+--   >>> let v1' = vec3d (-1, 2, 0) :: Col3 Double
+--   >>> let v1  = (1 / (norm v1') :: Double) *> v1'
+--   >>> let v2' = vec3d (-4, -2, 5) :: Col3 Double
+--   >>> let v2  = (1 / (norm v2') :: Double) *> v2'
+--   >>> frobenius_norm ((column' vecs 0) - v0) < 1e-12
+--   True
+--   >>> frobenius_norm ((column' vecs 1) - v1) < 1e-12
+--   True
+--   >>> frobenius_norm ((column' vecs 2) - v2) < 1e-12
+--   True
+--
+eigenvectors_symmetric :: forall m a. (Arity m, Algebraic.C a, Eq a)
+                       => Int
+                       -> Mat (S m) (S m) a
+                       -> (Col (S m) a, Mat (S m) (S m) a)
+eigenvectors_symmetric iterations matrix
+  | iterations <  0 = error "negative iterations requested"
+  | iterations == 0 = (diagonal matrix, identity_matrix)
+  | not $ symmetric matrix = error "argument is not symmetric"
+  | otherwise =
+      (values, vectors)
+      where
+        -- | We think of \"T\" as an approximation to an
+        -- upper-triangular matrix from which we get our
+        -- eigenvalues. The matrix \"P\" is the product of all
+        -- previous \"Q\"s and its columns approximate the
+        -- eigenvectors.
+        tp_pair :: Int -> (Mat (S m) (S m) a, Mat (S m) (S m) a)
+        tp_pair 0 = (matrix, identity_matrix)
+        tp_pair k = (tk,pk)
+                    where
+                      (t_prev, p_prev) = tp_pair (k-1)
+                      (qk,rk) = qr t_prev
+                      pk = p_prev*qk
+                      tk = rk*qk
 
+        (values, vectors) = (first diagonal) (tp_pair iterations)