1+1=10

记记笔记，放松一下...

矩阵的计算程序与库小记

接前面的矩阵知识小记，稍微整理一下相关的计算程序或库（不考虑商业版）。

对比

涉及矩阵的各种程序和库有很多，暂时关注一下如下4个：

  • Maxima 是一个计算机代数系统，主要用于符号计算。
  • NumPy 是 Python 中的科学计算库，广泛用于数值计算。
  • OpenCV 是 C++ 用于计算机视觉的库，支持一些矩阵操作，但主要用于图像处理。他有python绑定，但python绑定的矩阵操作主要使用NumPy。
  • Eigen 是 C++ 的矩阵和线性代数库，广泛用于数值计算。

矩阵常规操作

矩阵的创建，矩阵的运算，矩阵的分解...

操作 Maxima NumPy OpenCV (C++) Eigen
矩阵创建 matrix([a,b], [c,d]) np.array([[a, b], [c, d]]) cv::Mat A = (cv::Mat_<double>(2,2) << a, b, c, d); Eigen::MatrixXd A(2, 2); A << a, b, c, d;
矩阵加法 A + B A + B cv::add(A, B, C); A + B
矩阵乘法 A . B np.dot(A, B) 或 A @ B C = A * B; 或 cv::gemm(A, B, 1.0, cv::Mat(), 0.0, C); A * B
矩阵转置 transpose(A) A.T cv::transpose(A, B); A.transpose()
矩阵行列式 determinant(A) np.linalg.det(A) cv::determinant(A); A.determinant()
矩阵求逆 invert(A) np.linalg.inv(A) cv::invert(A, B); A.inverse()
矩阵求迹 trace(A) np.trace(A) cv::trace(A); A.trace()
矩阵的LU分解 lu_factor(A) scipy.linalg.lu(A) cv::LU(A, B); (部分功能) Eigen::FullPivLU
矩阵的QR分解 N/A np.linalg.qr(A) N/A Eigen::HouseholderQR
矩阵的特征值分解 eigenvalues(A) np.linalg.eig(A) cv::eigenNonSymmetric() Eigen::EigenSolver
矩阵的SVD分解 N/A np.linalg.svd(A) cv::SVDecomp(A, W, U, Vt); Eigen::JacobiSVD
矩阵与标量相乘 k * A k * A A = A * k; k * A
矩阵拼接 augment(A, B)（水平），stackmatrix(A, B)（竖直） np.hstack((A, B)) (水平)，np.vstack((A, B)) （竖直） cv::hconcat(A, B, C); (水平)，cv::vconcat(A, B, C); （竖直） Eigen::MatrixXd::hstack(A, B) ，Eigen::MatrixXd::vstack(A, B)
矩阵求范数 N/A np.linalg.norm(A) cv::norm(A); A.norm()
单位矩阵 ident(n) np.eye(n) cv::Mat::eye(n, n, CV_64F); Eigen::MatrixXd::Identity(n, n)
零矩阵 zeromatrix(m, n) np.zeros((m, n)) cv::Mat::zeros(m, n, CV_64F); Eigen::MatrixXd::Zero(m, n)
矩阵元素逐项相乘 A .* B A * B cv::multiply(A, B, C); A.cwiseProduct(B)
矩阵元素逐项取平方 N/A np.square(A) cv::pow(A, 2, B); A.array().square()
矩阵对角线元素 diagmatrix(A) np.diag(A) cv::Mat diag = cv::Mat::diag(A); A.diagonal()

存储结构

计算机内存是线性的，即内存地址是一个连续的线性空间，但矩阵是二维（或多维）的数据结构。为了将这些多维数据存储在线性内存中，就需要决定如何对数据进行线性化存储，即按行还是按列的顺序存储。

  • 行主序（Row-major order） 是大多数编程语言（如 C、C++、Python（NumPy 默认））的默认存储方式。这种存储方式更符合线性代数中常见的运算模式，例如矩阵乘法中的行向量和列向量操作。
  • 列主序（Column-major order） 是 Fortran 和一些数学工具（如 MATLAB 和 Eigen（默认））的选择。这是因为 Fortran 早期广泛用于科学计算，使用列主序可以更好地优化向量化操作和循环性能。

比如：对于 2x2 矩阵：

A = \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix}

按 行主序 存储的内存排列为：[1, 2, 3, 4]；按 列主序 存储的内存排列则为：[1, 3, 2, 4]。

库/工具 默认存储顺序 如何选择行/列主序
Maxima 无特别存储顺序 无法选择，列表的列表
NumPy 行主序（可选列主序） 使用 order='C' 选择行主序，<br>使用 order='F' 选择列主序
OpenCV 行主序 无法选择，默认使用行主序
Eigen 列主序（可选行主序） 使用 Eigen::RowMajor 选择行主序，<br>默认列主序

例子-特征值计算

以如下矩阵为例，计算其特征值

A = \begin{pmatrix} 4 & 0 & 1 \\ 2 & 3 & 2 \\ 1 & 0 & 4 \end{pmatrix}

maxima执行的是符号运算，我们可以使用它模拟手动计算过程：

“手动”计算

  • 使用charpoly获取特征多项式，也就是\mathbf{A}-t\mathbf{I}的行列式值
  • 使用factor对特征多项进行分解，得出 -(t-5)*(t-3)^2，可以肉眼看出特征值的结果
  • 使用solve解方程，得出特征值
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
(%i1) A: matrix([4,0,1],[2,3,2],[1,0,4]);
                                  [ 4  0  1 ]
                                  [         ]
(%o1)                             [ 2  3  2 ]
                                  [         ]
                                  [ 1  0  4 ]
(%i2) f: charpoly(A, t);
                                              2
(%o2)                      t + (3 - t) (4 - t)  - 3
(%i3) factor(f);
                                               2
(%o3)                         - (t - 5) (t - 3)
(%i4) solve(f=0, t);
(%o4)                           [t = 5, t = 3]

如何求解特征向量呢？？

比如，接上面，我们要继续求解特征值 3 的特征向量，需要解如下方程 (A - 3I) \mathbf{v} = 0

1
2
(%i12) linsolve([v3+v1=0, 2*v3+2*v1=0, v3+v1=0], [v1,v2,v3]);
(%o12)                 [v1 = - %r1, v2 = %r2, v3 = %r1]

解出来的解是参数形式。据此，将参数用（1,0）和（0,1）带入，可以构造出两个特征向量

1
[1, 0, −1], [0, 1, 0]

同样，对于特征值5，需要解 (A - 5I) \mathbf{v} = 0

1
2
(%i15) linsolve([v3-v1=0, 2*v3-2*v2+2*v1=0, v1-v3=0], [v1,v2,v3]);
(%o15) [v1 = %r3, v2 = 2 %r3, v3 = %r3]

接出来的解依然是参数形式，将参数用1带入，可以构造出特征向量

1
[1,2,1]

maxima自动计算

让 maxima 直接计算特征值也很简单：

1
2
3
4
5
6
7
8
(%i1) A: matrix([4,0,1],[2,3,2],[1,0,4]);
                                  [ 4  0  1 ]
                                  [         ]
(%o1)                             [ 2  3  2 ]
                                  [         ]
                                  [ 1  0  4 ]
(%i2) eigenvalues(A);
(%o2)                          [[5, 3], [1, 2]]

真正的特征值是放在结果的第一个list 当中，也就是5 和3。那第二个list 代表什么呢？代表的就是每个特征值的几何重数, 也就是每个特征值对应的特征向量空间之维度。

如果同时要求出特征向量，可以直接用eigenvectors函数：

1
2
(%i3) eigenvectors(A);
(%o3) [[[5, 3], [1, 2]], [1, 2, 1], [1, 0, −1], [0, 1, 0]]

numpy

Python下功能都是现成的：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
import numpy as np

A = np.array([[4, 0, 1],
              [2, 3, 2],
              [1, 0, 4]])

eigenvalues = np.linalg.eigvals(A)

# A 的特征值： [3. 5. 3.]
print("A 的特征值：", eigenvalues)

如果要同时求出特征向量，可以直接用

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
eigenvalues, eigenvectors = np.linalg.eig(A)

print("A 的特征值：", eigenvalues)
print("A 的特征向量：\n", eigenvectors)

#A 的特征值： [3. 5. 3.]
#A 的特征向量：
# [[ 0.          0.40824829 -0.70710678]
# [ 1.          0.81649658  0.        ]
# [ 0.          0.40824829  0.70710678]]

特征向量要竖着看（列向量）。

Eigen

Eigen下使用 EigenSolver 计算矩阵 A 的特征值和特征向量

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
#include <iostream>
#include <Eigen/Dense>

int main() {
    Eigen::Matrix3d A;
    A << 4, 0, 1,
         2, 3, 2,
         1, 0, 4;

    Eigen::EigenSolver<Eigen::Matrix3d> solver(A);

    std::cout  << solver.eigenvalues() << std::endl;
    std::cout  << solver.eigenvectors() << std::endl;

    return 0;
}

输出结果

1
2
3
4
5
6
(5,0)
(3,0)
(3,0)
 (0.408248,0)  (0.408248,0) (-0.268517,0)
 (0.816497,0) (-0.816497,0) (-0.925093,0)
 (0.408248,0) (-0.408248,0)  (0.268517,0)

特征向量要竖着看（列向量）。

OpenCV

在OpenCV中，使用 cv::eigen 或 cv::eigenNonSymmetric 计算特征值和特征向量

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
#include <opencv2/opencv.hpp>
#include <iostream>

int main() {
    cv::Mat A = (cv::Mat_<double>(3, 3) << 4, 0, 1,
                 2, 3, 2,
                 1, 0, 4);

    cv::Mat eigenvalues, eigenvectors;

    cv::eigenNonSymmetric(A, eigenvalues, eigenvectors);

    std::cout << eigenvalues << std::endl;
    std::cout << eigenvectors << std::endl;

    return 0;
}

输出结果

1
2
3
4
5
6
[5;
 3;
 3]
[0.547722557505166, 1.095445115010332, 0.5477225575051659;
 -0.4082482904638631, 0.8164965809277259, 0.4082482904638629;
 -0.5412318472679982, -1.153604282963793, 0.5412318472679984]

特征向量要竖着看（列向量）。

注意，上面用的cv::eigenNonSymmetric函数，尽管cv::eigen也能用于特征值求解，但是它只适用于对称矩阵。 如果不小心用错了函数，大致会得到下面这种错误的结果：

1
2
3
4
5
6
[5.903211925911553;
 3.80606343352537;
 1.290724640563077]
[0.3971125497870069, 0.5206573684395941, 0.7557893406837773;
 0.8876503388204475, -0.4271322870657471, -0.1721478589408799;
 0.233191978407506, 0.7392387395392244, -0.6317812811178029]

Octave

Octave：Octave 是一个开源的高级编程语言，主要用于数值计算，语法和功能与 MATLAB 高度兼容。有不少在线网站可以用于练习。

在一定程度上，Octave和Python的Numpy可以类比：

特性 Octave Python (NumPy + SciPy)
矩阵支持 原生支持矩阵类型，默认情况下一切操作都是基于矩阵 通过 NumPy 提供矩阵和数组支持，默认是数组（ndarray）
矩阵操作的简洁性 矩阵运算非常简洁，默认处理为二维矩阵，类似 MATLAB NumPy 的数组操作非常灵活，但 Python 本身没有内置矩阵支持
矩阵索引 1 基索引（从 1 开始） 0 基索引（从 0 开始），符合大多数编程语言的习惯
线性代数库 内置矩阵运算与线性代数函数，如 inv、eig、det 等 NumPy 和 SciPy 提供了丰富的线性代数函数，如 inv、eig、det 等
稀疏矩阵支持 支持稀疏矩阵，通过 sparse 函数进行稀疏矩阵运算 SciPy 提供了强大的稀疏矩阵支持，适合处理大规模稀疏矩阵
广播机制 无自动广播机制，矩阵的维度必须匹配 NumPy 支持广播机制，可在不同维度的数组之间进行操作
矩阵乘法 * 表示矩阵乘法，.* 表示元素逐个相乘 @ 表示矩阵乘法，* 表示元素逐个相乘（从 Python 3.5+ 开始）
向量化计算 默认支持向量化计算，所有操作都是矩阵或向量的 NumPy 默认支持向量化计算，避免显式的循环
性能 针对矩阵运算进行了优化，矩阵操作性能较好 NumPy 和 SciPy 的矩阵运算经过高度优化，性能非常好，部分操作可以与 Octave 相媲美
矩阵转置 使用 ' 符号进行矩阵转置 使用 T 属性进行矩阵转置（如 A.T）
高级矩阵分解 内置支持常见的矩阵分解如 QR、LU、SVD SciPy 提供了丰富的矩阵分解功能，如 QR、LU、SVD、Cholesky 等
矩阵求解 使用 A\b 语法求解线性方程组 使用 np.linalg.solve(A, b) 进行线性方程组求解
性能优化 支持与 BLAS/LAPACK 库集成，提升大规模矩阵运算性能 NumPy 和 SciPy 也可以与 BLAS/LAPACK 集成，提供类似的性能优化
GPU 支持 无原生支持，需要额外配置（如借助 OpenCL） 可通过 CuPy 等库轻松扩展到 GPU 支持

其他

  • Armadillo：Armadillo 是一个 C++ 线性代数库，提供类似 MATLAB 的简洁语法，适合科学计算和与 BLAS/LAPACK 集成的高效运算。
  • Blaze：Blaze 是一个高度优化的 C++ 线性代数库，专注于性能，通过自动向量化和并行化提供极致的矩阵运算效率。
  • CuPy 提供了与 NumPy 类似的接口，适合 Python 用户进行 GPU 加速。
  • MAGMA 支持混合 CPU 和 GPU 的矩阵运算，适合需要在多核 CPU 和 GPU 上并行执行的高性能任务。
  • TensorFlow 和 PyTorch 虽然是机器学习框架，但也可以用于 GPU 上的高效矩阵运算。
  • ...

参考

  • https://maxima.sourceforge.io
  • https://numpy.org/
  • https://eigen.tuxfamily.org
  • https://docs.opencv.org
  • https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors
  • https://octave.org/
  • https://arma.sourceforge.net/
  • https://bitbucket.org/blaze-lib/blaze
  • https://octave-online.net/
  • https://www.mycompiler.io/new/octave
  • https://docs.cupy.dev/en/stable/
  • https://icl.utk.edu/magma/

Posted by Debao Zhang Tools