Add tolerant versions of is_{upper,lower}_triangular.

[numerical-analysis.git] / src / Linear / Matrix.hs

diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs
index d166e44367e626b9d87979a4b0e9ab7b2b24e7d9..c6f4a83c71af1c5231113ec555fa79f2dcfcc038 100644 (file)
--- a/src/Linear/Matrix.hs
+++ b/src/Linear/Matrix.hs
@@ -2,6 +2,7 @@
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}

+{-# LANGUAGE NoMonomorphismRestriction #-}

 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE RebindableSyntax #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE RebindableSyntax #-}


@@ -25,6 +26,7 @@ import Data.Vector.Fixed (
   N5,
   S,
   Z,
   N5,
   S,
   Z,

+  generate,

   mk1,
   mk2,
   mk3,
   mk1,
   mk2,
   mk3,


@@ -42,16 +44,17 @@ import qualified Data.Vector.Fixed as V (
   toList,
   zipWith
   )
   toList,
   zipWith
   )

-import Data.Vector.Fixed.Boxed (Vec)

 import Data.Vector.Fixed.Cont (Arity, arity)
 import Linear.Vector
 import Normed
 
 import Data.Vector.Fixed.Cont (Arity, arity)
 import Linear.Vector
 import Normed
 

-import NumericPrelude hiding ((*), abs)
-import qualified NumericPrelude as NP ((*))
+import NumericPrelude hiding ( (*), abs )
+import qualified NumericPrelude as NP ( (*) )
+import qualified Algebra.Absolute as Absolute ( C )
+import Algebra.Absolute ( abs )
+import qualified Algebra.Additive as Additive

 import qualified Algebra.Algebraic as Algebraic
 import Algebra.Algebraic (root)
 import qualified Algebra.Algebraic as Algebraic
 import Algebra.Algebraic (root)

-import qualified Algebra.Additive as Additive

 import qualified Algebra.Ring as Ring
 import qualified Algebra.Module as Module
 import qualified Algebra.RealRing as RealRing
 import qualified Algebra.Ring as Ring
 import qualified Algebra.Module as Module
 import qualified Algebra.RealRing as RealRing


@@ -198,8 +201,6 @@ symmetric m =
 --   entries in the matrix. The i,j entry of the resulting matrix will
 --   have the value returned by lambda i j.
 --
 --   entries in the matrix. The i,j entry of the resulting matrix will
 --   have the value returned by lambda i j.
 --

---   TODO: Don't cheat with fromList.
---

 --   Examples:
 --
 --   >>> let lambda i j = i + j
 --   Examples:
 --
 --   >>> let lambda i j = i + j


@@ -208,15 +209,25 @@ symmetric m =
 --
 construct :: forall m n a. (Arity m, Arity n)
           => (Int -> Int -> a) -> Mat m n a
 --
 construct :: forall m n a. (Arity m, Arity n)
           => (Int -> Int -> a) -> Mat m n a

-construct lambda = Mat rows
+construct lambda = Mat $ generate make_row

   where
   where

-    -- The arity trick is used in Data.Vector.Fixed.length.
-    imax = (arity (undefined :: m)) - 1
-    jmax = (arity (undefined :: n)) - 1
-    row' i = V.fromList [ lambda i j | j <- [0..jmax] ]
-    rows = V.fromList [ row' i | i <- [0..imax] ]
+    make_row :: Int -> Vec n a
+    make_row i = generate (lambda i)

 
 
 
 

+-- | Create an identity matrix with the right dimensions.
+--
+--   Examples:
+--
+--   >>> identity_matrix :: Mat3 Int
+--   ((1,0,0),(0,1,0),(0,0,1))
+--   >>> identity_matrix :: Mat3 Double
+--   ((1.0,0.0,0.0),(0.0,1.0,0.0),(0.0,0.0,1.0))
+--
+identity_matrix :: (Arity m, Ring.C a) => Mat m m a
+identity_matrix =
+  construct (\i j -> if i == j then (fromInteger 1) else (fromInteger 0))
+

 -- | Given a positive-definite matrix @m@, computes the
 --   upper-triangular matrix @r@ with (transpose r)*r == m and all
 --   values on the diagonal of @r@ positive.
 -- | Given a positive-definite matrix @m@, computes the
 --   upper-triangular matrix @r@ with (transpose r)*r == m and all
 --   values on the diagonal of @r@ positive.


@@ -241,21 +252,26 @@ cholesky m = construct r
 
 
 -- | Returns True if the given matrix is upper-triangular, and False
 
 
 -- | Returns True if the given matrix is upper-triangular, and False

---   otherwise.
+--   otherwise. The parameter @epsilon@ lets the caller choose a
+--   tolerance.

 --
 --   Examples:
 --
 --
 --   Examples:
 --

---   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
+--   >>> let m = fromList [[1,1],[1e-12,1]] :: Mat2 Double

 --   >>> is_upper_triangular m
 --   False
 --   >>> is_upper_triangular m
 --   False

---
---   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
---   >>> is_upper_triangular m
+--   >>> is_upper_triangular' 1e-10 m

 --   True
 --
 --   True
 --

-is_upper_triangular :: (Eq a, Ring.C a, Arity m, Arity n)
-                    => Mat m n a -> Bool
-is_upper_triangular m =
+--   TODO:
+--
+--     1. Don't cheat with lists.
+--
+is_upper_triangular' :: (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
   where
     results = [[ test i j | i <- [0..(nrows m)-1]] | j <- [0..(ncols m)-1] ]
   and $ concat results
   where
     results = [[ test i j | i <- [0..(nrows m)-1]] | j <- [0..(ncols m)-1] ]


@@ -263,11 +279,36 @@ is_upper_triangular m =
     test :: Int -> Int -> Bool
     test i j
       | i <= j = True
     test :: Int -> Int -> Bool
     test i j
       | i <= j = True

-      | otherwise = m !!! (i,j) == 0
+      -- use "less than or equal to" so zero is a valid epsilon
+      | otherwise = abs (m !!! (i,j)) <= epsilon
+
+
+-- | Returns True if the given matrix is upper-triangular, and False
+--   otherwise. A specialized version of 'is_upper_triangular\'' with
+--   @epsilon = 0@.
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
+--   >>> is_upper_triangular m
+--   False
+--
+--   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
+--   >>> 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)
+                    => Mat m n a -> Bool
+is_upper_triangular = is_upper_triangular' 0

 
 
 -- | Returns True if the given matrix is lower-triangular, and False
 
 
 -- | Returns True if the given matrix is lower-triangular, and False

---   otherwise.
+--   otherwise. This is a specialized version of 'is_lower_triangular\''
+--   with @epsilon = 0@.

 --
 --   Examples:
 --
 --
 --   Examples:
 --


@@ -279,8 +320,9 @@ is_upper_triangular m =
 --   >>> is_lower_triangular m
 --   False
 --
 --   >>> is_lower_triangular m
 --   False
 --

-is_lower_triangular :: (Eq a,
+is_lower_triangular :: (Ord a,

                         Ring.C a,
                         Ring.C a,

+                        Absolute.C a,

                         Arity m,
                         Arity n)
                     => Mat m n a
                         Arity m,
                         Arity n)
                     => Mat m n a


@@ -288,6 +330,29 @@ is_lower_triangular :: (Eq a,
 is_lower_triangular = is_upper_triangular . transpose
 
 
 is_lower_triangular = is_upper_triangular . transpose
 
 

+-- | Returns True if the given matrix is lower-triangular, and False
+--   otherwise. The parameter @epsilon@ lets the caller choose a
+--   tolerance.
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1,1e-12],[1,1]] :: Mat2 Double
+--   >>> is_lower_triangular m
+--   False
+--   >>> is_lower_triangular' 1e-12 m
+--   True
+--
+is_lower_triangular' :: (Ord a,
+                         Ring.C a,
+                         Absolute.C a,
+                         Arity m,
+                         Arity n)
+                    => a -- ^ The tolerance @epsilon@.
+                    -> Mat m n a
+                    -> Bool
+is_lower_triangular' epsilon = (is_upper_triangular' epsilon) . transpose
+
+

 -- | Returns True if the given matrix is triangular, and False
 --   otherwise.
 --
 -- | Returns True if the given matrix is triangular, and False
 --   otherwise.
 --


@@ -305,8 +370,9 @@ is_lower_triangular = is_upper_triangular . transpose
 --   >>> is_triangular m
 --   False
 --
 --   >>> is_triangular m
 --   False
 --

-is_triangular :: (Eq a,
+is_triangular :: (Ord a,

                   Ring.C a,
                   Ring.C a,

+                  Absolute.C a,

                   Arity m,
                   Arity n)
               => Mat m n a
                   Arity m,
                   Arity n)
               => Mat m n a


@@ -344,8 +410,9 @@ class (Eq a, Ring.C a) => Determined p a where
 instance (Eq a, Ring.C a) => Determined (Mat (S Z) (S Z)) a where
   determinant (Mat rows) = (V.head . V.head) rows
 
 instance (Eq a, Ring.C a) => Determined (Mat (S Z) (S Z)) a where
   determinant (Mat rows) = (V.head . V.head) rows
 

-instance (Eq a,
+instance (Ord a,

           Ring.C a,
           Ring.C a,

+          Absolute.C a,

           Arity n,
           Determined (Mat (S n) (S n)) a)
          => Determined (Mat (S (S n)) (S (S n))) a where
           Arity n,
           Determined (Mat (S n) (S n)) a)
          => Determined (Mat (S (S n)) (S (S n))) a where


@@ -419,8 +486,8 @@ instance (Algebraic.C a,
           ToRational.C a,
           Arity m)
          => Normed (Mat (S m) N1 a) where
           ToRational.C a,
           Arity m)
          => Normed (Mat (S m) N1 a) where

-  -- | Generic p-norms. The usual norm in R^n is (norm_p 2). We treat
-  --   all matrices as big vectors.
+  -- | Generic p-norms for vectors in R^n that are represented as nx1
+  --   matrices.

   --
   --   Examples:
   --
   --
   --   Examples:
   --


@@ -448,7 +515,26 @@ instance (Algebraic.C a,
     fromRational' $ toRational $ V.maximum $ V.map V.maximum rows
 
 
     fromRational' $ toRational $ V.maximum $ V.map V.maximum rows
 
 

-
+-- | Compute the Frobenius norm of a matrix. This essentially treats
+--   the matrix as one long vector containing all of its entries (in
+--   any order, it doesn't matter).
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1, 2, 3],[4,5,6],[7,8,9]] :: Mat3 Double
+--   >>> frobenius_norm m == sqrt 285
+--   True
+--
+--   >>> let m = fromList [[1, -1, 1],[-1,1,-1],[1,-1,1]] :: Mat3 Double
+--   >>> frobenius_norm m == 3
+--   True
+--
+frobenius_norm :: (Algebraic.C a, Ring.C a) => Mat m n a -> a
+frobenius_norm (Mat rows) =
+  sqrt $ element_sum $ V.map row_sum rows
+  where
+    -- | Square and add up the entries of a row.
+    row_sum = element_sum . V.map (^2)

 
 
 -- Vector helpers. We want it to be easy to create low-dimension
 
 
 -- Vector helpers. We want it to be easy to create low-dimension