Commit c83114bf authored by Christoph Ruegge's avatar Christoph Ruegge
Browse files

Add sparse matrix section

parent 9d977504
%% Cell type:markdown id: tags:
# Linear algebra
%% Cell type:code id: tags:
``` haskell
{-# LANGUAGE TypeOperators #-}
import Data.Array.Repa
```
%% Cell type:code id: tags:
``` haskell
-- data for examples
m = fromListUnboxed (ix2 3 3) [1..9]
v = fromListUnboxed (ix1 3) [1,2,3]
```
%% Cell type:markdown id: tags:
## basic elementwise operations
```haskell
*^ :: (Num c, Shape sh, Source r1 c, Source r2 c) => Array r1 sh c -> Array r2 sh c -> Array D sh c
+^ :: (Num c, Shape sh, Source r1 c, Source r2 c) => Array r1 sh c -> Array r2 sh c -> Array D sh c
-^ :: (Num c, Shape sh, Source r1 c, Source r2 c) => Array r1 sh c -> Array r2 sh c -> Array D sh c
/^ :: (Fractional c, Shape sh, Source r1 c, Source r2 c) => Array r1 sh c -> Array r2 sh c -> Array D sh c
```
%% Cell type:code id: tags:
``` haskell
computeS (m +^ m) :: Array U DIM2 Double
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
We can use these e.g. to build an inner product:
%% Cell type:code id: tags:
``` haskell
dotProductS :: Shape sh => Array U sh Double -> Array U sh Double -> Double
dotProductS u v = sumAllS $ u *^ v
```
%% Cell type:code id: tags:
``` haskell
dotProductS v v
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
exercise: how to prevent matrices as arguments? can we also prevent vectors of different size being multiplied?
%% Cell type:markdown id: tags:
## matrix operations in repa
%% Cell type:code id: tags:
``` haskell
import Data.Array.Repa.Algorithms.Matrix
```
%% Cell type:markdown id: tags:
### matrix product
```haskell
mmultS :: Array U DIM2 Double -> Array U DIM2 Double -> Array U DIM2 Double
mmultP :: Monad m => Array U DIM2 Double -> Array U DIM2 Double -> m (Array U DIM2 Double)
```
%% Cell type:code id: tags:
``` haskell
m `mmultS` m
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
we could thus do matrix-vector multiplications as well: -- TODO ugly and possibly inefficient
%% Cell type:code id: tags:
``` haskell
m `mmultS` computeS (reshape (ix2 3 1) v)
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
### matrix trace
```haskell
trace2S :: Array U DIM2 Double -> Double
trace2P :: Monad m => Array U DIM2 Double -> m Double
```
### transpose matrix
```haskell
transpose2S :: Array U DIM2 Double -> Array U DIM2 Double
transpose2P :: Monad m => Array U DIM2 Double -> m (Array U DIM2 Double)
```
more generic transposition, for arbitrary representations and higher order arrays:
```haskell
transpose :: (Shape sh, Source r a) => Array r (sh :. Int :. Int) a -> Array D (sh :. Int :. Int) a
```
%% Cell type:markdown id: tags:
## hmatrix bindings
- not many linear algebra tools in *repa* itself
- use the functionality of *hmatrix* via _Numeric.LinearAlgebra.Repa_
- representations used:
- `D` - delayed representation
- `F` - foreign representation (buffers in c heap)
%% Cell type:code id: tags:
``` haskell
import Numeric.LinearAlgebra.Repa
```
%% Cell type:markdown id: tags:
- convenience helper functions (not part of `repa-linear-algebra`, developed by us)
%% Cell type:code id: tags:
``` haskell
import Numeric.LinearAlgebra.Helpers
```
%% Cell type:markdown id: tags:
- provides type aliases
```haskell
type Vec a = Array F DIM1 a
type Mat a = Array F DIM2 a
```
%% Cell type:markdown id: tags:
### matrices
```haskell
(><) :: Storable a => Int -> Int -> a -> Mat a -- the Cartman operator
mat :: Storable a => [[a]] -> Mat a
```
%% Cell type:code id: tags:
``` haskell
(3><2) [1..]
```
%%%% Output: display_data
%% Cell type:code id: tags:
``` haskell
mat [[1,2,3], [4,5,6]]
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
### vectors
```haskell
(|>) :: Storable a => Int -> [a] -> Vec a
vec :: Storable a => [a] -> Vec a
```
%% Cell type:code id: tags:
``` haskell
4 |> [1..]
```
%%%% Output: display_data
%% Cell type:code id: tags:
``` haskell
vec [2,3,6,9]
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
```haskell
asColumn :: Num a => Vec a -> Mat a
asRow :: Num a => Vec a -> Mat a
```
%% Cell type:code id: tags:
``` haskell
asColumn $ vec [1, 2, 3] :: Mat Double
```
%%%% Output: display_data
%% Cell type:code id: tags:
``` haskell
asRow $ vec [1, 2, 3] :: Mat Double
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
### special matrices
**diag** <br /> diagonal matrix
```haskell
diag :: (Num a, Storable a, Source r a) => Array r DIM1 a -> Array D DIM2 a
```
%% Cell type:code id: tags:
``` haskell
computeS $ diag (4 |> [1,2,3,4]) :: Mat Double
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**ident** <br /> identity matrix
```haskell
ident :: (Num a, Storable a) => Int -> Array D DIM2 a
```
%% Cell type:code id: tags:
``` haskell
computeS $ ident 4 :: Mat Double
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**ones** <br/> 1-array
```haskell
ones :: (Num a, Shape sh) => sh -> Array D sh a
```
%% Cell type:code id: tags:
``` haskell
computeS $ ones (Z:.3:.3) :: Mat Double
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**zeros** <br/> 0-array
```haskell
zeros :: (Num a, Shape sh) => sh -> Array D sh a
```
%% Cell type:code id: tags:
``` haskell
computeS $ zeros (Z:.3:.3) :: Mat Double
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
### storage
*hmatrix* functions work on `F` arrays. Each function in _Numeric.LinearAlgebra.Repa_ is available in 5 variations. E.g. for **dot**:
- `dot` takes `F` arrays as input.
- `dotS` takes `D` arrays as input and computes them sequentially.
- `dotSIO` takes `D` arrays as input and computes them sequentially in the `IO` monad.
- `dotP` takes `D` arrays as input and computes them parallelly in an arbitrary monad (*remember*: parallel computations in *repa* always work in a monad).
- `dotPIO` takes `D` arrays as input and computes them parallelly in the `IO` monad.
The `IO` versions can be faster than the non-`IO` versions.
The `F` type constructor is defined in `Data.Array.Repa.Repr.ForeignPtr`
**Note: the `P` functions compute the *argument arrays* in parallel. The algorithms themselves are done by *hmatrix* and are always single-threaded!**
#### Conversion
- **`delay`** converts to `D` in O(1)
- **`computeS`**, **`computeP`** convert from `D` in O(n)
- **`copyS`**, **`copyP`** convert arbitrary representation in O(n) (`copyS = computeS . delay`)
%% Cell type:markdown id: tags:
### matrix product
```haskell
app :: Mat Double -> Mat Double -> Mat Double
```
%% Cell type:code id: tags:
``` haskell
m = (3><3) [1..9]
v = 3 |> [1..3]
```
%% Cell type:code id: tags:
``` haskell
m `mul` m
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
### matrix-vector product
```haskell
app :: Mat Double -> Vec Double -> Vec Double
```
%% Cell type:code id: tags:
``` haskell
m `app` v
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
### inner product
```haskell
dot :: Numeric t => Vec t -> Vec t -> t
```
%% Cell type:code id: tags:
``` haskell
v `dot` v
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
### product with a scalar
- pure repa version
```haskell
scale :: (Num a, Shape sh, Source r a) => a -> Array r sh a -> Array D sh a
scale c = map (* c)
```
- hmatrix version(s)
```haskell
hscale :: (HShape sh, Numeric a) => a -> Array F sh a -> Array F sh a
```
`HShape` essentially means `DIM1` or `DIM2`.
%% Cell type:markdown id: tags:
### other products
- outer product $(x, y) \mapsto xy^T$
```haskell
outer :: Vec Double -> Vec Double -> Mat Double
```
- cross product in $\mathbb{R}^3$
```haskell
cross :: Vec Double -> Vec Double -> Vec Double
```
- kronecker product between matrices
```haskell
kronecker :: Mat Double -> Mat Double -> Mat Double
```
%% Cell type:markdown id: tags:
### sum and product of elements
```haskell
sumElements :: (Numeric t, HShape sh) => Array F sh t -> t
prodElements :: (Numeric t, HShape sh) => Array F sh t -> t
```
%% Cell type:markdown id: tags:
### transpose matrix
```haskell
htrans :: Numeric t => Mat t -> Mat t
```
%% Cell type:markdown id: tags:
### norms
- Frobenius norm of matrix (analogous to 2-norm)
```haskell
norm_Frob :: (Normed (Vector t), Element t) => Array F DIM2 t -> Double
```
- Nuclear norm of matrix (sum of singular values)
```haskell
norm_nuclear :: (Field t, Numeric t) => Array F DIM2 t -> Double
```
- 2-norm of vector "by hand"
%% Cell type:code id: tags:
``` haskell
norm2 :: (Floating t, Numeric t) => Vec t -> t
norm2 x = sqrt $ dot x x
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
- infinity norm
exercise...
%% Cell type:markdown id: tags:
## solving linear systems
%% Cell type:markdown id: tags:
### inverse and pseudoinverse
```haskell
inv :: Field t => Mat t -> Mat t
pinv :: Field t => Mat t -> Mat t
```
%% Cell type:code id: tags:
``` haskell
inv . computeS . diag $ vec [0,1,2]
```
%%%% Output: display_data
%% Cell type:code id: tags:
``` haskell
pinv . computeS . diag $ vec [0,1,2]
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
### condition number
The condition number
$$\mathrm{cond}(A) = \|A^{-1} \| \cdot \| A \|$$
is an upper bound on how much the relative error if the solution $x$ of $A x = y$ increases over the relative error in the RHS $y$.
```haskell
rcond :: Field t => Mat t -> Double
```
computes the **inverse** of the condition number.
%% Cell type:code id: tags:
``` haskell
hilbert :: Int -> Mat Double
hilbert n = computeS . fromFunction (ix2 n n) $ \(Z:.i:.j) -> 1 / fromIntegral (i+j+1)
```
%% Cell type:code id: tags:
``` haskell
hilbert 3
```
%%%% Output: display_data
%% Cell type:code id: tags:
``` haskell
m = hilbert 10
rcond m
```
%%%% Output: display_data
%% Cell type:code id: tags:
``` haskell
x = 10 |> [1..]
y = m `app` x
-- noisy RHS, first component perturbed by 0.001
y' = computeS $ y +^ (10 |> (0.001 : repeat 0)) :: Vec Double
-- solve noisy equation (see below)
m <\> y'
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
### "one-off" solving
#### the "backslash" operator
Solve $A x = y$ using SVD. $y$ can be a vector or a matrix (i.e. solve multiple equation at once).
```haskell
(<\>), solve :: (Field t, HShape sh) => Mat t -> Array F sh t -> Array F sh t
```
%% Cell type:code id: tags:
``` haskell
m = mat [[1, 1, 0], [0, 1, 1], [1, 1, 1]] -- arbitrary non-singular matrix
x = 3 |> [1..]
y = m `app` x
m <\> y
```
%%%% Output: display_data
%% Cell type:code id: tags:
``` haskell
m = computeS . diag $ vec [0,1,2] -- arbitrary singular matrix
y = m `app` x
m <\> y
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
#### other solvers
For these, RHS is always a matrix.
##### least squares
```haskell
linearSolveLS :: Field t => Mat t -> Mat t -> Mat t
```
##### using LU decomposition
```haskell
linearSolve :: Field t => Mat t -> Mat t -> Maybe (Mat t)
```
returns `Nothing` for singular matrix
%% Cell type:markdown id: tags:
### decompositions
When repeatedly solving equations with the same matrix, it is more efficient to use a precomputed decomposition.
#### example: Cholesky decomposition
For positive definite matrices, the Cholesky decomposition $$A = R^T R,$$ with $R$ upper triangular, can be used.
```haskell
chol :: Field t => Herm t => Mat t
```
`Herm t` is a type for Hermitian matrices. Can be obtained by