From e91c536e40c56cd12a2fbb06ec86c729253017a8 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Tue, 11 Feb 2014 13:15:15 -0500 Subject: [PATCH] Add another test for 'cholesky'. Improve implementation of upper/lower triangular tests. Drop some typeclass constraints on upper/lower triangular tests. --- src/Linear/Matrix.hs | 91 +++++++++++++++++++++++++++----------------- 1 file changed, 57 insertions(+), 34 deletions(-) diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs index cfa8380..e4acc9a 100644 --- a/src/Linear/Matrix.hs +++ b/src/Linear/Matrix.hs @@ -297,10 +297,34 @@ identity_matrix = -- Examples: -- -- >>> let m1 = fromList [[20,-1], [-1,20]] :: Mat2 Double --- >>> cholesky m1 --- ((4.47213595499958,-0.22360679774997896),(0.0,4.466542286825459)) --- >>> (transpose (cholesky m1)) * (cholesky m1) --- ((20.000000000000004,-1.0),(-1.0,20.0)) +-- >>> let r = cholesky m1 +-- >>> frobenius_norm ((transpose r)*r - m1) < 1e-10 +-- True +-- >>> is_upper_triangular r +-- True +-- +-- >>> import Naturals ( N7 ) +-- >>> let k1 = [6, -3, 0, 0, 0, 0, 0] :: [Double] +-- >>> let k2 = [-3, 10.5, -7.5, 0, 0, 0, 0] :: [Double] +-- >>> let k3 = [0, -7.5, 12.5, 0, 0, 0, 0] :: [Double] +-- >>> let k4 = [0, 0, 0, 6, 0, 0, 0] :: [Double] +-- >>> let k5 = [0, 0, 0, 0, 6, 0, 0] :: [Double] +-- >>> let k6 = [0, 0, 0, 0, 0, 6, 0] :: [Double] +-- >>> let k7 = [0, 0, 0, 0, 0, 0, 15] :: [Double] +-- >>> let big_K = fromList [k1,k2,k3,k4,k5,k6,k7] :: Mat N7 N7 Double +-- +-- >>> let e1 = [2.449489742783178,0,0,0,0,0,0] :: [Double] +-- >>> let e2 = [-1.224744871391589,3,0,0,0,0,0] :: [Double] +-- >>> let e3 = [0,-5/2,5/2,0,0,0,0] :: [Double] +-- >>> let e4 = [0,0,0,2.449489742783178,0,0,0] :: [Double] +-- >>> let e5 = [0,0,0,0,2.449489742783178,0,0] :: [Double] +-- >>> let e6 = [0,0,0,0,0,2.449489742783178,0] :: [Double] +-- >>> let e7 = [0,0,0,0,0,0,3.872983346207417] :: [Double] +-- >>> let expected = fromList [e1,e2,e3,e4,e5,e6,e7] :: Mat N7 N7 Double +-- +-- >>> let r = cholesky big_K +-- >>> frobenius_norm (r - (transpose expected)) < 1e-12 +-- True -- cholesky :: forall m n a. (Algebraic.C a, Arity m, Arity n) => (Mat m n a) -> (Mat m n a) @@ -325,29 +349,26 @@ cholesky m = construct r -- >>> is_upper_triangular' 1e-10 m -- True -- --- TODO: --- --- 1. Don't cheat with lists. --- -is_upper_triangular' :: (Ord a, Ring.C a, Absolute.C a, Arity m, Arity n) +is_upper_triangular' :: forall m n a. + (Ord a, Ring.C a, Absolute.C a, Arity m, Arity n) => a -- ^ The tolerance @epsilon@. -> Mat m n a -> Bool -is_upper_triangular' epsilon m = - and $ concat results +is_upper_triangular' epsilon matrix = + ifoldl2 f True matrix where - results = [[ test i j | i <- [0..(nrows m)-1]] | j <- [0..(ncols m)-1] ] - - test :: Int -> Int -> Bool - test i j + f :: Int -> Int -> Bool -> a -> Bool + f _ _ False _ = False + f i j True x | i <= j = True -- use "less than or equal to" so zero is a valid epsilon - | otherwise = abs (m !!! (i,j)) <= epsilon + | otherwise = abs x <= epsilon -- | Returns True if the given matrix is upper-triangular, and False --- otherwise. A specialized version of 'is_upper_triangular\'' with --- @epsilon = 0@. +-- otherwise. We don't delegate to the general +-- 'is_upper_triangular'' here because it imposes additional +-- typeclass constraints throughout the library. -- -- Examples: -- @@ -359,18 +380,22 @@ is_upper_triangular' epsilon m = -- >>> is_upper_triangular m -- True -- --- TODO: --- --- 1. The Ord constraint is too strong here, Eq would suffice. --- -is_upper_triangular :: (Ord a, Ring.C a, Absolute.C a, Arity m, Arity n) +is_upper_triangular :: forall m n a. + (Eq a, Ring.C a, Arity m, Arity n) => Mat m n a -> Bool -is_upper_triangular = is_upper_triangular' 0 +is_upper_triangular matrix = + ifoldl2 f True matrix + where + f :: Int -> Int -> Bool -> a -> Bool + f _ _ False _ = False + f i j True x + | i <= j = True + | otherwise = x == 0 + -- | Returns True if the given matrix is lower-triangular, and False --- otherwise. This is a specialized version of 'is_lower_triangular\'' --- with @epsilon = 0@. +-- otherwise. -- -- Examples: -- @@ -382,9 +407,8 @@ is_upper_triangular = is_upper_triangular' 0 -- >>> is_lower_triangular m -- False -- -is_lower_triangular :: (Ord a, +is_lower_triangular :: (Eq a, Ring.C a, - Absolute.C a, Arity m, Arity n) => Mat m n a @@ -638,9 +662,9 @@ vec5d (v,w,x,y,z) = Mat (mk5 (mk1 v) (mk1 w) (mk1 x) (mk1 y) (mk1 z)) scalar :: a -> Mat1 a scalar x = Mat (mk1 (mk1 x)) -dot :: (RealRing.C a, n ~ N1, m ~ S t, Arity t) - => Mat m n a - -> Mat m n a +dot :: (RealRing.C a, m ~ S t, Arity t) + => Col m a + -> Col m a -> a v1 `dot` v2 = ((transpose v1) * v2) !!! (0, 0) @@ -656,12 +680,11 @@ v1 `dot` v2 = ((transpose v1) * v2) !!! (0, 0) -- angle :: (Transcendental.C a, RealRing.C a, - n ~ N1, m ~ S t, Arity t, ToRational.C a) - => Mat m n a - -> Mat m n a + => Col m a + -> Col m a -> a angle v1 v2 = acos theta -- 2.44.2