哈佛大学｜丢掉神经网络：快速求解物理学逆问题!(附代码算例讲解)

优化离散损失来解决物理学中的逆问题：无需神经网络的快速准确学习

Solving inverse problems in physics by optimizing a discrete loss: Fast and accurate learning without neural networks

在物理学中，逆问题一直是科研人员关注的核心课题之一。逆问题通常涉及通过观测数据来推导出未知的系统参数，这些问题广泛存在于流体动力学、光学、热力学等众多科学与工程领域。然而，逆问题的求解极具挑战，尤其是在需要同时满足物理和几何约束的情况下。随着机器学习的兴起，物理信息神经网络（PINNs）成为了一种热门的方法，但其计算开销巨大且求解速度有限，这也使得人们不断探索更高效的替代方案。

一种全新的逆问题求解框架——ODIL（Optimizing Discrete Loss，优化离散损失）。与传统依赖神经网络进行求解的方法不同，ODIL通过优化离散化后的损失函数，使用传统的偏微分方程（PDE）离散方法，并结合现代机器学习工具，如自动微分，实现了数百倍的速度提升，同时保持甚至超越了PINNs的精度。

PINNs VS ODIL：不同路径的逆问题求解

物理信息神经网络（PINNs）是一种通过神经网络来学习PDE解的方法，自2019年被重新提出后，迅速成为学术界的研究热点。PINNs通过将未知的物理场映射到神经网络的权重上，然后通过最小化物理残差来训练神经网络。这种方法虽然在处理一些复杂的、难以通过传统数值方法求解的逆问题时表现出色，但也面临一些难以忽视的局限：

1. 计算成本高昂：PINNs的核心在于需要不断评估每个训练点上的损失函数，其计算成本随着网络深度和采样点的增加而大幅上升。
2. 收敛性差：由于神经网络的参数众多，PINNs的训练过程往往会遇到收敛性问题，需要大量调参，且优化过程容易陷入局部最优。
3. 难以应用于高阶偏导数方程：在求解高阶导数的物理问题时，PINNs的计算复杂度往往呈指数级增长，导致求解过程极为缓慢。

而ODIL的提出则从根本上改变了这一现状。ODIL通过优化离散化的PDE求解，采用基于网格的离散化方法，直接在离散后的网格上计算损失并进行优化，不再依赖神经网络对解进行全局拟合，这种方法不仅减少了计算资源的消耗，还保持了传统PDE离散方法的高精度、稳定性和收敛性。

ODIL的工作原理及其优势

ODIL的核心在于通过离散化PDE并将其损失函数最小化来求解逆问题。与神经网络不同，ODIL使用传统的数值方法对PDE进行离散，将物理约束直接嵌入到优化问题中。这种方法让ODIL具有以下显著优势：

1. 速度提升五个数量级：ODIL通过使用基于网格的离散化方法并结合梯度优化和牛顿法等经典数值优化方法，相较于PINNs实现了五个数量级的加速。例如，在求解具有19,000个参数的波动方程时，PINNs需要耗时13小时，而ODIL只需0.5秒！
2. 继承经典数值方法的优势：ODIL继承了网格离散化方法的高精度、数值稳定性和收敛性。在涉及复杂流场重建和高阶偏导数问题时，ODIL展现了比PINNs更优越的数值表现。
3. 与机器学习工具无缝集成：ODIL还利用了自动微分等现代机器学习工具，使得求解过程中的梯度计算变得更加方便，这种融合传统数值方法和机器学习技术的方式，为ODIL带来了强大的灵活性和适应性。

实验与应用：五个数量级的提升带来无限可能

ODIL已经在多个物理和工程领域中展现了其强大的应用潜力。例如：

• PDE约束优化：在典型的PDE约束优化问题中，ODIL通过优化离散化的PDE，展现了高效、快速的求解能力，并超越了目前流行的PINNs方法。
• 流动重建：ODIL在重建流动场时，能够以更高的精度、更低的计算成本解决逆问题，特别是在处理Navier-Stokes方程时，ODIL比PINNs表现出更好的收敛性和计算稳定性。
• 光流计算和数据同化：ODIL通过优化离散化损失函数，还能够处理如光流计算、系统识别等复杂任务，适应范围广泛。

在一系列实验中，ODIL的表现都让人惊叹。例如，在解决二维泊松方程时，PINNs需要运行500,000次迭代，而ODIL仅需8次迭代即可达到相同的精度。此外，ODIL在重建流动和复杂数据同化任务中，也展现了其高效且稳健的特性。

ODIL的具体应用案例

ODIL不仅在理论上具有显著的优势，在实际的应用场景中也展现了极强的实践价值。以下是一些ODIL的具体应用案例，帮助大家更好地理解其潜力：

1. 波动方程的求解在解决波动方程时，ODIL展示了其强大的求解能力。与传统PINNs相比，ODIL的求解速度提高了数个数量级，尤其是在涉及复杂边界条件和初始条件的情况下。通过在网格上进行离散化，ODIL避免了大量的全局优化计算，从而显著降低了计算成本。

2. 泊松方程的高效求解在求解泊松方程时，ODIL通过将方程离散化并使用牛顿法进行优化，仅需少量迭代就达到了传统方法的精度。这使得ODIL在工程应用中，尤其是需要快速计算的场景中，具有不可替代的优势。

3. 光流计算中的应用光流计算是一种在计算机视觉中广泛应用的技术，用于估计图像序列中每个像素的运动。ODIL通过将光流问题中的偏微分方程进行离散化，能够快速、精确地计算出运动场，从而在图像处理和视频分析中提供了强有力的工具。

4. 流动重建与系统识别在流体动力学中，流动重建和系统识别是两个重要的问题。ODIL通过对Navier-Stokes方程进行离散化，能够有效地重建复杂流场，并对系统参数进行识别。这对于航空航天、汽车制造等领域的工程应用具有重要意义。

5. 热传导问题中的参数反演在热传导问题中，ODIL被用于反演材料的导热系数等参数。通过对热传导方程进行离散化，并结合实际的温度测量数据，ODIL能够准确地求解出材料的导热特性，为材料科学研究提供了新的工具。

ODIL的局限与未来发展方向

虽然ODIL在许多方面优于PINNs，但其也并非毫无局限。目前，ODIL的应用主要集中在使用网格离散化的方法上，对于那些网格难以定义或计算代价极高的问题，ODIL的表现可能不如PINNs灵活。此外，ODIL的求解速度在很大程度上依赖于离散化的精度和网格的划分质量，这使得它在某些复杂几何形状的求解中需要更多的前处理工作。

未来，ODIL可以向以下几个方向发展：

1. 结合自适应网格：通过引入自适应网格划分技术，ODIL可以在求解过程中根据误差大小自动调整网格密度，从而进一步提高计算效率和精度。
2. 拓展到无网格方法：对于一些难以进行网格划分的问题，未来ODIL可以探索与无网格方法的结合，以提升其在复杂几何形状下的适应性。
3. 与深度学习的融合：虽然ODIL目前不依赖于神经网络，但未来可以探索与深度学习的结合，例如在初始解的预测中引入神经网络，从而减少优化迭代次数，加快收敛速度。

泊松方程案例代码

仅以此代码为例，求解多维立方体上的，带有零Dirichlet边界条件。泊松方程是一种偏微分方程，用于描述在某些领域中的散度平衡问题，适用于诸如电势和温度分布等领域。

1. 泊松方程描述和求解域定义

   """
   Solves the Poisson equation in a multi-dimensional cube
   with zero Dirichlet boundary conditions.
   """
   

   代码求解泊松方程，使用零Dirichlet边界条件，这意味着在立方体边界上，函数值为零。

2. 定义参考解 get_ref_u 函数

   def get_ref_u(name, args, domain):
       xw = domain.points()
       ndim = len(xw)
       if name == 'hat':
           p = 5
           u = np.prod([(1 - x) * x * 5 for x in xw], axis=0)
           u = (u**p / (1 + u**p))**(1 / p)
       elif name == 'osc':
           pi = np.pi
           k = args.osc_k
           x, y = xw
           u = np.sin(pi * (k * x)**2) * np.sin(pi * y)
       else:
           raise ValueError("Unknown name=" + name)
       return u
   

   get_ref_u 用于生成泊松方程的参考解（通常是作为比较的解析解或用于测试的解）。name 参数决定使用哪种类型的解，比如 hat 或 osc（振荡解）。

• 对于 hat 解：该函数定义了一个有界函数，其数学形式类似于 (1-x)*x 的多项式。
• 对于 osc 解：使用了正弦函数来定义振荡解，其特征参数为 osc_k。正弦函数的引入表明解是振荡型的。

def get_ref_rhs(name, args, domain):
    xw = domain.points()
    ndim = len(xw)
    if name == 'osc':
        pi, cos, sin = np.pi, np.cos, np.sin
        k = args.osc_k
        x, y = xw
        fu = (((-4 * k**4 * pi**2 * x**2) - pi**2) * sin(k**2 * pi * x**2) +
              2 * k**2 * pi * cos(k**2 * pi * x**2)) * sin(pi * y)
    else:
        raise ValueError("Unknown name=" + name)
    return fu

def split_wm_wp(st, dirs):
    q = st[0]
    qwm = [st[2 * i + 1] for i in dirs]
    qwp = [st[2 * i + 2] for i in dirs]
    return q, qwm, qwp

def apply_bc_u_mod(st, iw, nw, dirs, mod):
    'Applies zero-Dirichlet boundary conditions.'
    q, qwm, qwp = split_wm_wp(st, dirs)
    zero = mod.cast(0, q.dtype)
    for i in dirs:
        extrap = odil.core.extrap_quadh
        qm = mod.where(iw[i] == 0, extrap(qwp[i], q, zero), qwm[i])
        qp = mod.where(iw[i] == nw[i] - 1, extrap(qwm[i], q, zero), qwp[i])
        qwm[i], qwp[i] = qm, qp
    for i in dirs:
        st[2 * i + 1] = qwm[i]
        st[2 * i + 2] = qwp[i]

def get_discrete_rhs(u, domain, mod):
    ndim = domain.ndim
    dirs = range(ndim)
    dw = domain.step()
    iw = domain.indices()
    nw = domain.size()
    u_st = [None] * (2 * ndim + 1)
    u_st[0] = u
    for i in dirs:
        u_st[2 * i + 1] = mod.roll(u, 1, i)
        u_st[2 * i + 2] = mod.roll(u, -1, i)
    apply_bc_u_mod(u_st, iw, nw, dirs, mod)
    u, uwm, uwp = split_wm_wp(u_st, dirs)
    u_ww = [(uwp[i] - 2 * u + uwm[i]) / dw[i]**2 for i in dirs]
    fu = sum(u_ww)
    return fu

if args.mgloss:
    from functools import partial
    restrict = partial(odil.core.restrict_to_coarser,
                       loc='c' * ndim,
                       mod=mod)
    for _ in range(args.mgloss):
        fu = restrict(fu)
        res += [fu]

def main():
    args = parse_args()
    odil.setup_outdir(args)
    problem, state = make_problem(args)
    callback = odil.make_callback(problem,
                                  args,
                                  plot_func=plot_func,
                                  history_func=history_func,
                                  report_func=report_func)
    odil.util.optimize(args, args.optimizer, problem, state, callback)
    plot_func(problem, state, 0, None, None)

if __name__ == "__main__":
    main()

创新点分析

• 泊松方程的离散化：泊松方程 Δu = f 在代码中通过有限差分方法进行离散化。每个方向的二阶导数通过中心差分 (u_{i+1} - 2*u_i + u_{i-1}) / Δx^2 进行近似。
• 零Dirichlet边界条件：通过 apply_bc_u_mod 函数实现，使得在边界上的解值为零，符合物理上常见的边界条件要求。
• 多维度处理：代码中对 ndim 做了处理，因此可以解决1到多维度的泊松问题。
• 多网格方法：通过将误差限制到更粗的网格来实现更快速的收敛，提高求解效率。
• 完整代码可查看：

https://github.com/cselab/odil/blob/main