Rust的一些科學計算相關經驗（稀疏矩陣計算的相關生态仍有很大欠缺）

大家好，之前在論壇裡問了不少有關線性代數計算庫的問題，現在姑且來交個作業，順便給出一些用

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 在求解大型線性方程組（系數為稀疏矩陣時）時仍有碾壓性的優勢。

這次我建構的摸準模型是微分方程的隐式動力學求解，差分格式使用的是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'}

。且

的值在每個時間步上需要用多個矩陣進行計算并求解。矩陣尺寸由模型分解出的單元數量決定。

Rust開了優化。Python使用scipy庫。

第一個方法由 sprs 接管所有的矩陣計算，在計算 [K']{u} = {F'} 時将所有矩陣轉化為 sparse21 的矩陣格式計算完後再轉化回 sprs 的矩陣格式。每個時間步都需要來回轉化，來回各一次，是以是2000次。
由于最開始使用了 ndarray 和 sprs 導緻積重難返。又不打算重構。是以沒有純 nalgebra 的實作。方法2的 Rust(ndarray + nalgebra) 意思為，所有計算由 ndarray 實作，除了在計算逆矩陣時。計算逆矩陣時先轉化為 nalgebra 的 DMatrix 并求逆，結果再轉化回 ndarray 的矩陣格式。逆矩陣在整個過程中隻計算一次。是以隻需要來回轉化一輪，來回各一次。
計算逆矩陣，然後每次都是用它計算 {u} = inv([K']){F'} 。
直接求解 [K']{u} = {F'} 。求解1000次。

顯然轉化為密集矩陣的方法在矩陣規模提高之後所使用的時間是不可接受的。但對于密集度在0.5%以上的矩陣，無論時間步的數量有多少，都有一定的實用價值。
sparse21 基本是個玩具庫。應該沒有對矩陣形式進行任何分類求解，而是用的通用方法。但它的計算有個很有意思的地方，在規模在 3499 和 3999 時，所用時間相比之前低了非常多。但是導出結果分析後又幾乎與Python給結果一緻，是以計算本身沒什麼問題。試了幾次，時間誤差也就正負一兩秒。是以大概是觸發了什麼奇怪的優化吧？
大概是五對角矩陣的逆矩陣仍有一定的稀疏性，或是Python求稀疏矩陣逆的疊代法速度過快，python使用逆矩陣法也有很高的速度優勢。
Python使用scipy的spsolve看來是觸發了對五對角矩陣的優化疊代法。計算耗時的增加相比于矩陣規模的增長幾乎可以忽略不計。scipy這個庫還是十分靠譜的。

Rust做為靜态語言，強制函數輸入輸出的類型和各類靜态檢查真的太香了。基本隻要不手誤寫錯公式，算出來答案都是對的。Python寫的時候心裡沒底，不報錯不一定代表結果是想要的結果。
Rust編譯器确實給力，通過編譯器檢查後，隻出現了一兩個小BUG，很容易定位，十幾分鐘就修好了。
Python模組化大概花了0.5~1秒，而Rust模組化時間幾乎可以忽略不計。純Rust的性能還是非常可靠的。Rust離動力學的基礎科學計算的距離其實就差了一個稀疏矩陣求解 Ax=B 。但這個确實又很難。 nalgebra 的庫如果能再給力一點支援稀疏矩陣求解那就真的太香了。
目前的生态來看，python還是科研的首選。