Data.Eigen.LA

The problem: You have a system of equations, that you have written as a single matrix equation

Ax = b

Where A and b are matrices (b could be a vector, as a special case). You want to find a solution x.

The solution: You can choose between various decompositions, depending on what your matrix A looks like, and depending on whether you favor speed or accuracy. However, let's start with an example that works in all cases, and is a good compromise:

import Data.Eigen.Matrix
import Data.Eigen.LA

main = do
    let
        a :: MatrixXd
        a = fromList [[1,2,3], [4,5,6], [7,8,10]]
        b = fromList [[3],[3],[4]]
        x = solve ColPivHouseholderQR a b
    putStrLn "Here is the matrix A:" >> print a
    putStrLn "Here is the vector b:" >> print b
    putStrLn "The solution is:" >> print x

produces the following output

Here is the matrix A:
Matrix 3x3
1.0 2.0 3.0
4.0 5.0 6.0
7.0 8.0 10.0

Here is the vector b:
Matrix 3x1
3.0
3.0
4.0

The solution is:
Matrix 3x1
-2.0000000000000004
1.0000000000000018
0.9999999999999989

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.

Basic linear solving

data Decomposition

Decomposition           Requirements on the matrix          Speed   Accuracy  Rank  Kernel  Image

PartialPivLU            Invertible                          ++      +         -     -       -
FullPivLU               None                                -       +++       +     +       +
HouseholderQR           None                                ++      +         -     -       -
ColPivHouseholderQR     None                                +       ++        +     -       -
FullPivHouseholderQR    None                                -       +++       +     -       -
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)

Constructors

PartialPivLU

LU decomposition of a matrix with partial pivoting.

FullPivLU

LU decomposition of a matrix with complete pivoting.

HouseholderQR

Householder QR decomposition of a matrix.

ColPivHouseholderQR

Householder rank-revealing QR decomposition of a matrix with column-pivoting.

FullPivHouseholderQR

Householder rank-revealing QR decomposition of a matrix with full pivoting.

LLT

Standard Cholesky decomposition (LL^T) of a matrix.

LDLT

Robust Cholesky decomposition of a matrix with pivoting.

JacobiSVD

Two-sided Jacobi SVD decomposition of a rectangular matrix.

Instances

Enum Decomposition

Eq Decomposition

Read Decomposition

Show Decomposition

solve :: Elem a b => Decomposition -> Matrix a b -> Matrix a b -> Matrix a b

x = solve d a b

finds a solution x of ax = b equation using decomposition d

relativeError :: Elem a b => Matrix a b -> Matrix a b -> Matrix a b -> a

e = relativeError x a b

computes norm (ax - b) / norm b where norm is L2 norm

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]] :: MatrixXd
    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 :: Elem a b => Decomposition -> Matrix a b -> Int

The rank of the matrix

kernel :: Elem a b => Decomposition -> Matrix a b -> Matrix a b

Return matrix whose columns form a basis of the null-space of A

image :: Elem a b => Decomposition -> Matrix a b -> Matrix a b

Return a matrix whose columns form a basis of the column-space of A

Multiple linear regression

A linear regression model that contains more than one predictor variable.

linearRegression :: [[Double]] -> ([Double], Double)

(coeffs, error) = linearRegression points

computes multiple linear regression y = a1 x1 + a2 x2 + ... + an xn + b using ColPivHouseholderQR decomposition

• point format is [y, x1..xn]
• coeffs format is [b, a1..an]
• error is calculated using relativeError

import Data.Eigen.LA
main = print $ linearRegression [
    [-4.32, 3.02, 6.89],
    [-3.79, 2.01, 5.39],
    [-4.01, 2.41, 6.01],
    [-3.86, 2.09, 5.55],
    [-4.10, 2.58, 6.32]]

produces the following output

([-2.3466569233817127,-0.2534897541434826,-0.1749653335680988],1.8905965120153139e-3)