From cb41ccadccd4305065c3576d63d505a5d35c5279 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Wed, 12 Feb 2014 01:52:35 -0500 Subject: [PATCH] Fix the incorrect definition of "minor". Add a "preminor" function that does what "minor" used to. Implement cofactors and a matrix inverse in terms of them. --- src/Linear/Matrix.hs | 85 +++++++++++++++++++++++++++++++++++--------- 1 file changed, 69 insertions(+), 16 deletions(-) diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs index 119d770..2fe5340 100644 --- a/src/Linear/Matrix.hs +++ b/src/Linear/Matrix.hs @@ -56,6 +56,7 @@ import Algebra.Absolute ( abs ) import qualified Algebra.Additive as Additive ( C ) import qualified Algebra.Algebraic as Algebraic ( C ) import Algebra.Algebraic ( root ) +import qualified Algebra.Field as Field ( C ) import qualified Algebra.Ring as Ring ( C ) import qualified Algebra.Module as Module ( C ) import qualified Algebra.RealRing as RealRing ( C ) @@ -73,20 +74,22 @@ type Mat3 a = Mat N3 N3 a type Mat4 a = Mat N4 N4 a type Mat5 a = Mat N5 N5 a +-- * Type synonyms for 1-by-n row "vectors". + -- | Type synonym for row vectors expressed as 1-by-n matrices. type Row n a = Mat N1 n a --- Type synonyms for 1-by-n row "vectors". type Row1 a = Row N1 a type Row2 a = Row N2 a type Row3 a = Row N3 a type Row4 a = Row N4 a type Row5 a = Row N5 a +-- * Type synonyms for n-by-1 column "vectors". + -- | Type synonym for column vectors expressed as n-by-1 matrices. type Col n a = Mat n N1 a --- Type synonyms for n-by-1 column "vectors". type Col1 a = Col N1 a type Col2 a = Col N2 a type Col3 a = Col N3 a @@ -467,35 +470,50 @@ is_triangular :: (Ord a, is_triangular m = is_upper_triangular m || is_lower_triangular m --- | Return the (i,j)th minor of m. +-- | Delete the @i@th row and @j@th column from the matrix. The name +-- \"preminor\" is made up, but is meant to signify that this is +-- usually used in the computationof a minor. A minor is simply the +-- determinant of a preminor in that case. -- -- Examples: -- -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int --- >>> minor m 0 0 :: Mat2 Int +-- >>> preminor m 0 0 :: Mat2 Int -- ((5,6),(8,9)) --- >>> minor m 1 1 :: Mat2 Int +-- >>> preminor m 1 1 :: Mat2 Int -- ((1,3),(7,9)) -- -minor :: (m ~ S r, - n ~ S t, - Arity r, - Arity t) - => Mat m n a +preminor :: (Arity m, Arity n) + => Mat (S m) (S n) a -> Int -> Int - -> Mat r t a -minor (Mat rows) i j = m + -> Mat m n a +preminor (Mat rows) i j = m where rows' = delete rows i m = Mat $ V.map ((flip delete) j) rows' +-- | Compute the i,jth minor of a @matrix@. +-- +-- Examples: +-- +-- >>> let m1 = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Double +-- >>> minor m1 1 1 +-- -12.0 +-- +minor :: (Arity m, Determined (Mat m m) a) + => Mat (S m) (S m) a + -> Int + -> Int + -> a +minor matrix i j = determinant (preminor matrix i j) + class (Eq a, Ring.C a) => Determined p a where determinant :: (p a) -> a instance (Eq a, Ring.C a) => Determined (Mat (S Z) (S Z)) a where - determinant (Mat rows) = (V.head . V.head) rows + determinant = unscalar instance (Ord a, Ring.C a, @@ -517,10 +535,8 @@ instance (Ord a, where m' i j = m !!! (i,j) - det_minor i j = determinant (minor m i j) - determinant_recursive = - sum [ (-1)^(toInteger j) NP.* (m' 0 j) NP.* (det_minor 0 j) + sum [ (-1)^(toInteger j) NP.* (m' 0 j) NP.* (minor m 0 j) | j <- [0..(ncols m)-1] ] @@ -1028,3 +1044,40 @@ set_idx matrix (i,j) newval = if k == i && l == j then newval else existing + + +-- | Compute the i,jth cofactor of the given @matrix@. This simply +-- premultiplues the i,jth minor by (-1)^(i+j). +cofactor :: (Arity m, Determined (Mat m m) a) + => Mat (S m) (S m) a + -> Int + -> Int + -> a +cofactor matrix i j = + (-1)^(toInteger i + toInteger j) NP.* (minor matrix i j) + + +-- | Compute the inverse of a matrix using cofactor expansion +-- (generalized Cramer's rule). +-- +-- Examples: +-- +-- >>> let m1 = fromList [[37,22],[17,54]] :: Mat2 Double +-- >>> let e1 = [54/1624, -22/1624] :: [Double] +-- >>> let e2 = [-17/1624, 37/1624] :: [Double] +-- >>> let expected = fromList [e1, e2] :: Mat2 Double +-- >>> let actual = inverse m1 +-- >>> frobenius_norm (actual - expected) < 1e-12 +-- True +-- +inverse :: (Arity m, + Determined (Mat (S m) (S m)) a, + Determined (Mat m m) a, + Field.C a) + => Mat (S m) (S m) a + -> Mat (S m) (S m) a +inverse matrix = + (1 / (determinant matrix)) *> (transpose $ construct lambda) + where + lambda i j = cofactor matrix i j + -- 2.44.2