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

Exercise 20.5

parent abd7e737
%% Cell type:markdown id: tags:
# Exercises: linear algebra
%% Cell type:markdown id: tags:
## problem 1
Implement the following vector norms:
- 1-norm
- p-norm
- infinitiy-norm
For simplicity use Signatures like
```haskell
norm1 :: Vec Double -> Double
```
%% Cell type:code id: tags:
``` haskell
import Data.Array.Repa hiding (map)
import qualified Data.Array.Repa as R
import Numeric.LinearAlgebra.Helpers (norm2, vec, Vec(..))
norm1 :: Vec Double -> Double
norm1 v = foldAllS (\acc -> (acc +) . abs) 0 v
normP :: Double -> Vec Double -> Double
normP p v
| p < 1 = error "normP: no such p-norm"
| otherwise = (**(1/p)) $ foldAllS (\acc -> (acc +) . (**p) . abs) 0 v
normInf :: Vec Double -> Double
normInf = maximum . map abs . toList
-- test
v = vec [1,2,-3]
norm1 v
norm2 v -- from the library
normP 2 v
normInf v
```
%% Cell type:markdown id: tags:
## problem 2
- Generate a $15 \times 15$ Hilbert matrix $H = (h_{ij})_{i,j=1}^{15}$ with $h_{ij} = \frac{1}{i+j-1}$.
- Compute the determinant of $H$.
- Generate the vector $r = (1, 1, \ldots, 1)^T$ and the product $y = H r$
- Solve the equation using
- $LU$ decomposition
- singular value decomposition
- `<\>`
%% Cell type:code id: tags:
``` haskell
import Numeric.LinearAlgebra.Helpers
import Numeric.LinearAlgebra.Repa
import Data.Array.Repa hiding (map)
import qualified Data.Array.Repa as Repa
computeS $ ones (ix1 2) :: Vec Double
```
%% Cell type:code id: tags:
``` haskell
n = 15
hilb :: Mat Double
hilb = computeS $ fromFunction (ix2 n n) $
\(Z:. i :. j) -> 1 / fromIntegral (i+j+1)
rhs :: Vec Double
rhs = hilb `app` computeS (ones (ix1 n))
linearSolve hilb (asColumn rhs)
linearSolveSVD hilb (asColumn rhs)
hilb <\> rhs
```
%% Cell type:markdown id: tags:
## problem 3
Which of the vectors
$$v_1=\left(\begin{array}{c} 1\\2\\3\\4\\5 \end{array}\right),\quad v_2=\left(\begin{array}{c}-1\\27\\26\\1\\-27 \end{array}\right),\quad v_3=\left(\begin{array}{c} 160\\-48\\112\\-160\\48 \end{array}\right),\quad v_4=\left(\begin{array}{c} 120\\234\\-23\\-43\\29 \end{array}\right)$$
are orthogonal?
<!--
v1 = vector([1,2,3,4,5])
v2 = vector([-1,27,26,1,-27])
v3 = vector([160,-48,112,-160,48])
v4 = vector([120,234,-23,-43,29])
-->
%% Cell type:code id: tags:
``` haskell
let vs = (4><5) [ 1, 2, 3, 4, 5
, -1, 27, 26, 1, -27
, 160, -48, 112, -160, 48
, 120, 234, -23, -43, 29
]
in vs `mul` htrans vs
```
%% Cell type:markdown id: tags:
## problem 4
The module `Data.Array.Repa.FFTW` provides functions for Fast Fourier Transforms (FFT) for Repa arrays of up to three dimensions. Use the 1-d version to implement a simple low-pass filter that works as follows:
- First, an FFT is applied to the input signal.
- The coefficients corresponding to frequencies over some pre-determined cutoff are set to 0, the other ones are left untouched.
- Then an inverse FFT is performed.
The map from index `k` of the resulting array to the corresponding frequency for an FFT of length `n` is
```haskell
frq n k
| 2*k < n = k
| otherwise = k - n
```
Use the low-pass filter to denoise a signal created by adding uniformly distributed noise between $-0.03$ and $0.03$ to the function $x \mapsto exp(-x^2)$ for $x$ sampled on 200 points between $-2$ and $2$ and plot the results.
Note:
- The function, its noisy version and the filtered version should be real, but the FFT functions require and return `Complex Double` arrays.
- The FFT functions work with arrays in representation `F` (foreign). To convert from and to `D` or `U`, `copyS` can be used.
- A simple way to generate a list of `n` random numbers without using `IO` is
```haskell
take n . randomRs (-0.03, 0.03) $ mkStdGen 0
```
%% Cell type:code id: tags:
``` haskell
{-# LANGUAGE FlexibleContexts #-}
import Data.Complex
import Data.Array.Repa hiding ((++), map)
import qualified Data.Array.Repa as R
import Data.Array.Repa.Eval
import Data.Array.Repa.FFTW
import Data.Array.Repa.Repr.ForeignPtr
import System.Random
```
%% Cell type:code id: tags:
``` haskell
frq n k
| 2*k < n = k
| otherwise = k - n
```
%% Cell type:code id: tags:
``` haskell
cutoff :: Int -> Array F DIM1 (Complex Double) -> Array F DIM1 (Complex Double)
cutoff cut x = ifft $ computeS $ fromFunction (Z:.n) getElem
where Z:.n = extent x
fx = fft x
getElem (Z:.k)
| abs (frq n k) < cut = index fx (Z:.k)
| otherwise = 0
cutoffReal :: Source r Double => Int -> Array r DIM1 Double -> Array D DIM1 Double
cutoffReal cut = R.map realPart . cutoff cut . computeS . R.map (:+ 0)
```
%% Cell type:code id: tags:
``` haskell
linspace :: Int -> Double -> Double -> Array D DIM1 Double
linspace n a b = fromFunction (Z:.n) $ \(Z:.i) -> a + fromIntegral i / (fromIntegral n - 1) * (b - a)
n :: Int
n = 200
xs :: Array U DIM1 Double
xs = computeS $ linspace n (-2) 2
func :: Array U DIM1 Double
func = computeS $ R.map (\x -> exp (-x^2)) xs
noise :: Array U DIM1 Double
noise = fromList (Z:.n) $ take n $ randomRs (-0.03, 0.03) $ mkStdGen 0
func' :: Array U DIM1 Double
func' = computeS $ func +^ noise
```
%% Cell type:code id: tags:
``` haskell
import Graphics.Rendering.Chart.Easy
import Graphics.Rendering.Chart.Backend.Cairo
toRenderable $ plot (line "exact" [zip (toList xs) (toList func)])
toRenderable $ plot (line "noisy" [zip (toList xs) (toList func')])
toRenderable $ plot (line "smoothed" [zip (toList xs) (toList $ cutoffReal 30 func')])
```
%% Cell type:markdown id: tags:
## problem 5
Implement a parallel Jacobi solver for linear systems using Repa arrays. Since Repa's `mmultP` (from `Data.Array.Repa.Algorithms.Matrix`) only works on `DIM2` arrays, the right hand side should be a matrix as well, i.e. you should effectively solve
$$A X = Y$$
for $X$, where $A$, $X$ and $Y$ are matrices. "Conventional" linear systems can then be solved by taking $X$ and $Y$ as matrices with only a single column.
In order to actually use the parallelisation, this exercise should not be done in the notebook interface. You should also test your program on the 1d Poisson example from lecture 21 and examine its behaviour in ThreadScope. HMatrix matrices can be converted to Repa `DIM2` arrays using `matrixToRepa` from `Data.Packed.Repa`.
In order to actually use the parallelisation, this exercise should not be done in the notebook interface. You should also test your program on the 1d Poisson example from lecture 22 and examine its behaviour in ThreadScope.
%% Cell type:code id: tags:
``` haskell
{-# LANGUAGE BangPatterns #-}
import Data.Array.Repa as R
import Data.Array.Repa.Algorithms.Matrix
import Numeric.LinearAlgebra.Helpers
jacobi :: Monad m => Array U DIM2 Double -> Array U DIM2 Double -> m (Array U DIM2 Double)
-- simple stopping rule: always do 500 iterations
jacobi a y = jacobi' 500 $ computeUnboxedS $ zeros (extent y)
where
d = fromFunction (extent y) $ \(Z:.(!j):._) -> unsafeIndex a (Z:.j:.j)
r = computeUnboxedS $ R.traverse a id $ \get idx@(Z:.(!j):.(!k)) ->
if j==k then 0 else get idx
jacobi' 0 x = return x
jacobi' n x = do
rx <- r `mmultP` x
jacobi' (n-1) =<< computeUnboxedP ((y -^ rx) /^ d)
```
%% Cell type:code id: tags:
``` haskell
a :: Array U DIM2 Double
a = fromListUnboxed (Z:.4:.4) [ 2, -1, 0, 0
, -1, 2, -1, 0
, 0, -1, 2, -1
, 0, 0, -1, 2
]
x :: Array U DIM2 Double
x = fromListUnboxed (Z:.4:.1) [1..4]
y :: Array U DIM2 Double
y = a `mmultS` x
jacobi a y
```
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment