Multigrid
This topic brings together many elements we've already talked about
One dimensional boundary value problem
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\)
One dimensional boundary value problem
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\)
One dimensional boundary value problem
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\)
\(\cdots\)
\(\cdots\)
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)\)
First we let:
Obtain an approximate solution by solving \(A\bm{v} = \bm{f}\) with:
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\)
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\)
Let's try this on our model problem
We'll take these as inital guesses
k = 1
k = 6
k = 11
\(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)\)
\(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\)
Vector Components
Error
k = 1
k = 6
k = 11
Moral?
Vector Components
Error
After a few iterations \(\omega=\frac{2}{3}\)
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
Iterations
\(\|e\|_\infty\)
Error reduction stalls, why?
Jacobi
Jacobi SOR
Jacobi
Gauss Siedel
\(x_0, x_1, x_2\)
\(x_j\)
\(x_N\)
\(\cdots\)
\(\cdots\)
\(\cdots\)
\(x_0, x_1, x_2\)
\(x_j\)
\(x_{N/2}\)
\(\cdots\)
\(\cdots\)
\(\Omega_h\)
\(\Omega_{2h}\)
\(\Omega_{2^kh}\)
\(\vdots\)
\(\vdots\)
\(x_0\)
\(x_{N/2^k}\)
\(\cdots\)
Move computation to coarser grids
Why?
\(x_0, x_1, x_2\)
\(x_j\)
\(x_N\)
\(\cdots\)
\(\cdots\)
\(\cdots\)
\(x_0, x_1, x_2\)
\(x_j\)
\(x_{N/2}\)
\(\cdots\)
\(\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
\(x_0, x_1, x_2\)
\(x_j\)
\(x_N\)
\(\cdots\)
\(\cdots\)
\(\cdots\)
\(x_0, x_1, x_2\)
\(x_j\)
\(x_{N/2}\)
\(\cdots\)
\(\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
A general idea
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?
Interpolate \(e\) on the finer grid, and perform residual correction (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\)
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
What about the high frequency error?
Iterating on a coarse grid can correct high frequency errors, the steps are
Similar issues to nested iteration. How do you move to a coarser grid, what's the operator there?
Given an approximation to \(x\), calculate the residual \(r\) then project onto a coarser grid (restriction).
Given an approximation to \(x\), calculate the residual \(r\) then project onto a coarser grid (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}\)
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.
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 is a combination of these idea
Correction scheme
Solve equation in step 4 recursively.
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}\)
Finest Grid
Coarsest Grid
Finest Grid
Storage requirements?
Computation requirements?
Coarsest Grid
Finest Grid
Parallel sparse GMRES solver:
Pulled matrices from the matrix market to test for correctness. Not quite "big data" but large enough to be non-trivial
Low Syncronization GMRES algorithms
A simple strategy for varying the restart parameter in GMRES(m)
A simple strategy for varying the restart parameter in GMRES(m)
Multigrid on non-Cartesian lattices
Ideas?
Mathematical Methods for Engineers II: Multigrid Methods - https://ocw.mit.edu/courses/mathematics/18-086-mathematical-methods-for-engineers-ii-spring-2006/readings/am63.pdf