Skip to content

GitLab

Commit 9d977504 authored by Christoph Ruegge's avatar Christoph Ruegge
Browse files

Diverses

parent c562b1e5
Hide whitespace changes
Inline Side-by-side
lecture/18 linear algebra-REPA.ipynb deleted 100644 → 0
View file @ c562b1e5
%% Cell type:markdown id: tags:
# Linear algebra
%% Cell type:code id: tags:
``` haskell
{-# LANGUAGE TypeOperators #-}
import Numeric.LinearAlgebra.Repa
import Data.Array.Repa
import Data.Array.Repa.Repr.ForeignPtr
import Data.Array.Repa.Algorithms.Matrix
```
%% Cell type:code id: tags:
``` haskell
-- data for examples
b = fromList (ix1 3) [1,2,3]
c = fromList (ix2 2 3) [1,2..6]
```
%% 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
computeP (c +^ c) :: IO (Array U DIM2 Double)
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
We can use these e.g. to build a scalar product:
%% Cell type:code id: tags:
``` haskell
dotProduct :: Array U DIM1 Double -> Array U DIM1 Double -> Array U DIM0 Double
dotProduct u v = sumS $ u *^ v
dotProduct b b
```
%%%% 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
from `Data.Array.Repa.Algorithms.Matrix`
**mmult** <br />
matrix multiplication
```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
mmultP c c
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
we could thus do matrix-vector multiplications as well: -- TODO ugly
%% Cell type:code id: tags:
``` haskell
b = fromList (ix1 3) [1,2,3] :: Array U DIM1 Double
b' = (computeS $ reshape (ix2 3 1) b) :: Array U DIM2 Double
mmultP c b'
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**trace2** <br />
get the trace
```haskell
trace2S :: Array U DIM2 Double -> Double
trace2P :: Monad m => Array U DIM2 Double -> m Double
```
**transpose2** <br />
transpose the matrix
```haskell
transpose2S :: Array U DIM2 Double -> Array U DIM2 Double
transpose2P :: Monad m => Array U DIM2 Double -> m (Array U DIM2 Double)
```
%% 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)
- functions available as -S and -P versions where applicable, working on delayed representation
%% Cell type:markdown id: tags:
### operators
**cross** <br /> cross product
```haskell
cross :: Array F DIM1 Double -> Array F DIM1 Double -> Array F DIM1 Double
```
%% Cell type:markdown id: tags:
**`dot`** <br /> scalar product
```haskell
dot :: Numeric t => Array F DIM1 t -> Array F DIM1 t -> t
```
%% Cell type:code id: tags:
``` haskell
dot b b
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**norm2** <br /> 2-Norm
we have to reinvent the wheel - for some reasons - however, it's easy:
%% Cell type:code id: tags:
``` haskell
norm2 :: (Floating t, Numeric t) => Array F DIM1 t -> t
norm2 v = sqrt $ dot v v
norm2 b
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**normInf** <br /> infinity norm
%% Cell type:markdown id: tags:
exercise...
%% Cell type:markdown id: tags:
# TODO ab hier überarbeiten
## Matrices
%% Cell type:code id: tags:
``` haskell
(3><4) [1..]
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**fromLists** <br /> creation from nested list
```haskell
fromLists :: Element t => [[t]] -> Matrix t
```
%% Cell type:code id: tags:
``` haskell
fromLists [[-1,-2],[3,4]]
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
### special matrices
**diag** <br /> diagonal matrix
```haskell
diag :: (Num a, Element a) => Vector a -> Matrix a
```
%% Cell type:code id: tags:
``` haskell
diag (4 |> [1,2,3,4])
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**ones** <br/> 1-matrix
```haskell
ones :: Int -> Int -> Matrix Double
```
%% Cell type:code id: tags:
``` haskell
ones 3 3
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**zeros** <br /> 0-matrix
```haskell
zeros :: Int -> Int -> Matrix Double
```
%% Cell type:code id: tags:
``` haskell
zeros 3 3
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**ident** <br/> identity matrix
```haskell
ident :: (Num a, Element a) => Int -> Matrix a
```
%% Cell type:code id: tags:
``` haskell
ident 2
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**(@@>)** <br/> index from tuple (i,j)
%% Cell type:code id: tags:
``` haskell
ident 3 @@> (1,1)
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**find** <br/> find entries from given predicate
```haskell
find :: (e -> Bool) -> c e -> [IndexOf c]
```
%% Cell type:code id: tags:
``` haskell
find (>0) (ident 3 :: Matrix Double)
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
### Sparse matrices (`Numeric.LinearAlgebra.Data`)
Sparse matrices are created via Lists of indices with the type
type AssocMatrix = [((row, column), value)]
which are then used via `mkSparse` to create the sparse matrix.
```haskell
mkSparse :: AssocMatrix -> GMatrix
```
where
```haskell
data GMatrix
```
is the polymorphic matrix for dense, sparse, diagonal, banded, and constant elements.
%% Cell type:code id: tags:
``` haskell
import Numeric.LinearAlgebra.Data
indexlist = [ ((i,j),fromIntegral $ (i+1)*(j+1)) | i<- [0,1..2::Int] , j <- [0,1..2] ]
print indexlist
mkSparse indexlist
```
%%%% Output: display_data
%%%% Output: display_data
%% Cell type:markdown id: tags:
There is also
```haskell
toDense :: AssocMatrix -> Matrix Double
```
and
```haskell
mkDense :: Matrix Double -> GMatrixxxxxx
```
A `GMatrix` only provides basic operations (from `Numeric.LinearAlgebra.HMatrix`):
- It can be applied to a vector using
```haskell
(!#>) :: GMatrix -> Vector Double -> Vector Double
```
%% Cell type:markdown id: tags:
- It can be used in linear systems using
```haskell
cgSolve :: Bool -> GMatrix -> Vector Double -> Vector Double
```
%% Cell type:markdown id: tags:
### matrix operations
**(<>)** <br/> matrix-matrix and matrix-vector multiplication
```haskell
(<>) :: (Product t, Mul a b c) => a t -> b t -> c t
```
**Note**: this is defined in `Numeric.LinearAlgebra.HMatrix` and `Numeric.Container`, with types
%% Cell type:code id: tags:
``` haskell
import Numeric.LinearAlgebra.HMatrix
:t (Numeric.LinearAlgebra.HMatrix.<>)
:t (Numeric.Container.<>)
```
%%%% Output: display_data
%%%% Output: display_data
%% Cell type:markdown id: tags:
the `Numeric.Container` version is more generic, so in most cases, one should
%% Cell type:code id: tags:
``` haskell
import Numeric.LinearAlgebra.HMatrix hiding ((<>))
```
%% Cell type:code id: tags:
``` haskell
c = (2><2) [1::Double .. ]
c
c <> c
```
%%%% Output: display_data
%%%% Output: display_data
%% Cell type:markdown id: tags:
**(*)** <br /> elementwise multiplication
%% Cell type:code id: tags:
``` haskell
c * c
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**trans** <br /> transpose
```haskell
trans :: Matrix t -> Matrix t
```
%% Cell type:code id: tags:
``` haskell
trans c
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
### algorithms
%% Cell type:code id: tags:
``` haskell
import Numeric.LinearAlgebra.Algorithms
```
%% Cell type:markdown id: tags:
**eig** <br /> eigenvalues $\lambda$ and eigenvectors $v$ of a Matrix $A$ in a tuple
$$A v = \lambda v $$
```haskell
eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
```
%% Cell type:code id: tags:
``` haskell
eig c
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**svd** <br /> singular value decomposition of a matrix $M \in K^{m\times n}$ is the decomposition
$$M = U \Sigma V $$
where $U^{m \times m}$, $V^{n \times n}$ are unitary matrix over $K$ , $\Sigma^{m \times n}$ diagonal matrix with non-negative real numbers on the diagonal.
```haskell
svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
```
%% Cell type:code id: tags:
``` haskell
svd c
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**pnorm** <br /> p-norm of a matrix
```haskell
pnorm :: Normed c t => NormType -> c t -> RealOf t
```
%% Cell type:markdown id: tags:
with `NormType`
- `Infinity`
- `PNorm1`
- `PNorm2`
- `Frobeniu`
%% Cell type:code id: tags:
``` haskell
pnorm Infinity c
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
## solve linear systems
**lu** <br /> LU-decomposition of a Matrix $A$
$$A = L U$$
with a lower triangular matrix $L$ and an upper triangular matrix $U$.
```haskell
lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)
```
%% Cell type:code id: tags:
``` haskell
lu $ (2><2) [1,2,3,4 ::Double]
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**rcond** <br /> reciprocal condition of a matrix in respect to a linear equation $Ax = b$ .
it is roughly how big the solution $x$ will change with respect to a change in $b$
$$ cond (A) = \|A^{-1} \| \cdot \| A \|$$
%% Cell type:code id: tags:
``` haskell
rcond c
```
%%%% Output: display_data
%% Cell type:code id: tags:
``` haskell
hilbert n = buildMatrix n n (\(i,j) -> 1 / fromIntegral (i+j +1))
print $ hilbert 3
1/rcond (hilbert 10)
```
%%%% Output: display_data
%%%% Output: display_data
%% Cell type:markdown id: tags:
**`<\>`** <br /> least-squares solution of a linear system of the form $A x = b$ with a SVD
```haskell
(<\>) :: Field t => Matrix t -> b t -> c t
```
%% Cell type:code id: tags:
``` haskell
a = (2><2) [1,2,3,4 ::Double]
b = 2 |> [1,3]
a <\> b
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
**linearSolve** <br /> solution of a linear system with LU decomposition
%% Cell type:code id: tags:
``` haskell
linearSolve a (asColumn b)
```
%%%% Output: display_data
%% Cell type:markdown id: tags: