Data.Matrix

Matrix datatype and operations.

Every provided example has been tested.

Matrix type

data Matrix a

Type of matrices.

Instances

Functor Matrix

Foldable Matrix

Traversable Matrix

Eq a => Eq (Matrix a)

Num a => Num (Matrix a)

Show a => Show (Matrix a)

NFData a => NFData (Matrix a)

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

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

nrows :: Matrix a -> Int

Number of rows.

ncols :: Matrix a -> Int

Number of columns.

forceMatrix :: Matrix a -> Matrix a

O(rows*cols). Similar to force, drop any extra memory.

Useful when using submatrix from a big matrix.

Builders

matrix

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 )

fromList

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

Create a matrix from an non-empty list of non-empty lists. Each list must have the same number of elements. For example:

 fromLists [ [1,2,3]      ( 1 2 3 )
           , [4,5,6]      ( 4 5 6 )
           , [7,8,9] ] =  ( 7 8 9 )

rowVector :: Vector a -> Matrix a

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

colVector :: Vector a -> Matrix a

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

Special matrices

zero

Arguments

:: Num a

=> Int

Rows

-> Int

Columns

-> Matrix a

O(rows*cols). The zero matrix of the given size.

 zero n m =
                 n
   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

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 )

permMatrix

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.

Accessing

getElem

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

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

Short alias for getElem.

safeGet :: Int -> Int -> Matrix a -> Maybe a

Safe variant of getElem.

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

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

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

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

getDiag :: Matrix a -> Vector a

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

Manipulating matrices

setElem

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.

O(1). Replace the value of a cell in a matrix.

transpose :: Matrix a -> Matrix a

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 )

extendTo

Arguments

:: Num a

=> Int

Minimal number of rows.

-> Int

Minimal number of columns.

-> Matrix a

-> Matrix a

Extend a matrix to a given size adding zeroes. 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 4 5 ( 7 8 9 ) = ( 0 0 0 0 0 )

mapRow

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 )

mapCol

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

submatrix

Arguments

:: Int

Starting row

-> Int

Ending row

-> Int

Starting column

-> Int

Ending column

-> Matrix a

-> Matrix a

O(subrows*subcols). 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 )

minorMatrix

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 )

splitBlocks

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)

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.

Implementation is done via slicing of vectors.

Joining blocks

(<|>) :: Matrix a -> Matrix a -> Matrix a

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

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

Join blocks of the form detailed in splitBlocks.

Matrix multiplication

About matrix multiplication

Three 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.
• 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 it is not recommended.
• multStrassenMixed: This function mixes the multStd and multStrassen 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, you should use multStd instead, since multStrassenMixed is going to switch to that function anyway.

Functions

multStd :: Num a => Matrix a -> Matrix a -> Matrix a

Standard matrix multiplication by definition.

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

Strassen's matrix multiplication.

multStrassenMixed :: Num a => Matrix a -> Matrix a -> Matrix a

Mixed Strassen's matrix multiplication.

Linear transformations

scaleMatrix :: Num a => a -> Matrix a -> Matrix a

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

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

Add to one row a scalar multiple of other row. Example:

                   ( 1 2 3 )   (  1  2  3 )
                   ( 4 5 6 )   (  6  9 12 )
 combineRows 2 2 1 ( 7 8 9 ) = (  7  8  9 )

switchRows

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 )

switchCols

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 -> (Matrix a, Matrix a, Matrix a, a)

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 )

luDecomp' :: (Ord a, Fractional a) => Matrix a -> (Matrix a, Matrix a, Matrix a, Matrix a, a, a)

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 is a permutation matrix.
• d,e is the determinant of P,Q.
• 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 )

cholDecomp :: Floating a => Matrix a -> Matrix a

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

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

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

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.

detLU :: (Ord a, Fractional a) => Matrix a -> a

Matrix determinant using LU decomposition.