X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FLinear%2FQR.hs;h=8d9ea42a1e19bd880ffb9211238e53ba3d7b8c5b;hb=a96be0c7818e0408b4aac9076df131c5b3ec6ff4;hp=283512a29494da41cb24db69ee2170cb213e53f1;hpb=ae914d13235a4582077a5cb2b1edd630d9c6ad62;p=numerical-analysis.git diff --git a/src/Linear/QR.hs b/src/Linear/QR.hs index 283512a..8d9ea42 100644 --- a/src/Linear/QR.hs +++ b/src/Linear/QR.hs @@ -4,22 +4,28 @@ -- | QR factorization via Givens rotations. -- 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 ( 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 ) @@ -133,7 +139,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 @@ -149,15 +155,134 @@ 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) --- 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\"). +-- +-- 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 ) +-- +-- >>> 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 + -> 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)