matrix-0.3.5.0: A native implementation of matrix operations.

Data.Matrix

Description

Matrix datatype and operations.

Every provided example has been tested. Run `cabal test` for further tests.

Synopsis

# Matrix type

data Matrix a Source

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`.

Instances

 Source Source Source Source Eq a => Eq (Matrix a) Source Num a => Num (Matrix a) Source Show a => Show (Matrix a) Source Monoid a => Monoid (Matrix a) Source NFData a => NFData (Matrix a) Source

prettyMatrix :: Show a => Matrix a -> String Source

Display a matrix as a `String` using the `Show` instance of its elements.

nrows :: Matrix a -> Int Source

Number of rows.

ncols :: Matrix a -> Int Source

Number of columns.

O(rows*cols). Similar to `force`. It copies the matrix content dropping any extra memory.

Useful when using `submatrix` from a big matrix.

# Builders

Arguments

 :: Int Rows -> Int Columns -> ((Int, Int) -> a) Generator function -> Matrix a

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 )```

rowVector :: Vector a -> Matrix a Source

O(1). Represent a vector as a one row matrix.

colVector :: Vector a -> Matrix a Source

O(1). Represent a vector as a one column matrix.

## Special matrices

Arguments

 :: Num a => Int Rows -> Int Columns -> Matrix a

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 )```

Arguments

 :: a Default element -> Vector a Diagonal vector -> Matrix a

Similar to `diagonalList`, but using `Vector`, which should be more efficient.

Arguments

 :: 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

Arguments

 :: Int Rows -> Int Columns -> [a] List of elements -> Matrix a

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

Arguments

 :: Int Row -> Int Column -> Matrix a Matrix -> a

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.

(!) :: Matrix a -> (Int, Int) -> a Source

Short alias for `getElem`.

Arguments

 :: Int Row -> Int Column -> Matrix a Matrix -> a

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.

getRow :: Int -> Matrix a -> Vector a Source

O(1). Get a row of a matrix as a vector.

getCol :: Int -> Matrix a -> Vector a Source

O(rows). Get a column of a matrix as a vector.

getDiag :: Matrix a -> Vector a Source

O(min rows cols). Diagonal of a not necessarily square matrix.

O(rows*cols). Transform a `Matrix` to a `Vector` of size rows*cols. This is equivalent to get all the rows of the matrix using `getRow` and then append them, but far more efficient.

# Manipulating matrices

Arguments

 :: 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.

Arguments

 :: 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 )```

Arguments

 :: a Default element. -> Int Number of rows. -> Int Number of columns. -> Matrix a -> Matrix a

Set the size of a matrix to given parameters. Use a default element for undefined entries if the matrix has been extended.

Arguments

 :: 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).

Arguments

 :: (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 )```

Arguments

 :: (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

Arguments

 :: Int Starting row -> Int Ending row -> Int Starting column -> Int Ending column -> Matrix a -> Matrix a

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 )```

Arguments

 :: Int Row `r` to remove. -> Int Column `c` to remove. -> Matrix a Original matrix. -> Matrix a Matrix with row `r` and column `c` removed.

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 )```

Arguments

 :: 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

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 from `multStd`. According to our benchmarks with this version, `multStd2` is around 3 times faster than `multStd`.
• `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 order `2^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 the `Num` 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 use `multStd` or `multStd2` instead, since `multStrassenMixed` 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.

multStrassen :: Num a => Matrix a -> Matrix a -> Matrix a Source

Strassen's matrix multiplication.

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 )```

Arguments

 :: Int Row 1. -> Int Row 2. -> Matrix a Original matrix. -> Matrix a Matrix with rows 1 and 2 switched.

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 )```

Arguments

 :: Int Col 1. -> Int Col 2. -> Matrix a Original matrix. -> Matrix a Matrix with cols 1 and 2 switched.

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.

flatten :: Matrix (Matrix a) -> Matrix a Source

Flatten a matrix of matrices. All sub matrices must have same dimensions This criteria is not checked.