Safe Haskell | None |
---|---|
Language | Haskell98 |
Matrix datatype and operations.
Every provided example has been tested.
Run cabal test
for further tests.
- data Matrix a
- prettyMatrix :: Show a => Matrix a -> String
- nrows :: Matrix a -> Int
- ncols :: Matrix a -> Int
- forceMatrix :: Matrix a -> Matrix a
- matrix :: Int -> Int -> ((Int, Int) -> a) -> Matrix a
- rowVector :: Vector a -> Matrix a
- colVector :: Vector a -> Matrix a
- zero :: Num a => Int -> Int -> Matrix a
- identity :: Num a => Int -> Matrix a
- diagonalList :: Int -> a -> [a] -> Matrix a
- diagonal :: a -> Vector a -> Matrix a
- permMatrix :: Num a => Int -> Int -> Int -> Matrix a
- fromList :: Int -> Int -> [a] -> Matrix a
- fromLists :: [[a]] -> Matrix a
- toList :: Matrix a -> [a]
- toLists :: Matrix a -> [[a]]
- getElem :: Int -> Int -> Matrix a -> a
- (!) :: Matrix a -> (Int, Int) -> a
- unsafeGet :: Int -> Int -> Matrix a -> a
- safeGet :: Int -> Int -> Matrix a -> Maybe a
- safeSet :: a -> (Int, Int) -> Matrix a -> Maybe (Matrix a)
- getRow :: Int -> Matrix a -> Vector a
- getCol :: Int -> Matrix a -> Vector a
- getDiag :: Matrix a -> Vector a
- getMatrixAsVector :: Matrix a -> Vector a
- setElem :: a -> (Int, Int) -> Matrix a -> Matrix a
- unsafeSet :: a -> (Int, Int) -> Matrix a -> Matrix a
- transpose :: Matrix a -> Matrix a
- setSize :: a -> Int -> Int -> Matrix a -> Matrix a
- extendTo :: a -> Int -> Int -> Matrix a -> Matrix a
- inverse :: (Fractional a, Eq a) => Matrix a -> Either String (Matrix a)
- rref :: (Fractional a, Eq a) => Matrix a -> Either String (Matrix a)
- mapRow :: (Int -> a -> a) -> Int -> Matrix a -> Matrix a
- mapCol :: (Int -> a -> a) -> Int -> Matrix a -> Matrix a
- submatrix :: Int -> Int -> Int -> Int -> Matrix a -> Matrix a
- minorMatrix :: Int -> Int -> Matrix a -> Matrix a
- splitBlocks :: Int -> Int -> Matrix a -> (Matrix a, Matrix a, Matrix a, Matrix a)
- (<|>) :: Matrix a -> Matrix a -> Matrix a
- (<->) :: Matrix a -> Matrix a -> Matrix a
- joinBlocks :: (Matrix a, Matrix a, Matrix a, Matrix a) -> Matrix a
- elementwise :: (a -> b -> c) -> Matrix a -> Matrix b -> Matrix c
- elementwiseUnsafe :: (a -> b -> c) -> Matrix a -> Matrix b -> Matrix c
- multStd :: Num a => Matrix a -> Matrix a -> Matrix a
- multStd2 :: Num a => Matrix a -> Matrix a -> Matrix a
- multStrassen :: Num a => Matrix a -> Matrix a -> Matrix a
- multStrassenMixed :: Num a => Matrix a -> Matrix a -> Matrix a
- scaleMatrix :: Num a => a -> Matrix a -> Matrix a
- scaleRow :: Num a => a -> Int -> Matrix a -> Matrix a
- combineRows :: Num a => Int -> a -> Int -> Matrix a -> Matrix a
- switchRows :: Int -> Int -> Matrix a -> Matrix a
- switchCols :: Int -> Int -> Matrix a -> Matrix a
- luDecomp :: (Ord a, Fractional a) => Matrix a -> Maybe (Matrix a, Matrix a, Matrix a, a)
- luDecompUnsafe :: (Ord a, Fractional a) => Matrix a -> (Matrix a, Matrix a, Matrix a, a)
- luDecomp' :: (Ord a, Fractional a) => Matrix a -> Maybe (Matrix a, Matrix a, Matrix a, Matrix a, a, a)
- luDecompUnsafe' :: (Ord a, Fractional a) => Matrix a -> (Matrix a, Matrix a, Matrix a, Matrix a, a, a)
- cholDecomp :: Floating a => Matrix a -> Matrix a
- trace :: Num a => Matrix a -> a
- diagProd :: Num a => Matrix a -> a
- detLaplace :: Num a => Matrix a -> a
- detLU :: (Ord a, Fractional a) => Matrix a -> a
- flatten :: Matrix (Matrix a) -> Matrix a
Matrix type
Type of matrices.
Elements can be of any type. Rows and columns
are indexed starting by 1. This means that, if m :: Matrix a
and
i,j :: Int
, then m ! (i,j)
is the element in the i
-th row and
j
-th column of m
.
prettyMatrix :: Show a => Matrix a -> String Source
forceMatrix :: Matrix a -> Matrix a Source
Builders
O(rows*cols). Generate a matrix from a generator function. Example of usage:
( 1 0 -1 -2 ) ( 3 2 1 0 ) ( 5 4 3 2 ) matrix 4 4 $ \(i,j) -> 2*i - j = ( 7 6 5 4 )
Special matrices
O(rows*cols). The zero matrix of the given size.
zero n m = m 1 ( 0 0 ... 0 0 ) 2 ( 0 0 ... 0 0 ) ( ... ) ( 0 0 ... 0 0 ) n ( 0 0 ... 0 0 )
identity :: Num a => Int -> Matrix a Source
O(rows*cols). Identity matrix of the given order.
identity n = n 1 ( 1 0 ... 0 0 ) 2 ( 0 1 ... 0 0 ) ( ... ) ( 0 0 ... 1 0 ) n ( 0 0 ... 0 1 )
diagonalList :: Int -> a -> [a] -> Matrix a Source
Diagonal matrix from a non-empty list given the desired size. Non-diagonal elements will be filled with the given default element. The list must have at least order elements.
diagonalList n 0 [1..] = n 1 ( 1 0 ... 0 0 ) 2 ( 0 2 ... 0 0 ) ( ... ) ( 0 0 ... n-1 0 ) n ( 0 0 ... 0 n )
Similar to diagonalList
, but using Vector
, which
should be more efficient.
:: Num a | |
=> Int | Size of the matrix. |
-> Int | Permuted row 1. |
-> Int | Permuted row 2. |
-> Matrix a | Permutation matrix. |
O(rows*cols). Permutation matrix.
permMatrix n i j = i j n 1 ( 1 0 ... 0 ... 0 ... 0 0 ) 2 ( 0 1 ... 0 ... 0 ... 0 0 ) ( ... ... ... ) i ( 0 0 ... 0 ... 1 ... 0 0 ) ( ... ... ... ) j ( 0 0 ... 1 ... 0 ... 0 0 ) ( ... ... ... ) ( 0 0 ... 0 ... 0 ... 1 0 ) n ( 0 0 ... 0 ... 0 ... 0 1 )
When i == j
it reduces to identity
n
.
List conversions
Create a matrix from a non-empty list given the desired size. The list must have at least rows*cols elements. An example:
( 1 2 3 ) ( 4 5 6 ) fromList 3 3 [1..] = ( 7 8 9 )
fromLists :: [[a]] -> Matrix a Source
Create a matrix from a non-empty list of non-empty lists. Each list must have at least as many elements as the first list. Examples:
fromLists [ [1,2,3] ( 1 2 3 ) , [4,5,6] ( 4 5 6 ) , [7,8,9] ] = ( 7 8 9 )
fromLists [ [1,2,3 ] ( 1 2 3 ) , [4,5,6,7] ( 4 5 6 ) , [8,9,0 ] ] = ( 8 9 0 )
toList :: Matrix a -> [a] Source
Get the elements of a matrix stored in a list.
( 1 2 3 ) ( 4 5 6 ) toList ( 7 8 9 ) = [1,2,3,4,5,6,7,8,9]
toLists :: Matrix a -> [[a]] Source
Get the elements of a matrix stored in a list of lists, where each list contains the elements of a single row.
( 1 2 3 ) [ [1,2,3] ( 4 5 6 ) , [4,5,6] toLists ( 7 8 9 ) = , [7,8,9] ]
Accessing
O(1). Get an element of a matrix. Indices range from (1,1) to (n,m).
It returns an error
if the requested element is outside of range.
O(1). Unsafe variant of getElem
, without bounds checking.
safeGet :: Int -> Int -> Matrix a -> Maybe a Source
Variant of getElem
that returns Maybe instead of an error.
safeSet :: a -> (Int, Int) -> Matrix a -> Maybe (Matrix a) Source
Variant of setElem
that returns Maybe instead of an error.
getDiag :: Matrix a -> Vector a Source
O(min rows cols). Diagonal of a not necessarily square matrix.
getMatrixAsVector :: Matrix a -> Vector a Source
Manipulating matrices
:: a | New value. |
-> (Int, Int) | Position to replace. |
-> Matrix a | Original matrix. |
-> Matrix a | Matrix with the given position replaced with the given value. |
Replace the value of a cell in a matrix.
:: a | New value. |
-> (Int, Int) | Position to replace. |
-> Matrix a | Original matrix. |
-> Matrix a | Matrix with the given position replaced with the given value. |
Unsafe variant of setElem
, without bounds checking.
transpose :: Matrix a -> Matrix a Source
O(rows*cols). The transpose of a matrix. Example:
( 1 2 3 ) ( 1 4 7 ) ( 4 5 6 ) ( 2 5 8 ) transpose ( 7 8 9 ) = ( 3 6 9 )
Set the size of a matrix to given parameters. Use a default element for undefined entries if the matrix has been extended.
:: a | Element to add when extending. |
-> Int | Minimal number of rows. |
-> Int | Minimal number of columns. |
-> Matrix a | |
-> Matrix a |
Extend a matrix to a given size adding a default element. If the matrix already has the required size, nothing happens. The matrix is never reduced in size. Example:
( 1 2 3 0 0 ) ( 1 2 3 ) ( 4 5 6 0 0 ) ( 4 5 6 ) ( 7 8 9 0 0 ) extendTo 0 4 5 ( 7 8 9 ) = ( 0 0 0 0 0 )
The definition of extendTo
is based on setSize
:
extendTo e n m a = setSize e (max n $ nrows a) (max m $ ncols a) a
inverse :: (Fractional a, Eq a) => Matrix a -> Either String (Matrix a) Source
O(rows*rows*rows) = O(cols*cols*cols). The inverse of a square matrix. Uses naive Gaussian elimination formula.
rref :: (Fractional a, Eq a) => Matrix a -> Either String (Matrix a) Source
O(rows*rows*cols). Converts a matrix to reduced row echelon form, thus solving a linear system of equations. This requires that (cols > rows) if cols < rows, then there are fewer variables than equations and the problem cannot be solved consistently. If rows = cols, then it is basically a homogenous system of equations, so it will be reduced to identity or an error depending on whether the marix is invertible (this case is allowed for robustness).
:: (Int -> a -> a) | Function takes the current column as additional argument. |
-> Int | Row to map. |
-> Matrix a | |
-> Matrix a |
O(rows*cols). Map a function over a row. Example:
( 1 2 3 ) ( 1 2 3 ) ( 4 5 6 ) ( 5 6 7 ) mapRow (\_ x -> x + 1) 2 ( 7 8 9 ) = ( 7 8 9 )
:: (Int -> a -> a) | Function takes the current row as additional argument. |
-> Int | Column to map. |
-> Matrix a | |
-> Matrix a |
O(rows*cols). Map a function over a column. Example:
( 1 2 3 ) ( 1 3 3 ) ( 4 5 6 ) ( 4 6 6 ) mapCol (\_ x -> x + 1) 2 ( 7 8 9 ) = ( 7 9 9 )
Submatrices
Splitting blocks
O(1). Extract a submatrix given row and column limits. Example:
( 1 2 3 ) ( 4 5 6 ) ( 2 3 ) submatrix 1 2 2 3 ( 7 8 9 ) = ( 5 6 )
:: Int | Row |
-> Int | Column |
-> Matrix a | Original matrix. |
-> Matrix a | Matrix with row |
O(rows*cols). Remove a row and a column from a matrix. Example:
( 1 2 3 ) ( 4 5 6 ) ( 1 3 ) minorMatrix 2 2 ( 7 8 9 ) = ( 7 9 )
:: Int | Row of the splitting element. |
-> Int | Column of the splitting element. |
-> Matrix a | Matrix to split. |
-> (Matrix a, Matrix a, Matrix a, Matrix a) | (TL,TR,BL,BR) |
O(1). Make a block-partition of a matrix using a given element as reference. The element will stay in the bottom-right corner of the top-left corner matrix.
( ) ( | ) ( ) ( ... | ... ) ( x ) ( x | ) splitBlocks i j ( ) = (-------------) , where x = a_{i,j} ( ) ( | ) ( ) ( ... | ... ) ( ) ( | )
Note that some blocks can end up empty. We use the following notation for these blocks:
( TL | TR ) (---------) ( BL | BR )
Where T = Top, B = Bottom, L = Left, R = Right.
Joining blocks
(<|>) :: Matrix a -> Matrix a -> Matrix a Source
Horizontally join two matrices. Visually:
( A ) <|> ( B ) = ( A | B )
Where both matrices A and B have the same number of rows. This condition is not checked.
(<->) :: Matrix a -> Matrix a -> Matrix a Source
Vertically join two matrices. Visually:
( A ) ( A ) <-> ( B ) = ( - ) ( B )
Where both matrices A and B have the same number of columns. This condition is not checked.
joinBlocks :: (Matrix a, Matrix a, Matrix a, Matrix a) -> Matrix a Source
Join blocks of the form detailed in splitBlocks
. Precisely:
joinBlocks (tl,tr,bl,br) = (tl <|> tr) <-> (bl <|> br)
Matrix operations
elementwise :: (a -> b -> c) -> Matrix a -> Matrix b -> Matrix c Source
Perform an operation element-wise.
The second matrix must have at least as many rows
and columns as the first matrix. If it's bigger,
the leftover items will be ignored.
If it's smaller, it will cause a run-time error.
You may want to use elementwiseUnsafe
if you
are definitely sure that a run-time error won't
arise.
elementwiseUnsafe :: (a -> b -> c) -> Matrix a -> Matrix b -> Matrix c Source
Unsafe version of elementwise
, but faster.
Matrix multiplication
About matrix multiplication
Four methods are provided for matrix multiplication.
multStd
: Matrix multiplication following directly the definition. This is the best choice when you know for sure that your matrices are small.multStd2
: Matrix multiplication following directly the definition. However, using a different definition frommultStd
. According to our benchmarks with this version,multStd2
is around 3 times faster thanmultStd
.multStrassen
: Matrix multiplication following the Strassen's algorithm. Complexity grows slower but also some work is added partitioning the matrix. Also, it only works on square matrices of order2^n
, so if this condition is not met, it is zero-padded until this is accomplished. Therefore, its use is not recommended.multStrassenMixed
: This function mixes the previous methods. It provides a better performance in general. Method(
*
)
of theNum
class uses this function because it gives the best average performance. However, if you know for sure that your matrices are small (size less than 500x500), you should usemultStd
ormultStd2
instead, sincemultStrassenMixed
is going to switch to those functions anyway.
We keep researching how to get better performance for matrix multiplication.
If you want to be on the safe side, use (*
).
Functions
multStd :: Num a => Matrix a -> Matrix a -> Matrix a Source
Standard matrix multiplication by definition.
multStd2 :: Num a => Matrix a -> Matrix a -> Matrix a Source
Standard matrix multiplication by definition.
multStrassenMixed :: Num a => Matrix a -> Matrix a -> Matrix a Source
Mixed Strassen's matrix multiplication.
Linear transformations
scaleMatrix :: Num a => a -> Matrix a -> Matrix a Source
Scale a matrix by a given factor. Example:
( 1 2 3 ) ( 2 4 6 ) ( 4 5 6 ) ( 8 10 12 ) scaleMatrix 2 ( 7 8 9 ) = ( 14 16 18 )
scaleRow :: Num a => a -> Int -> Matrix a -> Matrix a Source
Scale a row by a given factor. Example:
( 1 2 3 ) ( 1 2 3 ) ( 4 5 6 ) ( 8 10 12 ) scaleRow 2 2 ( 7 8 9 ) = ( 7 8 9 )
combineRows :: Num a => Int -> a -> Int -> Matrix a -> Matrix a Source
Add to one row a scalar multiple of another row. Example:
( 1 2 3 ) ( 1 2 3 ) ( 4 5 6 ) ( 6 9 12 ) combineRows 2 2 1 ( 7 8 9 ) = ( 7 8 9 )
Switch two rows of a matrix. Example:
( 1 2 3 ) ( 4 5 6 ) ( 4 5 6 ) ( 1 2 3 ) switchRows 1 2 ( 7 8 9 ) = ( 7 8 9 )
Switch two coumns of a matrix. Example:
( 1 2 3 ) ( 2 1 3 ) ( 4 5 6 ) ( 5 4 6 ) switchCols 1 2 ( 7 8 9 ) = ( 8 7 9 )
Decompositions
luDecomp :: (Ord a, Fractional a) => Matrix a -> Maybe (Matrix a, Matrix a, Matrix a, a) Source
Matrix LU decomposition with partial pivoting. The result for a matrix M is given in the format (U,L,P,d) where:
- U is an upper triangular matrix.
- L is an unit lower triangular matrix.
- P is a permutation matrix.
- d is the determinant of P.
- PM = LU.
These properties are only guaranteed when the input matrix is invertible. An additional property matches thanks to the strategy followed for pivoting:
- L_(i,j) <= 1, for all i,j.
This follows from the maximal property of the selected pivots, which also leads to a better numerical stability of the algorithm.
Example:
( 1 2 0 ) ( 2 0 2 ) ( 1 0 0 ) ( 0 0 1 ) ( 0 2 1 ) ( 0 2 -1 ) ( 1/2 1 0 ) ( 1 0 0 ) luDecomp ( 2 0 2 ) = ( ( 0 0 2 ) , ( 0 1 1 ) , ( 0 1 0 ) , 1 )
Nothing
is returned if no LU decomposition exists.
luDecompUnsafe :: (Ord a, Fractional a) => Matrix a -> (Matrix a, Matrix a, Matrix a, a) Source
Unsafe version of luDecomp
. It fails when the input matrix is singular.
luDecomp' :: (Ord a, Fractional a) => Matrix a -> Maybe (Matrix a, Matrix a, Matrix a, Matrix a, a, a) Source
Matrix LU decomposition with complete pivoting. The result for a matrix M is given in the format (U,L,P,Q,d,e) where:
- U is an upper triangular matrix.
- L is an unit lower triangular matrix.
- P,Q are permutation matrices.
- d,e are the determinants of P and Q respectively.
- PMQ = LU.
These properties are only guaranteed when the input matrix is invertible. An additional property matches thanks to the strategy followed for pivoting:
- L_(i,j) <= 1, for all i,j.
This follows from the maximal property of the selected pivots, which also leads to a better numerical stability of the algorithm.
Example:
( 1 0 ) ( 2 1 ) ( 1 0 0 ) ( 0 0 1 ) ( 0 2 ) ( 0 2 ) ( 0 1 0 ) ( 0 1 0 ) ( 1 0 ) luDecomp' ( 2 1 ) = ( ( 0 0 ) , ( 1/2 -1/4 1 ) , ( 1 0 0 ) , ( 0 1 ) , -1 , 1 )
Nothing
is returned if no LU decomposition exists.
luDecompUnsafe' :: (Ord a, Fractional a) => Matrix a -> (Matrix a, Matrix a, Matrix a, Matrix a, a, a) Source
Unsafe version of luDecomp'
. It fails when the input matrix is singular.
cholDecomp :: Floating a => Matrix a -> Matrix a Source
Simple Cholesky decomposition of a symmetric, positive definite matrix. The result for a matrix M is a lower triangular matrix L such that:
- M = LL^T.
Example:
( 2 -1 0 ) ( 1.41 0 0 ) ( -1 2 -1 ) ( -0.70 1.22 0 ) cholDecomp ( 0 -1 2 ) = ( 0.00 -0.81 1.15 )
Properties
trace :: Num a => Matrix a -> a Source
Sum of the elements in the diagonal. See also getDiag
.
Example:
( 1 2 3 ) ( 4 5 6 ) trace ( 7 8 9 ) = 15
diagProd :: Num a => Matrix a -> a Source
Product of the elements in the diagonal. See also getDiag
.
Example:
( 1 2 3 ) ( 4 5 6 ) diagProd ( 7 8 9 ) = 45
Determinants
detLaplace :: Num a => Matrix a -> a Source
Matrix determinant using Laplace expansion.
If the elements of the Matrix
are instance of Ord
and Fractional
consider to use detLU
in order to obtain better performance.
Function detLaplace
is extremely slow.
detLU :: (Ord a, Fractional a) => Matrix a -> a Source
Matrix determinant using LU decomposition. It works even when the input matrix is singular.