From: Michael Orlitzky Date: Tue, 4 Feb 2014 17:16:52 +0000 (-0500) Subject: Add the eigenvectors_symmetric function and begin writing tests for it. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=8ba752ee2e6a9e4f80843c21b52e2d534fde01c0;p=numerical-analysis.git Add the eigenvectors_symmetric function and begin writing tests for it. --- diff --git a/src/Linear/QR.hs b/src/Linear/QR.hs index 79ca5f1..f9c5e30 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 Control.Arrow ( first ) import Data.Vector.Fixed ( N1, S, ifoldl ) import Data.Vector.Fixed.Cont ( Arity ) import NumericPrelude hiding ( (*) ) import Linear.Matrix ( + Col, Mat(..), (*), (!!!), construct, diagonal, identity_matrix, + symmetric, transpose ) @@ -193,7 +197,7 @@ 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 + -> Col (S m) a eigenvalues iterations matrix = diagonal (ut_approximation iterations) where @@ -201,3 +205,64 @@ eigenvalues iterations matrix = ut_approximation 0 = matrix ut_approximation k = rk*qk where (qk,rk) = qr (ut_approximation (k-1)) + +-- | 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. +-- +-- Examples: +-- +-- >>> import Linear.Matrix ( Col2, Col3, Mat2, Mat3 ) +-- >>> import Linear.Matrix ( frobenius_norm, fromList, identity_matrix ) +-- +-- >>> 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 (vals - expected_vals) +-- +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 + tk = rk*qk + pk = p_prev*qk + + + (values, vectors) = (first diagonal) (tp_pair iterations)