折腾笔记[5]-使用rust解矩阵方程
摘要
使用rust的nalgebra库解矩阵方程,稠密矩阵的代数运算(逆,特征值等).
关键信息
[package]
name = "exp65-rust-ziglang-slambook2"
version = "0.1.0"
edition = "2021"
[dependencies]
nalgebra = { version = "0.33.2",features = ["rand"]}
rand = "0.8.5"
原理简介
rust的类似Eigen线代库
[https://rustcc.cn/article?id=15a3a248-d6da-4094-8bc9-3411faf79668]
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 在求解大型线性方程组（系数为稀疏矩阵时）时仍有碾压性的优势。
C++ 的eigen/sparse库或许是比 scipy 更优的选择。而且比较方便 FFI 至 Rust 使用。
nalgebra库简介
[https://nalgebra.org/docs/user_guide/vectors_and_matrices]
nalgebra 是一个使用 Rust 编程语言编写的线性代数库。它适用于需要高效向量、矩阵操作以及几何变换的应用程序。该库旨在实现快速且全面的功能。
nalgebra 是免费且开源的，遵循 Apache 2.0 许可证发布。
nalgebra is a linear algebra library written using the Rust programming language. It targets applications requiring efficient vector and matrix operations, as well as geometric transformations. It is designed to be fast and comprehensive.
nalgebra is free and open-source, released under the Apache 2.0 license. It is developed by the Dimforge open-source company. You can support us by sponsoring us on GitHub sponsor.
解矩阵方程几种方式简介
[https://szup.github.io/2021/02/21/0221-ju-zhen-yu-xian-xing-fang-cheng-zu/]
1. 直接求逆
直接求逆法适用于系数矩阵是方阵且可逆的情况。给定一个线性方程组 \(Ax = b\)，其中 \(A\) 是一个 \(n \times n\) 的可逆矩阵，\(x\) 和 \(b\) 是 \(n\) 维列向量。解这个方程的一种方法是直接计算矩阵 \(A\) 的逆矩阵 \(A^{-1}\)，然后通过左乘 \(A^{-1}\) 来求解 \(x\)，即：
\[x = A^{-1}b
\]
优点：计算过程简单直接。
缺点：
- 当矩阵 \(A\) 很大时，计算逆矩阵非常耗时。
- 如果 \(A\) 接近奇异或者条件数很大，计算逆矩阵可能不准确。
2. QR分解方式
QR分解将矩阵 \(A\) 分解为一个正交矩阵 \(Q\) 和一个上三角矩阵 \(R\)，即 \(A = QR\)。解方程 \(Ax = b\) 的步骤如下：
- 计算 \(A\) 的 QR 分解。
- 将方程 \(Ax = b\) 转化为 \(Rx = Q^Tb\)。
- 解上三角方程组 \(Rx = Q^Tb\) 得到 \(x\)。
优点：
- 不需要计算逆矩阵，因此比直接求逆法更稳定。
- 适用于非方阵和方阵。
缺点： - 计算复杂度较高，尤其是对于大型矩阵。
3. Cholesky 分解方式
Cholesky 分解是 QR 分解的一个特例，它只适用于对称正定矩阵 \(A\)。分解将 \(A\) 表示为 \(A = LL^T\)，其中 \(L\) 是一个下三角矩阵。解方程 \(Ax = b\) 的步骤如下：
- 计算 \(A\) 的 Cholesky 分解 \(A = LL^T\)。
- 解两个三角方程组：先解 \(Ly = b\) 得到 \(y\)，然后解 \(L^Tx = y\) 得到 \(x\)。
优点：
- 计算效率高，比 QR 分解通常要快。
- 精度高，因为只涉及实数运算。
缺点： - 只适用于对称正定矩阵。
每种方法都有其适用场景和优缺点，选择哪种方法取决于具体问题的性质以及系数矩阵的特点。在实际应用中，还可能使用迭代方法等其他数值解法。
实现
原始C++代码:
[https://github.com/gaoxiang12/slambook2/blob/master/ch3/useEigen/eigenMatrix.cpp]
#include <iostream>
using namespace std;
#include <ctime>
// Eigen 核心部分
#include <Eigen/Core>
// 稠密矩阵的代数运算（逆，特征值等）
#include <Eigen/Dense>
using namespace Eigen;
#define MATRIX_SIZE 50
/****************************
* 本程序演示了 Eigen 基本类型的使用
****************************/
int main(int argc, char **argv) {
// Eigen 中所有向量和矩阵都是Eigen::Matrix，它是一个模板类。它的前三个参数为：数据类型，行，列
// 声明一个2*3的float矩阵
Matrix<float, 2, 3> matrix_23;
// 同时，Eigen 通过 typedef 提供了许多内置类型，不过底层仍是Eigen::Matrix
// 例如 Vector3d 实质上是 Eigen::Matrix<double, 3, 1>，即三维向量
Vector3d v_3d;
// 这是一样的
Matrix<float, 3, 1> vd_3d;
// Matrix3d 实质上是 Eigen::Matrix<double, 3, 3>
Matrix3d matrix_33 = Matrix3d::Zero(); //初始化为零
// 如果不确定矩阵大小，可以使用动态大小的矩阵
Matrix<double, Dynamic, Dynamic> matrix_dynamic;
// 更简单的
MatrixXd matrix_x;
// 这种类型还有很多，我们不一一列举
// 下面是对Eigen阵的操作
// 输入数据（初始化）
matrix_23 << 1, 2, 3, 4, 5, 6;
// 输出
cout << "matrix 2x3 from 1 to 6: \n" << matrix_23 << endl;
// 用()访问矩阵中的元素
cout << "print matrix 2x3: " << endl;
for (int i = 0; i < 2; i++) {
for (int j = 0; j < 3; j++) cout << matrix_23(i, j) << "\t";
cout << endl;
}
// 矩阵和向量相乘（实际上仍是矩阵和矩阵）
v_3d << 3, 2, 1;
vd_3d << 4, 5, 6;
// 但是在Eigen里你不能混合两种不同类型的矩阵，像这样是错的
// Matrix<double, 2, 1> result_wrong_type = matrix_23 * v_3d;
// 应该显式转换
Matrix<double, 2, 1> result = matrix_23.cast<double>() * v_3d;
cout << "[1,2,3;4,5,6]*[3,2,1]=" << result.transpose() << endl;
Matrix<float, 2, 1> result2 = matrix_23 * vd_3d;
cout << "[1,2,3;4,5,6]*[4,5,6]: " << result2.transpose() << endl;
// 同样你不能搞错矩阵的维度
// 试着取消下面的注释，看看Eigen会报什么错
// Eigen::Matrix<double, 2, 3> result_wrong_dimension = matrix_23.cast<double>() * v_3d;
// 一些矩阵运算
// 四则运算就不演示了，直接用+-*/即可。
matrix_33 = Matrix3d::Random(); // 随机数矩阵
cout << "random matrix: \n" << matrix_33 << endl;
cout << "transpose: \n" << matrix_33.transpose() << endl; // 转置
cout << "sum: " << matrix_33.sum() << endl; // 各元素和
cout << "trace: " << matrix_33.trace() << endl; // 迹
cout << "times 10: \n" << 10 * matrix_33 << endl; // 数乘
cout << "inverse: \n" << matrix_33.inverse() << endl; // 逆
cout << "det: " << matrix_33.determinant() << endl; // 行列式
// 特征值
// 实对称矩阵可以保证对角化成功
SelfAdjointEigenSolver<Matrix3d> eigen_solver(matrix_33.transpose() * matrix_33);
cout << "Eigen values = \n" << eigen_solver.eigenvalues() << endl;
cout << "Eigen vectors = \n" << eigen_solver.eigenvectors() << endl;
// 解方程
// 我们求解 matrix_NN * x = v_Nd 这个方程
// N的大小在前边的宏里定义，它由随机数生成
// 直接求逆自然是最直接的，但是求逆运算量大
Matrix<double, MATRIX_SIZE, MATRIX_SIZE> matrix_NN
= MatrixXd::Random(MATRIX_SIZE, MATRIX_SIZE);
matrix_NN = matrix_NN * matrix_NN.transpose(); // 保证半正定
Matrix<double, MATRIX_SIZE, 1> v_Nd = MatrixXd::Random(MATRIX_SIZE, 1);
clock_t time_stt = clock(); // 计时
// 直接求逆
Matrix<double, MATRIX_SIZE, 1> x = matrix_NN.inverse() * v_Nd;
cout << "time of normal inverse is "
<< 1000 * (clock() - time_stt) / (double) CLOCKS_PER_SEC << "ms" << endl;
cout << "x = " << x.transpose() << endl;
// 通常用矩阵分解来求，例如QR分解，速度会快很多
time_stt = clock();
x = matrix_NN.colPivHouseholderQr().solve(v_Nd);
cout << "time of Qr decomposition is "
<< 1000 * (clock() - time_stt) / (double) CLOCKS_PER_SEC << "ms" << endl;
cout << "x = " << x.transpose() << endl;
// 对于正定矩阵，还可以用cholesky分解来解方程
time_stt = clock();
x = matrix_NN.ldlt().solve(v_Nd);
cout << "time of ldlt decomposition is "
<< 1000 * (clock() - time_stt) / (double) CLOCKS_PER_SEC << "ms" << endl;
cout << "x = " << x.transpose() << endl;
return 0;
}
使用nalgebra重构代码:
// ch3-useNalgebra
// 稠密矩阵的代数运算(逆,特征值等)
// 测试Matrix,Vector3,Vector2,OVector基本类型的使用
extern crate nalgebra as na;
extern crate rand; // 确保在Cargo.toml中添加了rand依赖
use rand::Rng; // 引入Rng trait以使用随机数生成功能
use na::{Matrix3, Matrix, Vector3, OVector, Dyn, SymmetricEigen, SMatrix, Const};
use std::time::Instant;
// 矩阵大小
const MATRIX_SIZE: usize = 50;
// examples单文件都要有main函数
fn main() {
// 声明一个 2x3 的 f32 矩阵
let matrix_23 = na::Matrix2x3::<f32>::new(1.0, 2.0, 3.0, 4.0, 5.0, 6.0);
println!("matrix 2x3 from 1 to 6:\n{}", matrix_23);
// 用索引访问矩阵中的元素
println!("print matrix 2x3:");
for i in 0..2 {
for j in 0..3 {
print!("{}\t", matrix_23[(i, j)]);
}
println!();
}
// 矩阵和向量相乘
let v_3d = Vector3::new(3.0, 2.0, 1.0);
let vd_3d = Vector3::new(4.0, 5.0, 6.0);
// 显式转换
let result: na::Vector2<f64> = matrix_23.map(|x| x as f64) * v_3d;
println!("[1,2,3;4,5,6]*[3,2,1]={}", result.transpose());
let result2 = matrix_23 * vd_3d;
println!("[1,2,3;4,5,6]*[4,5,6]: {}", result2.transpose());
// 生成3阶整数随机方阵
// pub fn new_random(nrows: usize, ncols: usize) -> Self
// Matrix3为3阶方阵
// let matrix_33 = Matrix3::<f64>::new_random(); // 生成的全是小数
let matrix_33 = Matrix3::from_fn(|_, _| rand::random::<i32>() as f64);
println!("random matrix:\n{}", matrix_33);
// 矩阵转置
println!("transpose:\n{}", matrix_33.transpose());
// 矩阵各元素之和
println!("sum: {}", matrix_33.sum());
// 矩阵的迹
println!("trace: {}", matrix_33.trace());
// 矩阵倍乘
println!("times 10:\n{}", 10.0 * matrix_33);
// 矩阵求逆
println!("inverse:\n{}", matrix_33.try_inverse().unwrap());
// 矩阵行列式
println!("det: {}", matrix_33.determinant());
// 特征值
// 实对称矩阵可以保证对角化成功
let eigen_solver = matrix_33.clone().transpose() * matrix_33.clone();
// 计算对称矩阵的特征值和特征向量
let eigen = eigen_solver.symmetric_eigen();
// 获取特征值和特征向量
let eigen_values = eigen.eigenvalues;
let eigen_vectors = eigen.eigenvectors;
println!("Eigen values = \n{:?}", eigen_values);
println!("Eigen vectors = \n{}", eigen_vectors);
// 解方程
// 求解matrix_NN * x = v_Nd方程
let mut matrix_NN = SMatrix::<f64, MATRIX_SIZE, MATRIX_SIZE>::from_fn(|_, _| rand::random::<i32>() as f64);
matrix_NN = matrix_NN * matrix_NN.transpose(); // 保证半正定
let v_Nd = OVector::<f64, Const<MATRIX_SIZE> >::new_random();
// 记录开始时间
let start_time = Instant::now();
// 直接求逆方式解方程
let x = matrix_NN.try_inverse().unwrap() * v_Nd;
println!("time of normal inverse is {}ms", start_time.elapsed().as_millis());
println!("x = {}", x.transpose());
// QR 分解方式解方程
let start_time = Instant::now();
let x = matrix_NN.qr().solve(&v_Nd).unwrap();
println!("time of QR decomposition is {}ms", start_time.elapsed().as_millis());
println!("x = {}", x.transpose());
// Cholesky 分解方式解方程
let start_time = Instant::now();
let x = matrix_NN.cholesky().unwrap().solve(&v_Nd);
println!("time of Cholesky decomposition is {}ms", start_time.elapsed().as_millis());
println!("x = {}", x.transpose());
}
效果
matrix 2x3 from 1 to 6:
┌ ┐
│ 1 2 3 │
│ 4 5 6 │
└ ┘
print matrix 2x3:
1 2 3
4 5 6
[1,2,3;4,5,6]*[3,2,1]=
┌ ┐
│ 10 28 │
└ ┘
[1,2,3;4,5,6]*[4,5,6]:
┌ ┐
│ 32 77 │
└ ┘
random matrix:
┌ ┐
│ 109250093 -668874409 654624440 │
│ -1807777515 -1241728103 -72709058 │
│ -1978173100 -1194711519 713117762 │
└ ┘
transpose:
┌ ┐
│ 109250093 -1807777515 -1978173100 │
│ -668874409 -1241728103 -1194711519 │
│ 654624440 -72709058 713117762 │
└ ┘
sum: -5486981409
trace: -419360248
times 10:
┌ ┐
│ 1092500930 -6688744090 6546244400 │
│ -18077775150 -12417281030 -727090580 │
│ -19781731000 -11947115190 7131177620 │
└ ┘
inverse:
┌ ┐
│ 0.0000000007724109318863737 0.0000000002423611740513393 -0.0000000006843431016305952 │
│ -0.000000001138314285420112 -0.0000000010905565874484238 0.0000000009337518218597835 │
│ 0.00000000023559269944238927 -0.0000000011547435865903803 0.0000000010682877124113894 │
└ ┘
det: -1258869695923911800000000000
Eigen values =
[[1.0509286097505001e19, 8.694050660616129e17, 1.7344675402389238e17]]
Eigen vectors =
┌ ┐
│ 0.8164472586801655 -0.4321788280858501 0.38292993399313163 │
│ 0.551153645539266 0.38554199648086795 -0.7399912350550671 │
│ -0.17217297349590083 -0.8152170444144846 -0.5529716427572867 │
└ ┘
time of normal inverse is 9ms
x =
┌ ┐
│ -0.0000000000000000015975491572516443 0.0000000000000000024747949268696407 -0.00000000000000000020437622451601015 -0.0000000000000000010414872227270577 0.00000000000000000029336960867197514 0.0000000000000000006126855626882509 -0.0000000000000000007575489972054902 -0.000000000000000000048833514427405895 0.00000000000000000032336382003806876 0.0000000000000000021711467586520475 -0.0000000000000000006126828266221918 -0.00000000000000000018806148861941312 0.000000000000000001250381959129717 -0.0000000000000000009231889760824598 -0.0000000000000000013415945337412047 -0.00000000000000000026767849668613976 0.0000000000000000009474064254082232 0.0000000000000000019723365476720838 -0.00000000000000000011561140190105468 0.0000000000000000004269021003751211 -0.0000000000000000010668971005859026 -0.0000000000000000019310850106422213 0.0000000000000000001729842627249784 -0.0000000000000000005254364541948351 -0.0000000000000000012829167572273559 0.0000000000000000012094907424250509 0.00000000000000000015103978456218436 0.0000000000000000011411146258393528 0.000000000000000001838145304205453 -0.000000000000000000966553680421858 0.0000000000000000014756282913129499 -0.0000000000000000017919398199453582 -0.00000000000000000009657807621849883 -0.0000000000000000004953694581313736 -0.00000000000000000017310632032040021 0.000000000000000003147443552706895 0.000000000000000001586089875584719 0.000000000000000001553159674141187 -0.0000000000000000015899748634351808 -0.000000000000000002345156993639221 0.0000000000000000005563889636772475 -0.0000000000000000011139953424593626 -0.00000000000000000013924429591404612 0.0000000000000000004650761368431657 0.0000000000000000021109645923366945 -0.0000000000000000003178292744147794 -0.00000000000000000011142539902368143 0.0000000000000000008642475451802308 0.00000000000000000031590499902692005 -0.0000000000000000010927512484315357 │
└ ┘
time of QR decomposition is 4ms
x =
┌ ┐
│ -0.0000000000000000015975491572518145 0.0000000000000000024747949268699423 -0.0000000000000000002043762245160872 -0.000000000000000001041487222727218 0.0000000000000000002933696086719844 0.0000000000000000006126855626882805 -0.0000000000000000007575489972056011 -0.000000000000000000048833514427440544 0.0000000000000000003233638200380944 0.000000000000000002171146758652289 -0.0000000000000000006126828266223485 -0.00000000000000000018806148861950063 0.0000000000000000012503819591299012 -0.000000000000000000923188976082632 -0.0000000000000000013415945337414408 -0.00000000000000000026767849668611935 0.0000000000000000009474064254083964 0.000000000000000001972336547672298 -0.00000000000000000011561140190111655 0.0000000000000000004269021003751309 -0.0000000000000000010668971005859839 -0.000000000000000001931085010642424 0.00000000000000000017298426272494693 -0.0000000000000000005254364541949053 -0.0000000000000000012829167572274816 0.0000000000000000012094907424251909 0.00000000000000000015103978456218205 0.0000000000000000011411146258394686 0.0000000000000000018381453042056424 -0.0000000000000000009665536804219449 0.0000000000000000014756282913130739 -0.0000000000000000017919398199455304 -0.00000000000000000009657807621851155 -0.0000000000000000004953694581313966 -0.00000000000000000017310632032046158 0.0000000000000000031474435527073572 0.0000000000000000015860898755850427 0.0000000000000000015531596741413143 -0.000000000000000001589974863435349 -0.0000000000000000023451569936394668 0.0000000000000000005563889636772757 -0.0000000000000000011139953424593934 -0.0000000000000000001392442959140611 0.0000000000000000004650761368431138 0.000000000000000002110964592336955 -0.00000000000000000031782927441485634 -0.00000000000000000011142539902367284 0.0000000000000000008642475451803701 0.00000000000000000031590499902690517 -0.0000000000000000010927512484316682 │
└ ┘
time of Cholesky decomposition is 1ms
x =
┌ ┐
│ -0.0000000000000000015975491572516027 0.0000000000000000024747949268695695 -0.000000000000000000204376224515967 -0.0000000000000000010414872227270263 0.00000000000000000029336960867197355 0.0000000000000000006126855626882356 -0.0000000000000000007575489972054472 -0.000000000000000000048833514427396904 0.0000000000000000003233638200380565 0.0000000000000000021711467586519655 -0.0000000000000000006126828266221657 -0.0000000000000000001880614886193993 0.0000000000000000012503819591296465 -0.0000000000000000009231889760823843 -0.0000000000000000013415945337411213 -0.0000000000000000002676784966861455 0.0000000000000000009474064254081526 0.0000000000000000019723365476720206 -0.00000000000000000011561140190101994 0.00000000000000000042690210037511846 -0.000000000000000001066897100585899 -0.0000000000000000019310850106421555 0.0000000000000000001729842627250054 -0.000000000000000000525436454194817 -0.0000000000000000012829167572273449 0.0000000000000000012094907424249956 0.00000000000000000015103978456218205 0.0000000000000000011411146258393295 0.0000000000000000018381453042053955 -0.0000000000000000009665536804218216 0.0000000000000000014756282913129221 -0.0000000000000000017919398199453212 -0.00000000000000000009657807621851943 -0.0000000000000000004953694581313577 -0.00000000000000000017310632032039097 0.00000000000000000314744355270675 0.000000000000000001586089875584638 0.000000000000000001553159674141151 -0.0000000000000000015899748634351375 -0.000000000000000002345156993639138 0.0000000000000000005563889636772334 -0.000000000000000001113995342459331 -0.00000000000000000013924429591407496 0.000000000000000000465076136843169 0.0000000000000000021109645923366117 -0.00000000000000000031782927441476857 -0.00000000000000000011142539902369797 0.0000000000000000008642475451801823 0.0000000000000000003159049990269215 -0.0000000000000000010927512484315045 │
└ ┘
【推荐】国内首个AI IDE，深度理解中文开发场景，立即下载体验Trae
【推荐】编程新体验，更懂你的AI，立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包，你的智能百科全书，全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell：AI 加持，快人一步
· 全程不用写代码，我用AI程序员写了一个飞机大战
· MongoDB 8.0这个新功能碉堡了，比商业数据库还牛
· 记一次.NET内存居高不下排查解决与启示
· 白话解读 Dapr 1.15：你的「微服务管家」又秀新绝活了
· DeepSeek 开源周回顾「GitHub 热点速览」