Multigrid Method - Python实现

  • 算法特征:
    ①. pre-smoothing提取低频残差; ②. 向下插值计算残差补偿; ③.向上插值填充残差补偿; ④. post-smoothing降低整体残差
  • 算法推导:
    Part Ⅰ: 算法原理
    考虑一般线性系统:
    \begin{equation}
    Ax = b
    \label{eq_1}
    \end{equation}
    给定某初始值$x^{0}$, 残差为:
    \begin{equation}
    r^{0} = b - Ax^{0}
    \label{eq_2}
    \end{equation}
    以该残差为基础, 计算残差补偿:
    \begin{equation}
    \delta x^{0} = A^{-1}r^{0}
    \label{eq_3}
    \end{equation}
    填充残差补偿, 式(\ref{eq_1})解表示如下:
    \begin{equation}
    x = x^{0} + \delta x^{0}
    \label{eq_4}
    \end{equation}
    Part Ⅱ: 误差(即残差补偿)特征
    现以Jacobi迭代为例, 分解系数矩阵$A$, 式(\ref{eq_1})转换为
    \begin{equation}
    (L + D + U)x = b
    \label{eq_5}
    \end{equation}
    其中, $L$、$D$、$U$分别为系数矩阵$A$的下三角矩阵、对角矩阵及上三角矩阵. Jacobi迭代公式表示如下:
    \begin{equation}
    x^{k+1} = Jx^k + D^{-1}b
    \label{eq_6}
    \end{equation}
    其中, $J = -D^{-1}(L + U)$, 代表Jacobi迭代矩阵. 定义第$k$步迭代误差:
    \begin{equation}
    e^k = x^k - x
    \label{eq_7}
    \end{equation}
    根据式(\ref{eq_6})可知, 迭代误差满足如下条件:
    \begin{equation}
    e^{k+1} = Je^k
    \label{eq_8}
    \end{equation}
    结合上式, 对矩阵$J$进行特征分解, 则有
    \begin{equation}
    e^k = V\Lambda^kV^{-1}e^0
    \label{eq_9}
    \end{equation}
    其中, 矩阵$V$和$\Lambda$分别由矩阵$J$的本征矢和本征值组成. 显然, Jacobi迭代的收敛性由相应迭代矩阵的本征值决定. 若矩阵$J$所有本征值绝对值均小于1, 则Jacobi迭代收敛; 否则, 不收敛.
    现以有限差分法离散化1维Poisson's equation($\Delta u = f$)为例, 采用Jacobi迭代求解该线性系统, 并分析误差成分:
    \begin{equation}
    \begin{split}
    A &= \begin{bmatrix}-2 & 1 &  &  \\  1 & \ddots & \ddots &  \\ & \ddots & \ddots & 1  \\ &  & 1 & -2 \end{bmatrix} \frac{1}{\Delta x^2}   \\
    D &= -2I\frac{1}{\Delta x^2}  \\
    L + U &= \begin{bmatrix} 0 & 1 &  & \\ 1 & \ddots & \ddots & \\ & \ddots & \ddots & 1 \\ & & 1 & 0 \end{bmatrix} \frac{1}{\Delta x^2}   \\
    J &= \begin{bmatrix} 0 & 1/2 &  & \\ 1/2 & \ddots & \ddots & \\ & \ddots & \ddots & 1/2 \\ & & 1/2 & 0 \end{bmatrix}
    \end{split}
    \label{eq_10}
    \end{equation}
    矩阵$J$之本征矢与本征值分别为:
    \begin{equation}
    \begin{split}
    V_{ij} &= \sin\frac{ij\pi}{n+1}   \\
    \lambda_j &= \cos\frac{j\pi}{n+1}
    \end{split}
    \label{eq_11}
    \end{equation}
    其中, $i$、$j$取$[1, n]$之间的整数. 显然, 矩阵$J$之本征值$\lambda \in (-1, 1)$, Jacobi迭代在该线性系统上收敛.
    观察本征矢取值, 当$j$逐渐增大时, 频率逐渐升高. 由此可将误差大致分为低频组分与高频组分. 低频组分, 误差分布平滑, 利于近似表达; 高频组分, 误差分布不平滑, 不利于近似表达.
    观察本征值取值, 当$j$逐渐增大时, 其变化趋势为: $1 \rightarrow 0 \rightarrow -1$. 低频部分, 误差收敛速度慢; 中频部分, 误差收敛速度快; 高频部分, 误差收敛速度慢.
    对于低频部分, 误差分布平滑, 可采用向下插值技术将离散点上的残差由精细网格过渡至粗糙网格(即降低$n + 1$的取值), 以提升误差组分的频率, 达到快速收敛之目的, 再由向上插值将离散点上的误差由粗糙网格过渡至精细网格, 获取一个可以接受的残差补偿; 对于中频部分, 误差收敛速度快, 无需特殊处理; 对于高频部分, 由于误差收敛速度慢且分布不平滑, 需要引入under-relaxation技术.
    Part Ⅲ: under-relaxation技术
    结合式(\ref{eq_6}), 采用如下凸组合改写Jacobi迭代:
    \begin{equation}
    x^{k+1} = (1 - \alpha)x^k + \alpha [-D^{-1}(L + U)x^k + D^{-1}b], \quad \text{where }\alpha \in (0, 1)
    \label{eq_12}
    \end{equation}
    根据式(\ref{eq_7})之定义, 误差方程转换如下:
    \begin{equation}
    e^{k+1} = [(1 - \alpha )I - \alpha D^{-1}(L + U)]e^k
    \label{eq_13}
    \end{equation}
    此时, 迭代矩阵及其本征值分别为:
    \begin{equation}
    \begin{cases}
    J_\alpha = (1 - \alpha)I - \alpha D^{-1}(L + U)   \\
    \lambda_\alpha = 1 - \alpha + \alpha\lambda
    \end{cases}
    \label{eq_14}
    \end{equation}
    根据式(\ref{eq_11}), 本征值具体取值为:
    \begin{equation}
    \lambda_j^\alpha = 1 - \alpha + \alpha \cdot \cos\frac{j\pi}{n+1} \in (1 - 2\alpha, 1)
    \label{eq_15}
    \end{equation}
    若$\alpha = 0.5$, 则本征值取值区间为$(0, 1)$. 此时已成功消除分布不平滑且收敛速度慢的误差组分, 使Jacobi迭代在多重网格上求解Poisson's equation时快速收敛.
  • 代码实现:
    Part Ⅰ:
    现以Poisson's equation为例进行算法实施, 原始图像、Laplace变换图像及边界条件图像分别如下方左、中、右三图所示:

    Part Ⅱ:
    Multigrid实现如下: 
      1 # 多重网格法求解Poisson's equation
  2 
  3 import copy
  4 import numpy
  5 from PIL import Image
  6 from matplotlib import pyplot as plt
  7 from matplotlib import animation
  8 
  9 
 10 
 11 def get_F(picName):
 12     '''
 13     根据图片数据获取Poisson's equation等式右侧值
 14     '''
 15     img = Image.open(picName)
 16     arr = numpy.array(img).astype(float)
 17     img.close()
 18     
 19     dx = 1 / (arr.shape[1] - 1)
 20     
 21     F = numpy.zeros(arr.shape)
 22     for i in range(1, F.shape[0]-1):
 23         for j in range(1, F.shape[1]-1):
 24             F[i, j] = ((arr[i-1, j] + arr[i+1, j] + arr[i, j-1] + arr[i, j+1]) - 4 * arr[i, j]) / dx ** 2
 25     return F
 26 
 27 
 28     
 29 def get_B(picName):
 30     '''
 31     根据图片数据获取Poisson's equation边界条件, 非边界值填充0
 32     '''
 33     img = Image.open(picName)
 34     B = numpy.array(img).astype(float)
 35     img.close()
 36     B[1:-1, 1:-1] = 0
 37     return B
 38 
 39 
 40 
 41 def jacobi(U, F):
 42     '''
 43     Jacobi迭代法求解Poisson's equation: ΔU = F
 44     U: 迭代解初始值, 含边界条件
 45     F: Poisson's equation右侧值, 边界值无效
 46     横轴长度固定取1
 47     '''
 48     dx = 1 / (U.shape[1] - 1)
 49     UNew = copy.deepcopy(U)
 50     for i in range(1, UNew.shape[0]-1):
 51         for j in range(1, UNew.shape[1]-1):
 52             UNew[i, j] = (U[i-1, j] + U[i+1, j] + U[i, j-1] + U[i, j+1] - dx ** 2 * F[i, j]) / 4
 53     return UNew
 54 
 55 
 56     
 57 class Multigrid(object):
 58     '''
 59     多重网格法之实现
 60     '''
 61     def __init__(self, B, F, alpha=0.5):
 62         self.__B = B                # 边界条件, 矩阵非边界值无效
 63         self.__F = F                # Poisson's equation等式右侧
 64         self.__alpha = alpha        # under-relaxation凸组合因子
 65         
 66         self.__UList = list()       # 记录每帧迭代解
 67         
 68         
 69     def get_solu(self):
 70         return self.__UList
 71         
 72         
 73     def optimize(self, maxIter=3000, epsilon=1.e-1):
 74         '''
 75         maxIter: 最大迭代次数
 76         epsilon: 收敛判据, 迭代解趋于稳定则收敛
 77         '''
 78         self.__UList.clear()
 79         
 80         U0 = copy.deepcopy(self.__B)
 81         F0 = copy.deepcopy(self.__F)
 82         alpha = self.__alpha
 83         
 84         self.__UList.append(numpy.rint(U0).astype("int"))
 85         for idx in range(maxIter):
 86             R0 = self.__calc_residual(U0, F0)
 87             R0Norm = numpy.linalg.norm(R0)
 88             print("iterCnt: {:3d},   ResidualNorm: {}".format(idx, R0Norm))
 89             if R0Norm <= epsilon:
 90                 break
 91             
 92             U0 = self.__do_multigrid(U0, F0, alpha)
 93             self.__UList.append(numpy.rint(U0).astype("int"))
 94             
 95         return U0
 96         
 97     
 98     def __do_multigrid(self, U0, F0, alpha):
 99         # pre-smoothing
100         Uf = self.__smoothing(U0, F0, alpha, 10)
101         Rf = self.__calc_residual(Uf, F0)                     # 残差计算
102         # downward interpolation - find grid to coarse grid
103         Rc = self.__downward_interpolate(Rf)                  # 向下插值 - 残差计算
104         deltaU0 = numpy.zeros(Rc.shape)
105         if (Rc.shape[0] >= 5 and Rc.shape[1] >= 5):
106             deltaUc = self.__do_multigrid(deltaU0, Rc, alpha)
107         else:
108             deltaUc = self.__smoothing(deltaU0, Rf, alpha, 10)
109         # upward interpolation - coarse grid to fine grid
110         deltaUf = self.__upward_interpolate(deltaUc, Uf)      # 向上插值 - 残差补偿计算
111         Uf += deltaUf
112         # post-smoothing
113         U0 = self.__smoothing(Uf, F0, alpha, 10)
114         return U0
115             
116             
117     def __upward_interpolate(self, deltaUc, Uf):
118         deltaUf = numpy.zeros(Uf.shape)
119         deltaUf[2:-2:2, 2:-2:2] = deltaUc[1:-1, 1:-1]
120         deltaUf[1:-1:2, 2:-2:2] = (deltaUf[0:-2:2, 2:-2:2] + deltaUf[2::2, 2:-2:2]) / 2
121         deltaUf[:, 1:-1:2] = (deltaUf[:, 0:-2:2] + deltaUf[:, 2::2]) / 2
122         return deltaUf
123             
124             
125     def __downward_interpolate(self, Rf):
126         '''
127         向下插值 - R中边界值无效
128         '''
129         RcOri = Rf[2:-2:2, 2:-2:2]
130         RcTra = numpy.zeros((RcOri.shape[0]+2, RcOri.shape[1]+2))
131         RcTra[1:-1, 1:-1] = RcOri
132         return RcTra
133     
134     
135     def __calc_residual(self, U, F):
136         '''
137         残差计算
138         '''
139         dx = 1 / (U.shape[1] - 1)
140         R = numpy.zeros(U.shape)
141         
142         for i in range(1, U.shape[0]-1):
143             for j in range(1, U.shape[1]-1):
144                 R[i, j] = F[i, j] - (U[i-1, j] + U[i+1, j] + U[i, j-1] + U[i, j+1] - 4 * U[i, j]) / dx ** 2
145         return R
146             
147             
148     def __smoothing(self, U, F, alpha, iterNum=5):
149         for idx in range(iterNum):
150             U = (1 - alpha) * U + alpha * jacobi(U, F)       # under-relaxation
151         return U
152         
153 
154 
155 class MultigridPlot(object):
156 
157     fig = None
158     ax1 = None
159     line = None
160     txt = None
161     UList = None
162     
163     @classmethod
164     def plot_ani(cls, mgObj):
165         mgObj.optimize()
166         
167         cls.UList = mgObj.get_solu()
168         cls.fig = plt.figure(figsize=(5, 4))
169         cls.ax1 = plt.subplot()
170         cls.line = cls.ax1.imshow(cls.UList[0], cmap="gray")
171         cls.txt = cls.ax1.text(x=100, y=100, s="")
172         
173         ani = animation.FuncAnimation(cls.fig, cls.update, frames=numpy.arange(len(cls.UList)), blit=True, interval=1000, repeat=True)
174         
175         cls.ax1.xaxis.set_visible(False)
176         cls.ax1.yaxis.set_visible(False)
177         cls.fig.tight_layout()
178         ani.save("plot_ani.gif", writer="PillowWriter", fps=3, dpi=1000)
179         plt.close()
180     
181     
182     @classmethod
183     def update(cls, frame):
184         cls.line.set_data(cls.UList[frame])
185         cls.txt.set_text("iterCnt: {}".format(frame))
186         return cls.line, cls.txt
187         
188         
189     @staticmethod
190     def plot_fig(mgObj):
191         U = mgObj.optimize()
192         img = Image.fromarray(U).convert("L")
193         img.save("plot_fig.png")
194         img.close()
195         
196 
197 
198 if __name__ == "__main__":
199     picName = "./telescope-gray.jpg"
200     F = get_F(picName)
201     B = get_B(picName)
202     
203     img = Image.fromarray(F).convert("L")
204     img.save("F.png")
205     img.close()
206     
207     img2 = Image.fromarray(B).convert("L")
208     img2.save("B.png")
209     img2.close()
210     
211     mgObj = Multigrid(B, F)
212     MultigridPlot.plot_fig(mgObj)
213     # MultigridPlot.plot_ani(mgObj)
    View Code

  • 结果展示:
    显然, 通过Multigrid方法, 根据中间Laplace变换图像及右侧边界条件图像即可实现对左侧原始图像的快速求解, 迭代速度较常规迭代方法(Jacobi迭代、Gauss-Seidel迭代等)提升明显.

  • 使用建议:
    ①. 适用于结构化网格, 计算过程在精细网格与粗糙网格间多次变换;
    ②. 适用于求解椭圆型方程, 尤其是线性椭圆型方程.

  • 参考文档: 
    Numerical Methods for PDE 2016 by Professor Qiqi Wang & Jacob White from MIT

posted @ 2020-12-24 21:58  LOGAN_XIONG  阅读(471)  评论(0)  编辑  收藏  举报