Multigrid method 下降速度的问题

通过课程提供的两个参考文献,自己实现了一个简陋的multigrid solver,并试图解了一下 1d poisson。虽然residual的下降是等比的,但下降倍率只有不到2。然而按照参考文献中的理论,应该至少每一步下降倍率是3才对。所以想请教一下,现在这个下降速度是合理的吗?

运行结果如下:

start err 0.009005
iter 1, err = 0.004613
iter 2, err = 0.002457
iter 3, err = 0.001332
iter 4, err = 0.000732
iter 5, err = 0.000407
iter 6, err = 0.000228
iter 7, err = 0.000129
iter 8, err = 0.000073
iter 9, err = 0.000041
iter 10, err = 0.000023
iter 11, err = 0.000013
iter 12, err = 0.000007
iter 13, err = 0.000004
iter 14, err = 0.000002
iter 15, err = 0.000001
iter 16, err = 0.000001

代码

import numpy as np
import scipy
from scipy.sparse.linalg import factorized
from scipy.sparse import dia_matrix, csr_matrix

class MultiGridSolver:
  def __init__(self, A, Is, w=2/3):
    """
    Inputs:
        A   sparse matrix
        Is  list of tuples (I_f2c, I_c2f), where I_f2c is restriction matrix from
            fine grid to coarse grid, I_c2f is interploation matrix from coarse
            grid to fine grid
        w   parameter for weighted Jacobian method, default 2/3
    """
    self.A = A
    self.Is = Is
    self.size = A.shape[0]
    self.w = w

    # build multi-level weighted Jacobian Rw = I - w D^{-1} A
    self.As = []
    self.Rws = []
    self.D_invs=[]

    A_tmp = A
    for I_f2c, I_c2f in Is:
      self.As.append(A_tmp)

      # build weighted Jacobian
      D = A_tmp.diagonal()
      m, n = A_tmp.shape
      D_inv = scipy.sparse.dia_matrix((1/D, 0), shape=(m, n))
      Rw = scipy.sparse.eye(m) - w * D_inv * A_tmp
      self.Rws.append(Rw)
      self.D_invs.append(1/D)

      # go down to coarse grid
      A_tmp = I_f2c * A_tmp * I_c2f

    # build solver for coarsest grid
    self.coarsest_solver = factorized(A_tmp)

  def solve(self, f, v_init=None, alpha=1, eps=1e-6):
    if v_init is None or v.size() != self.size:
      v = np.zeros(self.size)
    else:
      v = v_init

    res = f - self.A * v

    err = np.linalg.norm(res)
    print("start err %f" % err)

    cnt = 0
    while err > eps:
      cnt+=1
      v_add = self._step(res, 0, alpha)
      v += v_add
      res = f - self.A * v
      err = np.linalg.norm(res)
      print("iter %d, err = %f" % (cnt, err))

    return v

  def _step(self, res, depth, alpha=1):
    if depth == len(self.Rws):
      v_out = self.coarsest_solver(res)
      return v_out

    # relax
    v_relax = np.zeros_like(res)
    for i in range(alpha):
      v_relax = self.Rws[depth] * v_relax + (1-self.w) * res * self.D_invs[depth]

    # to coarse grid and back
    res_relax = res - self.As[depth] * v_relax
    r_coarse = self.Is[depth][0] * res_relax
    v_coarse = self._step(r_coarse, depth+1)
    v_relax = v_relax + self.Is[depth][1] * v_coarse

    # relax again
    for i in range(alpha):
      v_relax = self.Rws[depth] * v_relax + (1-self.w) * res * self.D_invs[depth]

    return v_relax

def poisson_1d(n, depth):
  size = n
  diag = np.ones(n)
  A = dia_matrix(([-diag, 2*diag, -diag], [-1, 0, 1]), shape=(n,n))

  Is = []
  for i in range(depth):
    if n < 100:  # coarse grid not need to be too small
      break

    l = n // 2
    I_tmp = np.zeros((l, 2*l+1))
    for j in range(l):
      I_tmp[j, 2*j] = 1
      I_tmp[j, 2*j+1] = 2
      I_tmp[j, 2*j+2] = 1
    I_f2c = I_tmp[:, :n] / 4
    I_c2f = I_f2c.transpose() * 2

    I_f2c = scipy.sparse.csr_matrix(I_f2c)
    I_c2f = scipy.sparse.csr_matrix(I_c2f)

    Is.append((I_f2c, I_c2f))

    n = n // 2

  multi_solver = MultiGridSolver(A, Is)

  rng = np.random.RandomState(1124)
  f = rng.rand(size)
  f = f - f.mean()
  v = multi_solver.solve(f * 1/size, alpha=1)

  # plot the solution
  import matplotlib.pyplot as plt
  plt.plot(v)
  plt.plot(f)
  plt.show()

if __name__=="__main__":
  poisson_1d(1024, 10)
2 个赞

感觉我们一般用Multigrid看到等比下降就很开心了,一般不太care到底比例是多少。这里面影响因素很多,比如说

  • 用啥smoother?
  • 每层做多少次smoothing?
  • 边界条件长啥样?

所以实际情况很少能精准达到理论推导的residual reduction速率。

1 个赞

感觉multigrid求解效率要在很高维的情况下才能体现优势,虽然比CG或者Jacobian好,但在小一点的问题上还不如预分解来的快。还有一点问题:

  1. 每层做多少次smoothing是试出来的吗?我尝试了一下感觉一次就够,再多也没有什么改善了
  2. 现在每层间的倍率是2,通常在做模拟的时候也是2倍的scale吗?如果是这样的话,大多数multigrid的迭代都得20+层起了吧
1 个赞

但在小一点的问题上还不如预分解来的快

是的…

嗯,我一般都是试的。2-3次往往足够了。如果你的边界条件比较奇特,有时候需要更多次,甚至在边界上多smooth几次。

一般是2。偶尔也会有人在第一层用4。

20+ 层不至于,因为在2D,3D的情况grid一般也就 1024^3 这个级别,一般可能会做个 6 \sim 7 层到 16^3 以下, 然后上PARDISO。

层数不是越多越好,因为coarsen以后问题问题会产生各种变形(比如说,本来不连在一起的两个cells,coarsen以后连在一起了,这样问题就变味了),求出来的解也不已经就能够给原问题带来很好的修正。

1 个赞

那么高分辨率下MGPCG大概能比MICPCG快多少呢?

如果是1D Possion 可以直接分析fourier mode来确定最优的weighnumber.
我记得jacobi smoother 最优的w是2/3. Guass-sidel也是有一个最佳值。
如果你想仅仅追求收敛速度可以试试W cycle 或者 F cycle。
但是相应每一步的计算时间就多了。

1 个赞

对,我用的是2/3。这个值保证了在jacobi迭代一次后高频的一半频率幅度降为1/3,再加上我在最底层是直接解,所以按理说应该是至少1/3的降幅。可能是在粗细格子间转换时有些错位,导致漏掉一些频率。

我记得jacobi smoother 最优的w是2/3. Guass-sidel也是有一个最佳值。

请问关于w的最优取值有相关资料的吗?能否给一个链接?谢谢

你可以在讲MGPCG那节的slide里找到那个reference。

1 个赞