Jul 24, 2015
The library includes several solvers for the linear system $AX=B$. The choice depends on the properties of the coefficient matrix $A$.
For well conditioned determined (square) systems we can use linearSolve, based on the LU decomposition.
import Numeric.LinearAlgebra a = (2><2) [ 1, 2 , 3, 4 ] :: Matrix R b = (2><1) [ 5 , 6 :: R ]
λ> linearSolve a b
Just (2><1) [ -3.9999999999999987 , 4.499999999999999 ]
This solver returns a Maybe type because it fails when $A$ is singular.
λ> linearSolve ((2><2) [1,2,3,6]) b
Nothing
We test a bigger system with random elements and multiple right-hand sides:
λ> a <- randn 100 100
λ> dispShort 10 10 2 a
100x100 -0.39 1.54 0.87 -0.11 0.06 0.60 -0.89 .. -0.71 0.38 0.49 0.82 -0.30 0.05 1.01 -0.43 -0.31 .. 1.19 1.64 -2.03 -1.40 1.73 0.45 0.32 -1.38 1.09 .. -0.67 0.40 -0.86 -1.47 1.38 0.08 -0.28 1.12 -0.88 .. -0.54 1.10 0.07 -0.04 0.12 -2.70 0.62 -0.41 -0.30 .. -1.71 1.46 1.46 -1.10 1.02 0.45 -0.82 0.04 0.08 .. -0.00 -0.44 -1.63 -0.18 -1.06 -0.07 -0.89 -0.11 0.18 .. -1.26 3.03 : : : : : : : :: : : 0.01 0.88 0.54 -2.06 1.62 -0.38 -1.66 .. 0.35 1.12 0.80 -0.66 2.54 0.59 -1.18 0.11 -1.11 .. 0.10 -0.33
λ> b <- randn 100 5
λ> b
100x5 0.09 0.27 1.16 -1.85 -0.71 -0.90 1.74 1.99 -0.81 1.79 0.45 -0.53 1.19 0.15 -0.30 -2.37 0.16 1.09 0.45 -1.51 -0.28 0.97 -0.62 -0.95 -1.21 0.21 -1.89 -0.67 0.27 0.38 0.80 -0.14 0.20 1.88 1.17 : : : : : -0.13 -0.03 -0.95 0.24 -2.36 0.66 0.38 0.86 -0.45 0.46
λ> let Just x = linearSolve a b
λ> x
100x5 1.65 1.02 0.28 0.16 0.65 -0.63 0.05 -0.36 -1.84 -0.33 2.13 3.02 6.76 11.46 3.84 9.87 8.74 21.81 34.80 11.61 9.31 9.11 17.88 30.35 10.41 -2.61 -2.23 -4.77 -8.13 -3.11 -2.30 -3.23 -4.88 -7.99 -3.09 : : : : : -5.69 -4.25 -10.52 -16.54 -5.61 1.58 2.11 4.73 9.19 2.86
λ> norm_Frob (a<>x-b)
4.048112530322873e-12
For overdetermined systems we can use linearSolveLS, which computes the least squares solution. We test it with a random system of
200 equations, 100 unknowns, and 10 right-hand sides.
λ> a <- randn 200 100
λ> b <- rand 200 10
λ> let x = linearSolveLS a b
λ> x
100x10 0.04 0.00 0.06 0.05 0.07 0.06 0.01 0.12 0.12 0.01 -0.02 -0.02 -0.08 -0.03 -0.03 -0.05 -0.05 -0.02 -0.06 -0.04 0.05 -0.01 0.01 -0.03 0.00 -0.05 0.04 -0.01 -0.07 -0.02 -0.03 -0.01 -0.05 -0.03 -0.01 0.01 -0.04 0.03 -0.01 -0.06 -0.00 -0.04 -0.02 0.00 -0.03 -0.04 -0.01 -0.00 -0.07 -0.06 -0.12 -0.15 -0.16 -0.10 -0.17 -0.14 -0.17 -0.13 -0.20 -0.15 -0.03 0.00 0.05 0.01 0.04 0.04 -0.01 0.03 0.04 0.01 : : : : : : : : : : -0.03 -0.03 -0.04 -0.00 -0.02 0.06 0.02 0.02 -0.06 -0.04 -0.09 -0.07 -0.09 -0.07 -0.08 -0.10 -0.07 -0.10 -0.11 -0.07
λ> norm_Frob (a<>x - b)
19.055875304607525
If the system is symmetric positive definite we can use cholSolve. As an example, we solve the above overconstrained system using the normal equations: $A^TA X = A^T B$.
λ> let Just ch = mbChol (mTm a)
The function mTm a computes $A^TA$ and stores it in the Herm wrapper requied by the Cholesky decomposition. The result is returned in a Maybe type because this function can fail if the matrix is not positive definite.
λ> let x = cholSolve ch (tr a <> b)
λ> x
100x10 0.04 0.00 0.06 0.05 0.07 0.06 0.01 0.12 0.12 0.01 -0.02 -0.02 -0.08 -0.03 -0.03 -0.05 -0.05 -0.02 -0.06 -0.04 0.05 -0.01 0.01 -0.03 0.00 -0.05 0.04 -0.01 -0.07 -0.02 -0.03 -0.01 -0.05 -0.03 -0.01 0.01 -0.04 0.03 -0.01 -0.06 -0.00 -0.04 -0.02 0.00 -0.03 -0.04 -0.01 -0.00 -0.07 -0.06 -0.12 -0.15 -0.16 -0.10 -0.17 -0.14 -0.17 -0.13 -0.20 -0.15 -0.03 0.00 0.05 0.01 0.04 0.04 -0.01 0.03 0.04 0.01 : : : : : : : : : : -0.03 -0.03 -0.04 -0.00 -0.02 0.06 0.02 0.02 -0.06 -0.04 -0.09 -0.07 -0.09 -0.07 -0.08 -0.10 -0.07 -0.10 -0.11 -0.07
We obtain the same solution and residual:
λ> norm_Frob (a<>x - b)
19.055875304607525
Finally, if the coefficient matrix is rank-deficient, we must use the most general (and slower) solver (<\>), based on the SVD, similar to the backslash operator of Matlab/Octave.
λ> let a = (5><3)[1..] :: Matrix R
λ> a
(5><3) [ 1.0, 2.0, 3.0 , 4.0, 5.0, 6.0 , 7.0, 8.0, 9.0 , 10.0, 11.0, 12.0 , 13.0, 14.0, 15.0 ]
λ> rank a
2
λ> let b = vector [1,3,-2,10,5]
The regular least squares solver (silently) fails
λ> linearSolveLS a (asColumn b)
(3><1) [ -2.1224492054274448e15 , 4.2448984108548815e15 , -2.122449205427437e15 ]
but the more powerful <\> obtains the correct solution:
λ> let x = a <\> b
λ> x
fromList [0.46666666666666734,0.16666666666666718,-0.13333333333333416]
λ> norm_2 (a#>x - b)
7.661592523751182
The operator <\> also admits multiple right-hand sides as columns of a matrix.
Subdetermined systems can also be solved with linearSolveLS (full rank) and <\> (general). In this case we obtain the minimum norm solution.
There are several variants of the singular value decomposition, depending on how many singular vectors must be computed.
λ> let a = (5><3) [1..] :: Matrix R
λ> a
(5><3) [ 1.0, 2.0, 3.0 , 4.0, 5.0, 6.0 , 7.0, 8.0, 9.0 , 10.0, 11.0, 12.0 , 13.0, 14.0, 15.0 ]
This version computes the full set of left and right singular vectors:
λ> let (u,s,v) = svd a
λ> disp 3 u
5x5 -0.101 0.768 -0.016 -0.418 -0.474 -0.249 0.488 0.534 0.277 0.581 -0.396 0.208 -0.815 0.110 0.353 -0.543 -0.072 0.090 0.622 -0.552 -0.690 -0.352 0.206 -0.591 0.092
λ> s
fromList [35.18264833189422,1.4769076999800912,9.904891284970545e-16]
λ> diagRect 0 s 5 3
5x3 35.183 0.000 0.000 0.000 1.477 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
λ> v
3x3 -0.519 -0.751 -0.408 -0.576 -0.046 0.816 -0.632 0.659 -0.408
λ> u <> diagRect 0 s 5 3 <> tr v
5x3 1.000 2.000 3.000 4.000 5.000 6.000 7.000 8.000 9.000 10.000 11.000 12.000 13.000 14.000 15.000
The "thin" svd computes only the min (rows, columns) singular vectors:
λ> let (u,s,v) = thinSVD a
λ> u
5x3 -0.101 0.768 -0.016 -0.249 0.488 0.534 -0.396 0.208 -0.815 -0.543 -0.072 0.090 -0.690 -0.352 0.206
λ> asRow s
1x3 35.183 1.477 0.000
λ> v
3x3 -0.519 -0.751 -0.408 -0.576 -0.046 0.816 -0.632 0.659 -0.408
λ> u <> diag s <> tr v
5x3 1.000 2.000 3.000 4.000 5.000 6.000 7.000 8.000 9.000 10.000 11.000 12.000 13.000 14.000 15.000
The "compact" svd returns the singular vectors corresponding to nonzero singular values:
λ> let (u,s,v) = compactSVD a
λ> u
5x2 -0.101 0.768 -0.249 0.488 -0.396 0.208 -0.543 -0.072 -0.690 -0.352
λ> asRow s
1x2 35.183 1.477
λ> v
3x2 -0.519 -0.751 -0.576 -0.046 -0.632 0.659
λ> u <> diag s <> tr v
5x3 1.000 2.000 3.000 4.000 5.000 6.000 7.000 8.000 9.000 10.000 11.000 12.000 13.000 14.000 15.000
There are also functions to compute only the left or the right singular vectors.
The eigenvalues of a general real matrix can be complex:
λ> let m = matrix 2 [1,-3,4,-1]
λ> m
2x2 1 -3 4 -1
λ> let (l,v) = eig m
λ> l
fromList [0.0 :+ 3.316624790355399,0.0 :+ (-3.316624790355399)]
λ> v
(2><2) [ 0.18898223650461363 :+ 0.6267831705280087, 0.18898223650461363 :+ (-0.6267831705280087) , 0.7559289460184544 :+ 0.0, 0.7559289460184544 :+ (-0.0) ]
If the matrix is symmetric the eigenvalues are real.
λ> let ms = (m + tr m)/2
λ> ms
(2><2) [ 1.0, 0.5 , 0.5, -1.0 ]
λ> eig ms
(fromList [1.118033988749895 :+ 0.0,(-1.118033988749895) :+ 0.0],(2><2) [ 0.9732489894677301 :+ 0.0, (-0.22975292054736118) :+ 0.0 , 0.22975292054736118 :+ 0.0, 0.9732489894677301 :+ 0.0 ])
In this case it is much better to use eigSH, which is more efficient and directly obtains a real result. To avoid symmetry tests, now this function only admits arguments contained in the Herm wrapper, constructed by functions like sym and mTm which always obtain real symmetric or complex Hermitian results.
λ> let ms' = sym m
λ> ms'
Herm (2><2) [ 1.0, 0.5 , 0.5, -1.0 ]
λ> eigSH ms'
(fromList [1.118033988749895,-1.118033988749895],(2><2) [ -0.9732489894677301, 0.22975292054736116 , -0.22975292054736116, -0.9732489894677301 ])
The Cholesky decomposition is used to efficiently solve symmetric positive definite linear systems. It also requires a Herm argument. Since it may fail (if the matrix is not positive definite), it produces a Maybe result.
λ> mbChol ms'
Nothing
λ> let mp = mTm m
λ> mp
Herm (2><2) [ 17.0, -7.0 , -7.0, 10.0 ]
λ> mbChol mp
Just (2><2) [ 4.123105625617661, -1.697749375254331 , 0.0, 2.6678918753996625 ]
λ> mTm <$> mbChol mp
Just (Herm (2><2) [ 17.0, -7.0 , -7.0, 10.0 ])
Let's test the qr decomposition on a random matrix.
λ> m <- randn 5 3
λ> m
5x3 2.992 1.041 0.596 1.036 1.698 -1.023 1.072 0.706 -0.777 0.654 -0.585 -0.604 -1.701 0.656 0.134
λ> let (q,r) = qr m
λ> q
5x5 -0.786 -0.093 0.548 0.175 0.206 -0.272 -0.696 -0.428 0.262 -0.435 -0.282 -0.199 -0.430 -0.646 0.528 -0.172 0.383 -0.528 0.627 0.390 0.447 -0.566 0.230 0.300 0.580
λ> r
5x3 -3.807 -1.085 0.192 0.000 -2.015 0.503 0.000 0.000 1.448 0.000 0.000 0.000 0.000 0.000 0.000
λ> q <> r - m
5x3 0.000 -0.000 -0.000 0.000 0.000 -0.000 0.000 -0.000 -0.000 0.000 -0.000 0.000 -0.000 0.000 -0.000
We can easily find $\mathrm{argmin}_x ||Cx||^2$, $||x||=1$.
c = (5><4) [ 10, 1, 1, 1 , 1, 10, 1, 1 , 1, 1, 10, 1 , 1, 1, 1, 1 , 1, 2, 2, 1] :: Matrix R
λ> nullspace c
(4><0) []
λ> singularValues c
fromList [12.632795935579493,9.03279188134824,9.0,0.9061664738090339]
The solution is the right singular vector corresponding to the smallest singular value:
λ> (last . toColumns . snd . rightSV) c
fromList [-8.934175780982971e-2,-9.575718118885315e-2,-9.575718118885444e-2,0.9867518304077171]
This is already available in the library as null1.
λ> null1 c
fromList [-8.934175780982971e-2,-9.575718118885315e-2,-9.575718118885444e-2,0.9867518304077171]
λ> norm_2 (a #> null1 a)
4.636427468134552e-15
We can use the QR decomposition or the singular vectors of a random matrix to create a "random rotation":
randRot n = snd . rightSV <$> randn n n
λ> disp 3 =<< randRot 5
5x5 -0.034 0.451 -0.325 0.274 -0.784 -0.332 -0.061 0.865 0.216 -0.303 -0.144 0.364 -0.046 0.769 0.503 0.375 -0.722 -0.110 0.535 -0.199 0.853 0.373 0.364 -0.011 0.023
λ> r <- randRot 5
λ> disp 3 r
5x5 0.200 0.880 0.216 -0.369 0.042 0.396 0.022 -0.862 -0.257 -0.183 -0.799 0.332 -0.319 0.129 -0.366 0.310 -0.030 0.283 0.161 -0.893 -0.263 -0.337 0.168 -0.869 -0.184
λ> det r
-1.0000000000000016
This can be used to create a random matrix of given rank and dimensions. We define the function using applicative style:
randRank r (m,n) = (<>) <$> a <*> b
where
a = takeColumns r <$> randRot m
b = takeRows r <$> randRot n
λ> r <- randRank 2 (6,4)
λ> disp 3 r
6x4 0.065 0.007 0.093 0.203 -0.671 0.084 -0.292 -0.574 -0.228 0.080 0.117 0.295 -0.490 0.168 0.237 0.604 -0.152 0.037 0.008 0.040 0.427 -0.099 -0.007 -0.074
λ> rank r
2
λ> rcond r
6.046179037234281e-18
λ> singularValues r
fromList [1.0000000000000002,0.9999999999999997,1.1012104895297248e-16,6.046179037234282e-18]
A similar result can be obtained from the left and right singular vectors of a single random matrix:
randRank' r (n,m) = cutsv r <$> randn n m
where
cutsv r x = u' <> tr v'
where
(u,s,v) = svd x
u' = takeColumns r u
v' = takeColumns r v
λ> r <- randRank' 2 (3,5)
λ> disp 3 r
3x5 -0.296 -0.323 0.292 -0.587 -0.576 0.061 0.101 -0.026 0.139 0.121 -0.077 0.600 0.765 0.183 -0.114
λ> singularValues r
fromList [0.9999999999999998,0.9999999999999998,5.201849001304443e-17]