diff --git a/Data/Eigen/Internal.hs b/Data/Eigen/Internal.hs --- a/Data/Eigen/Internal.hs +++ b/Data/Eigen/Internal.hs @@ -54,3 +54,9 @@ foreign import ccall "eigen-proxy.h eigen_hypotNorm" c_hypotNorm :: Ptr CDouble -> CInt -> CInt -> IO CDouble foreign import ccall "eigen-proxy.h eigen_determinant" c_determinant :: Ptr CDouble -> CInt -> CInt -> IO CDouble +foreign import ccall "eigen-proxy.h eigen_rank" c_rank :: CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString +foreign import ccall "eigen-proxy.h eigen_image" c_image :: CInt -> Ptr (Ptr CDouble) -> Ptr CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString +foreign import ccall "eigen-proxy.h eigen_kernel" c_kernel :: CInt -> Ptr (Ptr CDouble) -> Ptr CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString +foreign import ccall "eigen-proxy.h free" c_free :: Ptr a -> IO () +foreign import ccall "eigen-proxy.h eigen_solve" c_solve :: CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString +foreign import ccall "eigen-proxy.h eigen_relativeError" c_relativeError :: Ptr CDouble -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString diff --git a/Data/Eigen/LA.hs b/Data/Eigen/LA.hs --- a/Data/Eigen/LA.hs +++ b/Data/Eigen/LA.hs @@ -48,25 +48,63 @@ Checking if a solution really exists: Only you know what error margin you want to allow for a solution to be considered valid. -You can compute relative error using @norm (ax - b) / norm b@ formula or use 'relativeError' function which provides the same calculation implemented slightly more efficient. +You can compute relative error using @'norm' (ax - b) / 'norm' b@ formula or use 'relativeError' function which provides the same calculation implemented slightly more efficient. -} module Data.Eigen.LA ( -- * Basic linear solving Decomposition(..), + solve, + relativeError, + -- * Rank-revealing decompositions + {- | +Certain decompositions are rank-revealing, i.e. are able to compute the 'rank' of a matrix. These are typically also the decompositions that behave best in the face of a non-full-rank matrix (which in the 'square' case means a singular matrix). + +@ +import Data.Eigen.Matrix +import Data.Eigen.LA + +main = do + let a = fromList [[1,2,5],[2,1,4],[3,0,3]] + putStrLn "Here is the matrix A:" >> print a + putStrLn "The rank of A is:" >> print (rank FullPivLU a) + putStrLn "Here is a matrix whose columns form a basis of the null-space of A:" >> print (kernel FullPivLU a) + putStrLn "Here is a matrix whose columns form a basis of the column-space of A:" >> print (image FullPivLU a) +@ + +produces the following output + +@ +Here is the matrix A: +Matrix 3x3 +1.0 2.0 5.0 +2.0 1.0 4.0 +3.0 0.0 3.0 + +The rank of A is: +2 +Here is a matrix whose columns form a basis of the null-space of A: +Matrix 3x1 +0.5000000000000001 +1.0 +-0.5 + +Here is a matrix whose columns form a basis of the column-space of A: +Matrix 3x2 +5.0 1.0 +4.0 2.0 +3.0 3.0 +@ + -} rank, kernel, image, - solve, - relativeError, -- * Multiple linear regression + {- | A linear regression model that contains more than one predictor variable. -} linearRegression ) where -import Foreign.Ptr -import Foreign.C.Types -import Foreign.C.String import Foreign.Storable import Foreign.Marshal.Alloc import qualified Foreign.Concurrent as FC @@ -76,15 +114,6 @@ import qualified Data.Eigen.Matrix.Mutable as M import qualified Data.Vector.Storable as VS -foreign import ccall "eigen-proxy.h eigen_solve" c_solve :: CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString -foreign import ccall "eigen-proxy.h eigen_relativeError" c_relativeError :: Ptr CDouble -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString -foreign import ccall "eigen-proxy.h eigen_rank" c_rank :: CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString -foreign import ccall "eigen-proxy.h eigen_image" c_image :: CInt -> Ptr (Ptr CDouble) -> Ptr CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString -foreign import ccall "eigen-proxy.h eigen_kernel" c_kernel :: CInt -> Ptr (Ptr CDouble) -> Ptr CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString -foreign import ccall "eigen-proxy.h free" c_free :: Ptr a -> IO () - - - {- | @ Decomposition Requirements on the matrix Speed Accuracy Rank Kernel Image @@ -97,9 +126,8 @@ LLT Positive definite +++ + - - - LDLT Positive or negative semidefinite +++ ++ - - - JacobiSVD None - +++ + - - - -The best way to do least squares solving for square matrices is with a SVD decomposition (JacobiSVD) @ +The best way to do least squares solving for square matrices is with a SVD decomposition ('JacobiSVD') -} data Decomposition @@ -120,12 +148,38 @@ -- | Two-sided Jacobi SVD decomposition of a rectangular matrix. | JacobiSVD deriving (Show, Enum) + +-- | [x = solve d a b] finds a solution @x@ of @ax = b@ equation using decomposition @d@ +solve :: Decomposition -> Matrix -> Matrix -> Matrix +solve d a b = I.performIO $ do + x <- M.new (m_cols a) 1 + M.unsafeWith x $ \x_vals x_rows x_cols -> + unsafeWith a $ \a_vals a_rows a_cols -> + unsafeWith b $ \b_vals b_rows b_cols -> + I.call $ I.c_solve (I.cast $ fromEnum d) + x_vals x_rows x_cols + a_vals a_rows a_cols + b_vals b_rows b_cols + unsafeFreeze x + +-- | [e = relativeError x a b] computes @norm (ax - b) / norm b@ where @norm@ is L2 norm +relativeError :: Matrix -> Matrix -> Matrix -> Double +relativeError x a b = I.performIO $ + unsafeWith x $ \x_vals x_rows x_cols -> + unsafeWith a $ \a_vals a_rows a_cols -> + unsafeWith b $ \b_vals b_rows b_cols -> + alloca $ \pe -> do + I.call $ I.c_relativeError pe + x_vals x_rows x_cols + a_vals a_rows a_cols + b_vals b_rows b_cols + I.cast <$> peek pe + +-- | The rank of the matrix rank :: Decomposition -> Matrix -> Int -rank d m = I.performIO $ - unsafeWith m $ \vals rows cols -> - alloca $ \pr -> do - I.call $ c_rank (I.cast $ fromEnum d) pr vals rows cols - I.cast <$> peek pr +rank d m = I.performIO $ alloca $ \pr -> do + I.call $ unsafeWith m $ I.c_rank (I.cast $ fromEnum d) pr + I.cast <$> peek pr -- | Return matrix whose columns form a basis of the null-space of @A@ kernel :: Decomposition -> Matrix -> Matrix @@ -134,13 +188,13 @@ alloca $ \prows -> alloca $ \pcols -> unsafeWith m1 $ \vals1 rows1 cols1 -> do - I.call $ c_kernel (I.cast $ fromEnum d) + I.call $ I.c_kernel (I.cast $ fromEnum d) pvals prows pcols vals1 rows1 cols1 vals <- peek pvals rows <- I.cast <$> peek prows cols <- I.cast <$> peek pcols - fp <- FC.newForeignPtr vals $ c_free vals + fp <- FC.newForeignPtr vals $ I.c_free vals return $ Matrix rows cols $ VS.unsafeFromForeignPtr0 fp $ rows * cols @@ -151,41 +205,16 @@ alloca $ \prows -> alloca $ \pcols -> unsafeWith m1 $ \vals1 rows1 cols1 -> do - I.call $ c_image (I.cast $ fromEnum d) + I.call $ I.c_image (I.cast $ fromEnum d) pvals prows pcols vals1 rows1 cols1 vals <- peek pvals rows <- I.cast <$> peek prows cols <- I.cast <$> peek pcols - fp <- FC.newForeignPtr vals $ c_free vals + fp <- FC.newForeignPtr vals $ I.c_free vals return $ Matrix rows cols $ VS.unsafeFromForeignPtr0 fp $ rows * cols --- | [x = solve d a b] finds a solution @x@ of @ax = b@ equation using decomposition @d@ -solve :: Decomposition -> Matrix -> Matrix -> Matrix -solve d a b = I.performIO $ do - x <- M.new (m_cols a) 1 - M.unsafeWith x $ \x_vals x_rows x_cols -> - unsafeWith a $ \a_vals a_rows a_cols -> - unsafeWith b $ \b_vals b_rows b_cols -> - I.call $ c_solve (I.cast $ fromEnum d) - x_vals x_rows x_cols - a_vals a_rows a_cols - b_vals b_rows b_cols - unsafeFreeze x --- | [e = relativeError x a b] computes @norm (ax - b) / norm b@ where @norm@ is L2 norm -relativeError :: Matrix -> Matrix -> Matrix -> Double -relativeError x a b = I.performIO $ - unsafeWith x $ \x_vals x_rows x_cols -> - unsafeWith a $ \a_vals a_rows a_cols -> - unsafeWith b $ \b_vals b_rows b_cols -> - alloca $ \pe -> do - I.call $ c_relativeError pe - x_vals x_rows x_cols - a_vals a_rows a_cols - b_vals b_rows b_cols - I.cast <$> peek pe - {- | [(coeffs, error) = linearRegression points] computes multiple linear regression @y = a1 x1 + a2 x2 + ... + an xn + b@ using 'ColPivHouseholderQR' decomposition @@ -203,7 +232,7 @@ [-4.01, 2.41, 6.01], [-3.86, 2.09, 5.55], [-4.10, 2.58, 6.32]] - @ +@ produces the following output diff --git a/Data/Eigen/Matrix.hs b/Data/Eigen/Matrix.hs --- a/Data/Eigen/Matrix.hs +++ b/Data/Eigen/Matrix.hs @@ -253,11 +253,11 @@ trace :: Matrix -> Double trace = _prop I.c_trace --- | Returns true if all of the coefficients in a given matrix evaluate to true +-- | Applied to a predicate and a matrix, all determines if all elements of the matrix satisfies the predicate all :: (Double -> Bool) -> Matrix -> Bool all f = VS.all (f . I.cast) . m_vals --- | Returns true if at least one of the coefficients in a given matrix evaluates to true +-- | Applied to a predicate and a matrix, any determines if any element of the matrix satisfies the predicate any :: (Double -> Bool) -> Matrix -> Bool any f = VS.any (f . I.cast) . m_vals @@ -287,22 +287,19 @@ | square m = _prop I.c_determinant m | otherwise = error "Matrix.determinant: non-square matrix" --- | Adding two matrices by adding the corresponding entries together +-- | Adding two matrices by adding the corresponding entries together. You can use @(+)@ function as well. add :: Matrix -> Matrix -> Matrix add m1 m2 | dims m1 == dims m2 = _binop const I.c_add m1 m2 | otherwise = error "Matrix.add: matrix should have the same size" --- | Return a + b +-- | Subtracting two matrices by subtracting the corresponding entries together. You can use @(-)@ function as well. sub :: Matrix -> Matrix -> Matrix sub m1 m2 | dims m1 == dims m2 = _binop const I.c_sub m1 m2 | otherwise = error "Matrix.add: matrix should have the same size" -{- | Matrix multiplication - -<> --} +-- | Matrix multiplication. You can use @(*)@ function as well. mul :: Matrix -> Matrix -> Matrix mul m1 m2 | m_cols m1 == m_rows m2 = _binop (\(rows, _) (_, cols) -> (rows, cols)) I.c_mul m1 m2 diff --git a/eigen.cabal b/eigen.cabal --- a/eigen.cabal +++ b/eigen.cabal @@ -1,5 +1,5 @@ name: eigen -version: 1.2.2 +version: 1.2.3 homepage: https://github.com/osidorkin/haskell-eigen synopsis: Eigen C++ library (linear algebra: matrices, vectors, numerical solvers). description: This module provides Haskell binding for Eigen C++ library. @@ -365,7 +365,8 @@ Data.Eigen.Matrix Data.Eigen.Matrix.Mutable Data.Eigen.Parallel - Data.Eigen.Internal + + other-modules: Data.Eigen.Internal ghc-options: -Wall -fno-warn-name-shadowing build-depends: base >= 3 && < 5,