大家好,之前在论坛里问了不少有关线性代数计算库的问题,现在姑且来交个作业,顺便给出一些用做科学计算的个人经验。结论我就直接放在开头了。
结论
因为现阶段Rust生态里没有什么靠谱的稀疏矩阵计算库,所以你的科学计算里包含稀疏矩阵求解形如或是需要求稀疏矩阵的逆矩阵,又不希望造轮子的话,我完全不推荐使用作为你的编程语言。
如果你硬要尝试的话,推荐使用库求解。需要注意的是,这个库只支持求解这个公式,这意味着稀疏矩阵本身如果需要通过数个矩阵进行组合的话,需要利用其它库,这里推荐使用库。
只需要进行加减乘除计算的话,可以使用库。但是它不支持形如的写法。而由于孤儿原则的存在,你没法对其直接进行乘号的重载。直接做法是使用库自带的函数,非常方便。我个人是使用包装了稀疏矩阵并重载了所有运算符。
进行密集矩阵的计算,请直接一步到位,使用库。因为库不支持求逆和求解。但是的中文资料非常少,请自行查阅英文官网。不要因为贪图中文资料多而使用这个库,血泪教训,除非你只需要加减乘除,或是你可以编译通过库。(反正我用WINDOWS系统各种失败。简直痛苦。)
目前来看,的在求解大型线性方程组(系数为稀疏矩阵时)时仍有碾压性的优势。
模型
这次我构建的摸准模型是微分方程的隐式动力学求解,差分格式使用的是Newmark-Beta法配合无条件稳定的参数。与显式动力学不同的是,隐式动力学通常要求解线性方程组,其中稀疏矩阵矩阵通常不为主对角矩阵,稀疏矩阵的逆矩阵通常是密集矩阵,导致计算量大增。直接求解可以利用矩阵的稀疏性进行迭代法求解,可以显著降低计算量。
模型原型为Shi et al. 2017描述的关于斜拉索-阻尼器系统的有限差分格式,考虑阻尼器刚度与拉索抗弯刚度的影响。刚度矩阵为五对角矩阵。五对角线上的元素均不为0。主对角线上除了首尾元素均相等,偏移量为1与-1的对角线上除了尾元素均相等,偏移量为2与-2的对角线上的元素均相等。时间增量是,时间步为,这意味着需要求解1000次。且的值在每个时间步上需要用多个矩阵进行计算并求解。矩阵尺寸由模型分解出的单元数量决定。
Rust开了优化。Python使用scipy库。
不同方法的计算耗时 (单位:秒)
分析
第一个方法由接管所有的矩阵计算,在计算时将所有矩阵转化为的矩阵格式计算完后再转化回的矩阵格式。每个时间步都需要来回转化,来回各一次,所以是2000次。
由于最开始使用了和导致积重难返。又不打算重构。所以没有纯的实现。方法2的意思为,所有计算由实现,除了在计算逆矩阵时。计算逆矩阵时先转化为的并求逆,结果再转化回的矩阵格式。逆矩阵在整个过程中只计算一次。所以只需要来回转化一轮,来回各一次。
计算逆矩阵,然后每次都是用它计算。
直接求解。求解1000次。
显然转化为密集矩阵的方法在矩阵规模提高之后所使用的时间是不可接受的。但对于密集度在0.5%以上的矩阵,无论时间步的数量有多少,都有一定的实用价值。
基本是个玩具库。应该没有对矩阵形式进行任何分类求解,而是用的通用方法。但它的计算有个很有意思的地方,在规模在和时,所用时间相比之前低了非常多。但是导出结果分析后又几乎与Python给结果一致,所以计算本身没什么问题。试了几次,时间误差也就正负一两秒。所以大概是触发了什么奇怪的优化吧?
大概是五对角矩阵的逆矩阵仍有一定的稀疏性,或是Python求稀疏矩阵逆的迭代法速度过快,python使用逆矩阵法也有很高的速度优势。
Python使用scipy的spsolve看来是触发了对五对角矩阵的优化迭代法。计算耗时的增加相比于矩阵规模的增长几乎可以忽略不计。scipy这个库还是十分靠谱的。
个人感想
Rust做为静态语言,强制函数输入输出的类型和各类静态检查真的太香了。基本只要不手误写错公式,算出来答案都是对的。Python写的时候心里没底,不报错不一定代表结果是想要的结果。
Rust编译器确实给力,通过编译器检查后,只出现了一两个小BUG,很容易定位,十几分钟就修好了。
Python建模大概花了0.5~1秒,而Rust建模时间几乎可以忽略不计。纯Rust的性能还是非常可靠的。Rust离动力学的基础科学计算的距离其实就差了一个稀疏矩阵求解。但这个确实又很难。的库如果能再给力一点支持稀疏矩阵求解那就真的太香了。
目前的生态来看,python还是科研的首选。
领取专属 10元无门槛券
私享最新 技术干货