update1(2021-01-11): C++ Eigen/sparse库

update2(2021-01-15): rust调用C++封装的动态库函数

大家好，之前在论坛里问了不少有关线性代数计算库的问题，现在姑且来交个作业，顺便给出一些用Rust做科学计算的个人经验。结论我就直接放在开头了。

结论

• 因为现阶段Rust生态里没有什么靠谱的稀疏矩阵计算库，所以你的科学计算里包含稀疏矩阵求解形如[A]{x} = {B}或是需要求稀疏矩阵[A]的逆矩阵，又不希望造轮子的话，我完全不推荐使用Rust作为你的编程语言。
• 如果你硬要尝试的话，[A]{x} = {B}推荐使用库sparse21求解。需要注意的是，这个库只支持求解这个公式，这意味着稀疏矩阵[A]本身如果需要通过数个矩阵进行组合的话，需要利用其它库，这里推荐使用库sprs。
• 只需要进行加减乘除计算的话，可以使用库sprs。但是它不支持形如f64 * [稀疏矩阵]的写法。而由于孤儿原则的存在，你没法对其直接进行乘号的重载。直接做法是使用库自带的map函数，非常方便。我个人是使用Enum包装了稀疏矩阵并重载了所有运算符。
• 进行密集矩阵的计算，请直接一步到位，使用库nalgebra。因为库ndarray不支持求逆和求解。但是nalgebra的中文资料非常少，请自行查阅英文官网。不要因为贪图ndarray中文资料多而使用这个库，血泪教训，除非你只需要加减乘除，或是你可以编译通过库ndarray-linalg。（反正我用WINDOWS系统各种失败。简直痛苦。）
• 目前来看，Python的Scipy在求解大型线性方程组（系数为稀疏矩阵时）时仍有碾压性的优势。
• C++的eigen/sparse库或许是比scipy更优的选择。而且比较方便FFI至rust使用。

模型

这次我构建的模型是微分方程的隐式动力学求解，差分格式使用的是Newmark-Beta法配合无条件稳定的参数。与显式动力学不同的是，隐式动力学通常要求解线性方程组[K']{u} = {F'}，其中稀疏矩阵矩阵[K]通常不为主对角矩阵，稀疏矩阵的逆矩阵通常是密集矩阵，导致计算量大增。直接求解{u}可以利用[k]矩阵的稀疏性进行迭代法求解，可以显著降低计算量。

模型原型为Shi et al. 2017描述的关于斜拉索-阻尼器系统的有限差分格式，考虑阻尼器刚度与拉索抗弯刚度的影响。刚度矩阵[K]为五对角矩阵。五对角线上的元素均不为0。主对角线上除了首尾元素均相等，偏移量为1与-1的对角线上除了尾元素均相等，偏移量为2与-2的对角线上的元素均相等。时间增量dt是0.005秒，时间步为1000步，这意味着需要求解1000次[K']{u} = {F'}。且F的值在每个时间步上需要用多个矩阵进行计算并求解。矩阵尺寸由模型分解出的单元数量决定。

Rust开了优化。Python使用scipy库。c++由VS编译，开O2优化。

不同方法的计算耗时 (单位：秒)

分析

1. 第一个方法，由sprs接管所有的矩阵计算，在计算[K']{u} = {F'}时将所有矩阵转化为sparse21的矩阵格式计算完后再转化回sprs的矩阵格式。每个时间步都需要来回转化，来回各一次，所以是2000次。
2. 第二个方法，由于最开始使用了ndarray和sprs导致积重难返。又不打算重构。所以没有纯nalgebra的实现。方法2的Rust(ndarray + nalgebra)意思为，所有计算由ndarray实现，除了在计算逆矩阵时。计算逆矩阵时先转化为nalgebra的DMatrix并求逆，结果再转化回ndarray的矩阵格式。逆矩阵在整个过程中只计算一次。所以只需要来回转化一轮，来回各一次。
3. 第三个方法，计算逆矩阵，然后每次都是用它计算{u} = inv([K']){F'}。
4. 第四个方法，直接求解[K']{u} = {F'}。求解1000次。
5. 第五个方法，c++，迭代法，共轭梯度法。
6. 第六个方法，c++，直接求解法，LU分解。
7. 第七个方法，rust调用使用C接口的C++动态库，C++/eigen/sparse直接LU分解求解（同第六个方法）。只抽象出一个求解器求解AX=B。这意味着，这个求解函数需要被使用上1000次。这1000次中，每次都需要让rust的稀疏矩阵退化为一个Vec，并在动态库中重建成一个Eigen式的稀疏矩阵，求解结果则需要反向转化一轮。这大大增加了消耗时间。时间上与python几乎一致，可以看出python大概也是走的这个流程。如果把循环1000次全托管于C++代码，则速度可以有质的提升，但代码将只能用在求解动力学问题中。

• 显然转化为密集矩阵的方法在矩阵规模提高之后所使用的时间是不可接受的。但对于密集度在0.5%以上的矩阵，无论时间步的数量有多少，都有一定的实用价值。
• sparse21基本是个玩具库。应该没有对矩阵形式进行任何分类求解，而是用的通用方法。但它的计算有个很有意思的地方，在规模在3499和3999时，所用时间相比之前低了非常多。但是导出结果分析后又几乎与Python给结果一致，所以计算本身没什么问题。试了几次，时间误差也就正负一两秒。所以大概是触发了什么奇怪的优化吧？
• 大概是五对角矩阵的逆矩阵仍有一定的稀疏性，或是Python求稀疏矩阵逆的迭代法速度过快，python使用逆矩阵法也有很高的速度优势。
• Python使用scipy的spsolve看来是触发了对五对角矩阵的优化迭代法。计算耗时的增加相比于矩阵规模的增长几乎可以忽略不计。scipy这个库还是十分靠谱的。
• 共轭梯度法在本例中效果不佳。
• 对于稀疏矩阵求解Ax = B。Eigen库采用的是定义求解器solver_sparse类型，再使用solver_sparse.compute(A)输入A参数，后再由solver_sparse.solve(B)求解出答案。显然由于本例的特殊性，LU分解过后仍有相当程度的稀疏性。且在本例中A在循环中是不变的刚度矩阵，而B为力向量，随时间变化，这种特有的求解方式出来的LU只需要分解一次却可以重复使用。大大减少了计算量。

个人感想

• Rust做为静态语言，强制函数输入输出的类型和各类静态检查真的太香了。基本只要不手误写错公式，算出来答案都是对的。Python写的时候心里没底，不报错不一定代表结果是想要的结果。
• Rust编译器确实给力，通过编译器检查后，只出现了一两个小BUG，很容易定位，十几分钟就修好了。
• Python建模大概花了0.5~1秒，而Rust建模时间几乎可以忽略不计。纯Rust的性能还是非常可靠的。Rust离动力学的基础科学计算的距离其实就差了一个稀疏矩阵求解Ax=B。但这个确实又很难。nalgebra的库如果能再给力一点支持稀疏矩阵求解那就真的太香了。
• 目前的生态来看，python还是科研的首选。
• 封装一个cdll之后，效率变得还不错。