【HPC 06】数值求解器

2024-1-15

矩阵直接求解, 迭代求解, lu分解求解, lu分解求det, householder过程求相似矩阵, qr分解求特征值

代码：https://github.com/mrzhuzhe/riven/solver

1. 概述

这次尝最简单的实现了一些常见矩阵求解器

直接法求解
LU分解求结果和行列式(det)
QR分解求特征值

本次大部分代码都来自这里

1. 大部分code 例子都来自这里 https://johnfoster.pge.utexas.edu/numerical-methods-book/LinearAlgebra_IterativeSolvers.html
2. 《Numerical Recipes in C》一本免费书网上很多
3. MIT 18.06 linear algebra
4. MIT numerical analysis 18.33

本次用了Eigen这个库作为baseline来对照结果正确性 https://eigen.tuxfamily.org/dox/GettingStarted.html

2.相关工作

2.1 联立线性方程组

fem

对于有限元和图形学中的弹簧质点系统

通过矩阵表示Vertex和Edge的邻近关系：

例如
上图的1 2 3 三个质点根据每个质点的位置x, y, 可以计算之间距离也就是连杆（弹簧）的长度
根据连杆的长度可以依据胡克定律计算相互之间的弹力
根据一个质点和其他质点的位置就可以得出质点间所有的相互作用力

通过守恒关系联立微分方程：

所有质点x y 方向受力的合力一定等于边界条件的支反力
例如上图在下坠过程中，内部合力为0，外力只有重力，质点之间合力为0
而接触地面后，质点之间合力之和为重力之和

通过动量守恒 速度加速度 角动量质量/密度 等守恒关系

我们可以根据这些关系建立线性方程，对系统状态进行求解

例如
已知位置求受力
再根据受力，弹性，阻尼等来求受力带来的加速度变化
就可以求得每个点接下来的速度和状态了

sph

2.2 线性代数和Eigen

因此实际工程中我们会经常碰到

需要求矩阵线性方程的解，分小规模和大规模的

其中大规模线性方程组求解非常非常的的慢，一般无法用直接法或者直接lu分解等O(n^3)方法

更致命的是，直接法,不仅计算量大，往往无法并行，矩阵消项过程每个都依赖于上一行的写回

而要采用一些基于数值迭代的方法

例如对近似正定矩阵可以用 雅各比迭代 或者共轭梯度

但是迭代法也不是万能的，迭代法一般都要满足一些Precondition

实际开发中，我们经常用到Eigen这个库来表示矩阵，用它类提供的矩阵相乘，点乘和小矩阵求解

eigen

3.实验

3.1 直接法求解

# cd /solver
# code /solver/test_simp.cpp
./build/bin/test_simp

打印了消项的过程

不可解的情况 singular

例如:

a + b = 3

a - b = 5

Rank is 2 equal to unknow Number solution is unique

Rank = number of independent 的限制条件(线性规划simplex)

但是如果我们再加一条independent的限制条件例如 3a - b = 1, 这和前两个式子组成的3a-b = 13矛盾，则无解

无穷解的情况 low rank

例如:

a + b = 3

2*a + 2*b = 6

Rank is lower than number of unknow solution is infinite

这个版本缺少一个矩阵重排，对应消项过程中pivot为0的情况

行与行之间存在依赖关系，无法并行化

3.2 LU分解

# cd /solver
# cmake -B build -S .
# cmake -build build
# code: /solver/test_lu.cpp
# run
./build/bin/test_lu

1. lu_factor的过程打印
2. plu_factor 的测试
3. plu_solver 和 eigen solver的结果对比

lu把原本一个线性方程组的求解，拆分成两个线性方程组的求解，
但是相对消项要对A和b都做消项, lu只需要对A做分解，注意这一步复杂度还是 O(n^3)，但是矩阵cols减少了
而两个现行方程组的求解只需在当前列回代消项复杂度从O(n^3)降低到O(n^2)

分解过程中还是行与行之间存在依赖

3.3 QR分解求特征值

# code: /solver/test_eigen.cpp
# run
./build/bin/test_eigen

迭代法求特征值和特征向量

power_method(mat, rows, cols) # in /solver/test_eigen.cpp

这个方法必须要有一个特别大的特征值的时候才会有用

householder过程求相似矩阵

householder(mat, rows, cols);

q*q in household

household 通过构造一个Q矩阵来A_similar = Q * A * Q，会将矩阵转化为一个三角化的相似矩阵，相似矩阵和原矩阵有相同的特征值

其中I=Q*Q，一般Q = I - 2w*w^T

idea household

理想对称正定状况下，转换完毕后的矩阵如上图

gram_schmidt过程求相似矩阵

gram_schmidt(A, Q, rows, cols);

q.tans*q

gram schmidt 过程跟之前类似构造一个 I = Q.transpose * Q

q就很容易得到R R = Q.transpose() * A，如下图

a=qr

QR迭代求特征值

来将 A_similar = Q.transpose() * A * Q 来将 A 对角化，最后pivot中轴线上的值就是特征值

qr_eigen(A, rows, cols)

Q迭代也对矩阵precondition有要求，矩阵为正定才一定能收敛

4. 结论

目前只是初步实现，还需要后续加深理解

后续要分析一下收敛成立的原因
还有precondition

另外还需分析一下并行化的方案

当然最重要的还是要和实际物理问题的求解联系起来

5. 后续工作

1. 实现 SVD 分解
2. 实现FFT
3. 实现GMRES
4. 在实际的物理问题场景中实现

HPC linear algebra

Starofus

你好呀，我是Sans