X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FLinear%2FQR.hs;h=a4fb325f96d4b510401f79d8517b6beb48b8fc79;hb=b32831b5dde3440b85cbef62f4c47fcce0ee974f;hp=f9c5e30294856de0f3f50d7c157067f1e0c9e48f;hpb=8ba752ee2e6a9e4f80843c21b52e2d534fde01c0;p=numerical-analysis.git

diff --git a/src/Linear/QR.hs b/src/Linear/QR.hs
index f9c5e30..a4fb325 100644
--- a/src/Linear/QR.hs
+++ b/src/Linear/QR.hs
@@ -13,7 +13,7 @@ where
 import qualified Algebra.Ring as Ring ( C )
 import qualified Algebra.Algebraic as Algebraic ( C )
 import Control.Arrow ( first )
-import Data.Vector.Fixed ( N1, S, ifoldl )
+import Data.Vector.Fixed ( S, ifoldl )
 import Data.Vector.Fixed.Cont ( Arity )
 import NumericPrelude hiding ( (*) )
 
@@ -140,11 +140,14 @@ givens_rotator i j xi xj =
 --   True
 --
 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)
+   => Mat (S m) (S n) a
+   -> (Mat (S m) (S m) a, Mat (S m) (S n) a)
 qr matrix =
   ifoldl col_function initial_qr columns
   where
     Mat columns = transpose matrix
+
+    initial_qr :: (Mat (S m) (S m) a, Mat (S m) (S n) a)
     initial_qr = (identity_matrix, matrix)
 
     -- | Process the column and spit out the current QR
@@ -155,14 +158,14 @@ 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
       | otherwise = (q*rotator, (transpose rotator)*r)
           where
             y = r !!! (idx, col_idx)
-            rotator :: Mat m m a
+            rotator :: Mat (S m) (S m) a
             rotator = givens_rotator col_idx idx (r !!! (col_idx, col_idx)) y
 
 
@@ -171,6 +174,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 )
@@ -198,12 +204,20 @@ eigenvalues :: forall m a. (Arity m, Algebraic.C a, Eq a)
              => Int
              -> Mat (S m) (S m) a
              -> Col (S m) 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))
+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
@@ -214,10 +228,15 @@ eigenvalues iterations matrix =
 --   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 ( frobenius_norm, fromList, identity_matrix )
+--   >>> 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
@@ -233,11 +252,16 @@ eigenvalues iterations matrix =
 --   >>> 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' = vec3d (-1, 2, 0) :: Col3 Double
 --   >>> let v1  = (1 / (norm v1') :: Double) *> v1'
---   >>> let v2' = vec3d (4, 2, 5) :: Col3 Double
+--   >>> let v2' = vec3d (-4, -2, 5) :: Col3 Double
 --   >>> let v2  = (1 / (norm v2') :: Double) *> v2'
---   >>> frobenius_norm (vals - expected_vals)
+--   >>> 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
@@ -261,8 +285,7 @@ eigenvectors_symmetric iterations matrix
                     where
                       (t_prev, p_prev) = tp_pair (k-1)
                       (qk,rk) = qr t_prev
-                      tk = rk*qk
                       pk = p_prev*qk
-
+                      tk = rk*qk
 
         (values, vectors) = (first diagonal) (tp_pair iterations)