Parallel Numerical Methods

Multigrid

Roadmap

Example problem
Iterative Methods
Residual Correction
Coarse grids
Multigrid

This topic brings together many elements we've already talked about

Example Problem

One dimensional boundary value problem

-u^{\prime\prime}(x) + \sigma u(x) = f(x)

We want to find a function \(u(x)\) satisfying the above.

Model Problem

For \(0 < x < 1\), \(\sigma > 0\) with \(u(0)=u(1)=0\)

Example Problem

One dimensional boundary value problem

-u^{\prime\prime}(x) + \sigma u(x) = f(x)

We want to find a function \(u(x)\) satisfying the above.

Model Problem

For \(0 < x < 1\), \(\sigma > 0\) with \(u(0)=u(1)=0\)

Example Problem

One dimensional boundary value problem

-u^{\prime\prime}(x) + \sigma u(x) = f(x)

Grid:

\(h:=\frac{1}{N}, x_j = jh, j = 0,1, \cdots N\)

\(x_0, x_1, x_2\)

\(x_j\)

\(x_N\)

Let \(f_i = f(x_i)\) and \(v_i \approx u(x_i)\), we want to recover \(v_i\)

\(\cdots\)

Example Problem

Taylor Series Approximation

\(u(x_{i+1}) = u(x_i) + hu^\prime(x_i) + \frac{h^2}{2}hu^{\prime\prime}(x_i) + \frac{h^3}{6}u^{\prime\prime\prime}(x_i) + O(h^4)\)

\(u(x_{i-1}) = u(x_i) - hu^\prime(x_i) + \frac{h^2}{2}hu^{\prime\prime}(x_i) - \frac{h^3}{6}u^{\prime\prime\prime}(x_i) + O(h^4)\)

Add these two together, then solve for \(u^{\prime\prime}\). This gives:

\(u^{\prime\prime}(x_i) = \frac{u(x_{i+1})-2u(x_i) + u(x_{i-1})}{h^2} + O(h^2)\)

Example Problem

First we let:

Obtain an approximate solution by solving \(A\bm{v} = \bm{f}\) with:

Residual Correction

If we let \(e = \tilde{x} - x \), that is, \(e\) is the error between an approximation \(\tilde{x}\) and ground truth, then \(A(\tilde{x} + e) = b\)

Let \(r = b-A\tilde{x}\). It follows that \(Ae = r\), so solve for \(e\) approximately, then add that to \(\tilde{x}\).

I.e Compute \(r\), then find \(e = A^{-1}r\) and compute \(x = \tilde{x} - r\)

Experiment

Let's try this on our model problem

\(v_j=\sin(k \pi j/N) \)

Choose \(0<k<N\)

\(f_j = 0\)

This problem has a trivial solution, but is numerically interesting

\(\sigma = 0.1\)

Error Reduction?

Let's try this on our model problem

We'll take these as inital guesses

k = 1

k = 6

k = 11

Jacobi and Gauss Siedel

\(x^*_i = \frac{b_i - \sum_{i\ne j} a_{i,j}x_j}{a_{ii}}\)

Beginning with an initial guess, both methods compute the next iterate by solving for each component of \(x\) in terms of other components of \(x\)

If D, L, and U are diagonal, strict lower triangular, and strict upper triangular portions of A, then Jacobi method can be written

\(x^* = D^{-1}(b-(L+U)x)\)

\(x^*_i = \frac{b_i - \sum_{j<i} a_{i,j}x^*_j -\sum_{j>i} a_{i,j}x_j}{a_{ii}}\)

\(x^* = (L+D)^{-1}(b-Ux)\)

Jacobi SOR

\(x^*_i = (1-\omega)x_i + \omega\frac{b_i - \sum_{i\ne j} a_{i,j}x_j}{a_{ii}}\)

Beginning with an initial guess, both methods compute the next iterate by solving for each component of \(x\) in terms of other components of \(x\)

Error Reduction?

Vector Components

Error

k = 1

k = 6

k = 11

Moral?

Vector Components

Error

After a few iterations \(\omega=\frac{2}{3}\)

Error Reduction?

Inital Guess (sum of our last inital guesses)

After 300 iterations

Jacobi is good at removing high-frequency error, but low frequency remains

Error

Error Reduction?

Iterations

\(\|e\|_\infty\)

Error reduction stalls, why?

Error Reduction?

Jacobi

Jacobi SOR

Error Reduction?

Jacobi

Gauss Siedel

Coarse Grids

\(x_0, x_1, x_2\)

\(x_j\)

\(x_N\)

\(\cdots\)

\(x_0, x_1, x_2\)

\(x_j\)

\(x_{N/2}\)

\(\cdots\)

\(\Omega_h\)

\(\Omega_{2h}\)

\(\Omega_{2^kh}\)

\(\vdots\)

\(x_0\)

\(x_{N/2^k}\)

\(\cdots\)

Move computation to coarser grids

Why?

Coarse Grids

\(x_0, x_1, x_2\)

\(x_j\)

\(x_N\)

\(\cdots\)

\(x_0, x_1, x_2\)

\(x_j\)

\(x_{N/2}\)

\(\cdots\)

\(\Omega_h\)

\(\Omega_{2h}\)

Move computation to coarser grids

Observation:

From wavelet theory, if we filter then downsample, we split the frequency content between high and low bands.

On the coarse grid, we have less work, and the low-frequency content that was on our original grid is now higher frequency

Coarse Grids

\(x_0, x_1, x_2\)

\(x_j\)

\(x_N\)

\(\cdots\)

\(x_0, x_1, x_2\)

\(x_j\)

\(x_{N/2}\)

\(\cdots\)

\(\Omega_h\)

\(\Omega_{2h}\)

Move computation to coarser grids

Observation:

From wavelet theory, if we filter then downsample, we split the frequency content between high and low bands.

On the coarse grid, we have less work, and the low-frequency content that was on our original grid is now higher frequency

Nested Iteration

A general idea

Solve equation on coarsest grid
Refine grid
If not on finest grid, go to 1

Problems?

How do we solve the problem on coarser grids. I.e. what is our operator? What if, at different levels of refinement, we introduce smooth error?

Prolongation

Interpolate \(e\) on the finer grid, and perform residual correction (prolongation).

Prolongation

Connection to wavelets?

Again, upsampling and filtering (linear b-spline)

This is a linear operator from the coarser grid onto the finer grid, we'll call it \(\Uparrow^{h}_{2h} : \Omega^{2h} \to \Omega^h\)

Prolongation

If we're solving \(Ax=b\) on a coarser grid, what's the linear operator \(A\) on that grid?

If we know that, then we just start with an initial guess vector of 0, use an iterative method for a few iterations, then "prolong" on to a finer grid.

We'll come back to this in a bit

Coarse Grid Correction

What about the high frequency error?

Iterating on a coarse grid can correct high frequency errors, the steps are

Solve the problem on the fine grid
Move to a coarser grid
Solve the problem on the coarser grid
Use this solution to update fine grid

Similar issues to nested iteration. How do you move to a coarser grid, what's the operator there?

Restriction

Given an approximation to \(x\), calculate the residual \(r\) then project onto a coarser grid (restriction).

Restriction

Given an approximation to \(x\), calculate the residual \(r\) then project onto a coarser grid (restriction).

Restriction

Connection to wavelets?

Filtering and downsampling

This is a linear operator from the finer grid onto the coarser grid, we'll call it \(\Downarrow^{h}_{2h} : \Omega^{h} \to \Omega^{2h}\)

Restriction

If we're solving \(Ax=b\) on a finer grid, what's the linear operator \(A\) on that grid?

If we know that, then we just start with an initial guess vector of 0, use an iterative method for a few iterations, then "restrict" to a coarser grid, solve for a few iterations, then project back up.

Linear operator levels

Let \(A_h\) be the operator of interest on a grid \(\Omega^h\)

\(A_{2h} = \Downarrow^{h}_{2h} A_h \Uparrow^{h}_{2h}\)

By noting the problem that must be solved on the coarser grids we can obtain the following relationship

Multigrid

Multigrid is a combination of these idea

Nested iteration fixes smooth error
Coarse grid correction fixes oscillatory error

Correction scheme

Iteratively solve \(A_hx=b\) on grid \(\Omega^h\) for a few iterations, with inital guess \(v\)
Compute \(r = b-A_hx\)
Compute \(r^\prime = \Downarrow^{h}_{2h}r \)
Solve \(A_{2h}e = r^\prime\) for \(e\) on \(\Omega^{2h}\)
Correct guess \(v=v + \Uparrow^{h}_{2h}e\)
Iterate \(A_hx=b\) on grid \(\Omega^h\) with guess \(v\) from previous step

Solve equation in step 4 recursively.

Coarse grid correction

Iterate \(A_h\tilde{x} = b\)

Compute \(r = b -A_h\tilde{x} \)

Restrict \(r^\prime = \Downarrow^h_{2h}r\)

Solve \(A_{2h}e =r^\prime\)

Prolong \(e= \Uparrow^h_{2h}e\)

Update \(\tilde{x} =\tilde{x} + e\)

Iterate \(A_h\tilde{x} = b\)

Storage requirement: \(r, b, \tilde{x}\)

V-cycle

Finest Grid

Coarsest Grid

Finest Grid

Storage requirements?

Computation requirements?

Full Multigrid

Coarsest Grid

Finest Grid

Multigrid

Final Project

Parallel sparse GMRES solver:

Currently implemented in CUDA (cuSPARSE and cuBLAS) and Eigen
- Hessenberg matrix and the orthogonal basis for the Krylov subspace are dense, so use cuBLAS to compute them
However, implementations already exist (in ViennaCL and MAGMA)

https://github.com/sheldonucr/GPU-GMRES

https://www.tandfonline.com/doi/full/10.1080/00207160.2017.1280156

Pulled matrices from the matrix market to test for correctness. Not quite "big data" but large enough to be non-trivial