From: Michael Orlitzky Date: Tue, 4 Feb 2014 07:19:04 +0000 (-0500) Subject: Implement the QR algorithm for computing eigenvalues. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=47c0f368bd1d6d1b279ea95e3dbf2ccabb879b75;p=numerical-analysis.git Implement the QR algorithm for computing eigenvalues. --- 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)) +