矩阵的计算程序与库小记

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

对比

涉及矩阵的各种程序和库有很多，暂时关注一下如下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解方程，得出特征值

(%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

(%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

(%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,2,1]

maxima自动计算

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

(%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函数：

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

numpy

Python下功能都是现成的：

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)

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

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 的特征值和特征向量

#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;
}

输出结果

(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 计算特征值和特征向量

#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;
}

输出结果

[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也能用于特征值求解，但是它只适用于对称矩阵。如果不小心用错了函数，大致会得到下面这种错误的结果：

[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/

1+1=10

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

矩阵的计算程序与库小记

对比

矩阵常规操作

存储结构

例子-特征值计算

“手动”计算

maxima自动计算

numpy

Eigen

OpenCV

Octave

其他

参考