Multigrid Method

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

  # 多重网格法求解Poisson's equation
  
  import copy
  import numpy
  from PIL import Image
  from matplotlib import pyplot as plt
  from matplotlib import animation
  
  
  
  def get_F(picName):
      '''
      根据图片数据获取Poisson's equation等式右侧值
      '''
      img = Image.open(picName)
      arr = numpy.array(img).astype(float)
      img.close()
      
      dx = 1 / (arr.shape[1] - 1)
      
      F = numpy.zeros(arr.shape)
      for i in range(1, F.shape[0]-1):
          for j in range(1, F.shape[1]-1):
              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
      return F
  
  
      
  def get_B(picName):
      '''
      根据图片数据获取Poisson's equation边界条件, 非边界值填充0
      '''
      img = Image.open(picName)
      B = numpy.array(img).astype(float)
      img.close()
      B[1:-1, 1:-1] = 0
      return B
  
  
  
  def jacobi(U, F):
      '''
      Jacobi迭代法求解Poisson's equation: ΔU = F
      U: 迭代解初始值, 含边界条件
      F: Poisson's equation右侧值, 边界值无效
      横轴长度固定取1
      '''
      dx = 1 / (U.shape[1] - 1)
      UNew = copy.deepcopy(U)
      for i in range(1, UNew.shape[0]-1):
          for j in range(1, UNew.shape[1]-1):
              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
      return UNew
  
  
      
  class Multigrid(object):
      '''
      多重网格法之实现
      '''
      def __init__(self, B, F, alpha=0.5):
          self.__B = B                # 边界条件, 矩阵非边界值无效
          self.__F = F                # Poisson's equation等式右侧
          self.__alpha = alpha        # under-relaxation凸组合因子
          
          self.__UList = list()       # 记录每帧迭代解
          
          
      def get_solu(self):
          return self.__UList
          
          
      def optimize(self, maxIter=3000, epsilon=1.e-1):
          '''
          maxIter: 最大迭代次数
          epsilon: 收敛判据, 迭代解趋于稳定则收敛
          '''
          self.__UList.clear()
          
          U0 = copy.deepcopy(self.__B)
          F0 = copy.deepcopy(self.__F)
          alpha = self.__alpha
          
          self.__UList.append(numpy.rint(U0).astype("int"))
          for idx in range(maxIter):
              R0 = self.__calc_residual(U0, F0)
              R0Norm = numpy.linalg.norm(R0)
              print("iterCnt: {:3d},   ResidualNorm: {}".format(idx, R0Norm))
              if R0Norm <= epsilon:
                  break
              
              U0 = self.__do_multigrid(U0, F0, alpha)
              self.__UList.append(numpy.rint(U0).astype("int"))
              
          return U0
          
      
      def __do_multigrid(self, U0, F0, alpha):
          # pre-smoothing
          Uf = self.__smoothing(U0, F0, alpha, 10)
          Rf = self.__calc_residual(Uf, F0)                     # 残差计算
          # downward interpolation - find grid to coarse grid
          Rc = self.__downward_interpolate(Rf)                  # 向下插值 - 残差计算
          deltaU0 = numpy.zeros(Rc.shape)
          if (Rc.shape[0] >= 5 and Rc.shape[1] >= 5):
              deltaUc = self.__do_multigrid(deltaU0, Rc, alpha)
          else:
              deltaUc = self.__smoothing(deltaU0, Rf, alpha, 10)
          # upward interpolation - coarse grid to fine grid
          deltaUf = self.__upward_interpolate(deltaUc, Uf)      # 向上插值 - 残差补偿计算
          Uf += deltaUf
          # post-smoothing
          U0 = self.__smoothing(Uf, F0, alpha, 10)
          return U0
              
              
      def __upward_interpolate(self, deltaUc, Uf):
          deltaUf = numpy.zeros(Uf.shape)
          deltaUf[2:-2:2, 2:-2:2] = deltaUc[1:-1, 1:-1]
          deltaUf[1:-1:2, 2:-2:2] = (deltaUf[0:-2:2, 2:-2:2] + deltaUf[2::2, 2:-2:2]) / 2
          deltaUf[:, 1:-1:2] = (deltaUf[:, 0:-2:2] + deltaUf[:, 2::2]) / 2
          return deltaUf
              
              
      def __downward_interpolate(self, Rf):
          '''
          向下插值 - R中边界值无效
          '''
          RcOri = Rf[2:-2:2, 2:-2:2]
          RcTra = numpy.zeros((RcOri.shape[0]+2, RcOri.shape[1]+2))
          RcTra[1:-1, 1:-1] = RcOri
          return RcTra
      
      
      def __calc_residual(self, U, F):
          '''
          残差计算
          '''
          dx = 1 / (U.shape[1] - 1)
          R = numpy.zeros(U.shape)
          
          for i in range(1, U.shape[0]-1):
              for j in range(1, U.shape[1]-1):
                  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
          return R
              
              
      def __smoothing(self, U, F, alpha, iterNum=5):
          for idx in range(iterNum):
              U = (1 - alpha) * U + alpha * jacobi(U, F)       # under-relaxation
          return U
          
  
  
  class MultigridPlot(object):
  
      fig = None
      ax1 = None
      line = None
      txt = None
      UList = None
      
      @classmethod
      def plot_ani(cls, mgObj):
          mgObj.optimize()
          
          cls.UList = mgObj.get_solu()
          cls.fig = plt.figure(figsize=(5, 4))
          cls.ax1 = plt.subplot()
          cls.line = cls.ax1.imshow(cls.UList[0], cmap="gray")
          cls.txt = cls.ax1.text(x=100, y=100, s="")
          
          ani = animation.FuncAnimation(cls.fig, cls.update, frames=numpy.arange(len(cls.UList)), blit=True, interval=1000, repeat=True)
          
          cls.ax1.xaxis.set_visible(False)
          cls.ax1.yaxis.set_visible(False)
          cls.fig.tight_layout()
          ani.save("plot_ani.gif", writer="PillowWriter", fps=3, dpi=1000)
          plt.close()
      
      
      @classmethod
      def update(cls, frame):
          cls.line.set_data(cls.UList[frame])
          cls.txt.set_text("iterCnt: {}".format(frame))
          return cls.line, cls.txt
          
          
      @staticmethod
      def plot_fig(mgObj):
          U = mgObj.optimize()
          img = Image.fromarray(U).convert("L")
          img.save("plot_fig.png")
          img.close()
          
  
  
  if __name__ == "__main__":
      picName = "./telescope-gray.jpg"
      F = get_F(picName)
      B = get_B(picName)
      
      img = Image.fromarray(F).convert("L")
      img.save("F.png")
      img.close()
      
      img2 = Image.fromarray(B).convert("L")
      img2.save("B.png")
      img2.close()
      
      mgObj = Multigrid(B, F)
      MultigridPlot.plot_fig(mgObj)
      # MultigridPlot.plot_ani(mgObj)

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

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

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