数值分析×机器学习使用CNN进行雅可比预条件子的生成（烦）

Posted 2021-07-31 囚生CY

• 英文标题：Machine Learning-Aided Numerical Linear Algebra: Convolutional Neural Networks for the Efficient Preconditioner Generation
• 中文标题：机器学习辅助数值线性代数：用于高效预处理器生成的卷积神经网络
• 论文下载链接：ResearchGate
• 论文项目代码：GitHub@block-prediction

Prescipt：

最近有点伤，不过还是及时调整过来，赶在五一结束把摘要赶出来了，稍微扩充了一些别的内容，觉得邓琪讲得东西还是有点意思的，要开始转变写论文摘要的风格了；

好马真的不吃回头草么？想不到有一天笔者也会困扰于这种问题，我也不知道该不该这样做，严格来说，之前根本不是开始，只是笔者一厢而为罢了，过了这么久了都，害；

谁能来督促我好好学完这学期啊… 就怕到头来啥事都做不成了，那就真的太失败了；

算了，找人去喝酒；

---

序言

刚好这次作业要用到预条件子，这里结合一个实际案例说明预条件子的概念与用途，因为笔者写了一点之后还是没有太弄明白预条件子的意义；

首先我们新建一个稠密度约为$0.5\%$的矩阵：

A = 0.6*speye(1000) + sprand(1000,1000,0.005,1/10000);

如果不使用预条件子，则$\mathrm{G}\mathrm{M}\mathrm{R}\mathrm{E}\mathrm{S}$算法将会需要非常多次的迭代：

b = rand(1000,1); 
[x,~,~,~,resid_plain] = gmres(A,b,50,1e-10,6);  % restart at 50
clf, semilogy(resid_plain,'-')
xlabel('iteration number'), ylabel('residual norm')    % ignore this line
title('Unpreconditioned GMRES')   % ignore this line

在绘制出的图像中可以发现基本上迭代到$300$次残差的二模才降到$10^{-2}$；

因此我们考虑先对矩阵A做个$\mathrm{L}\mathrm{U}$分解：

This version of incomplete LU factorization simply prohibits fill-in for the factors, freezing the sparsity pattern of the approximate factors.

[L,U] = ilu(A);
subplot(121), spy(A)
title('A')   % ignore this line
subplot(122), spy(L)
title('L')   % ignore this line

注意到上面的LU方法是近似的，考察下面的结果就会发现这是非零的，不过差得也不多：

norm(full(A - L*U))

于是事实上的预调件矩阵就是$M=LU$，事实上GMRES方法中是可以直接调用的：

[x,~,~,~,resid_prec] = gmres(A,b,[],1e-10,300,L,U);

最后就可以再做一次图了

clf, semilogy(resid_prec,'-')
xlabel('iteration number'), ylabel('residual norm')    % ignore this line
title('Precondtioned GMRES ')   % ignore this line
axis tight    % ignore this line

总之笔者现在理解的预条件子就是对输入矩阵做预处理，一些常见的机器学习对输入的正则化处理都可以理解为预条件子；

当然上面只是一种基于$\mathrm{L}\mathrm{U}$分解的预条件子，也可以基于$\mathrm{C}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{s}\mathrm{k}\mathrm{y}$分解做一个预条件子；

L = ichol(A+0.05*speye(1000));
[x,~,~,~,residPrec] = minres(A,b,1e-10,400,L,L');
hold on, semilogy(residPrec,'-')
title('Precondtioned MINRES')   % ignore this line
legend('no prec.','with prec.')    % ignore this line

另外关于各种不同的gmres重载方法的参数可以参考博客：GMRES在matlab中的描述

---

文章目录

• 序言
• • @[toc]
• 摘要 Abstract
• 1 引入 Introduction
• 2 背景与相关工作 Background and Related Work
• • 2.1 块雅可比预处理 Block-Jacobi preconditioning
  • 2.2 深度卷积神经网络 Deep convolutional neural networks
• 3 基于神经网络的预条件子模式预测 Preconditioner Pattern Prediction Using Neural Networks
• 4 实验评估 Experimental Evaluation
• • 4.1 数据集 Dataset
  • 4.2 实验环境 Experimental environment
  • 4.3 网络训练 Network training
• 5 总结与展望 Summary and Outlook
• 参考文献
• 图表汇总

摘要 Abstract

• 为高效块雅可比预条件子（effective block-Jacobi preconditioners）稀疏模式（sparsity patterns）是一项具有挑战性且耗用大量计算资源的任务，尤其在一些带有未知源（unknown origin）的问题中会更加困难；
• 本文中设计了一种卷积神经网络（convolutional neural network，下简称为CNN）来检测矩阵稀疏模式（matrix sparsity patterns）中的自然块结构（natural block structures）；
• 对于那些自然块结构（natural block structures）是基于随机分布的非零噪声（nonzero noise）构成的测试矩阵，本文证明了一个训练好的网络能够以超过$95\%$的精确度成功识别强连通分量（strongly connected components），并且由网络得出的块雅可比预条件子（block-Jacobi preconditioners）能够有效地加速迭代的（iterative）$\mathrm{G}\mathrm{M}\mathrm{R}\mathrm{E}\mathrm{S}$求解器；
• 将矩阵分解为若干大小为$128\times 128$的对角块（tile），针对每个对角块，在GPU上生成高效块雅可比预条件子的运行时间低于一毫秒；
• 关键词：
  • Block-Jacobi Preconditioning（块雅可比预条件子）；
  • Convolutional Neural Networks（卷积神经网络）；
  • Multi-label Classification（多标签分类）；

---

1 引入 Introduction

• 在科学计算中，解线性方程系统的迭代过程往往受益于采用适应（adapts to）线性系统性质的复杂预条件子（sophisticated preconditioner）；

• 如果对迭代求解器（iterative solver）进行的收敛性改进可以补偿预条件子的派生（derivation），则预条件子是高效的（efficient）。

  原文：

  A preconditioner is efficient if the convergence improvement rendered to the iterative solver compensates for the derivation of the preconditioner.

• 具体而言，在高效计算的环境下，预条件子的效率取决于两者的并行可延展性（scalability）：

  1. 迭代求解器启动之前的预条件子生成；
  2. 迭代求解器运行的每个步骤中预条件子应用程序（application），我理解就是调用预条件子的算法每一步；

• 基于雅可比（diagonal scaling，即对角化缩放）和块雅可比（block-diagonal scaling，即块对角缩放）可以一定程度提升（moderate improvements）迭代求解器的收敛效果（参考文献[18]），然而这依然是非常具有吸引力的，因为（块）对角化缩放只需要很小的算力开销。

  • 备注：

    1. 定义（雅可比矩阵）：

       设函数$f:{\mathbb{R}}^m\to {\mathbb{R}}^n$，即$\vec{y}=f(\vec{x})$，其中$\vec{y}=(y_1,y_2,...,y_n),\vec{x}=(x_1,x_2,...,x_m)$，则称矩阵$J$为雅可比矩阵：
       J = [ ∂ y 1 ∂ x 1 ∂ y 2 ∂ x 1 . . . ∂ y n ∂ x 1 ∂ y 1 ∂ x 2 ∂ y 2 ∂ x 2 . . . ∂ y n ∂ x 2 . . . . . . . . . . . . ∂ y 1 ∂ x m ∂ y 2 ∂ x m . . . ∂ y n ∂ x m ] m × n J=\\left[\\begin{matrix}\\frac{\\partial y_1}{\\partial x_1}&\\frac{\\partial y_2}{\\partial x_1}&...&\\frac{\\partial y_n}{\\partial x_1}\\\\\\frac{\\partial y_1}{\\partial x_2}&\\frac{\\partial y_2}{\\partial x_2}&...&\\frac{\\partial y_n}{\\partial x_2}\\\\...&...&...&...\\\\\\frac{\\partial y_1}{\\partial x_m}&\\frac{\\partial y_2}{\\partial x_m}&...&\\frac{\\partial y_n}{\\partial x_m}\\end{matrix}\\right]_{m\\times n} J=⎣⎢⎢⎢⎡​∂x1​∂y1​​∂x2​∂y1​​...∂xm​∂y1​​​∂x1​∂y2​​∂x2​∂y2​​...∂xm​∂y2​​​............​∂x1​∂yn​​∂x《数值分析》-- 雅可比迭代法高斯—塞德尔迭代法

       《数值分析》-- 雅可比迭代法高斯—塞德尔迭代法

       《数值分析》-- 雅可比迭代法高斯—塞德尔迭代法

       雅可比矩阵的作用

       雅可比矩阵

       一元函数的梯度和雅可比矩阵是否想用