矩阵线性代数¶

算法流程¶

Gauss-Jordan 消元法求逆¶

目标：计算方阵 \(A\) 的逆矩阵 \(A^{-1}\)。

算法：

1. 构造增广矩阵 \([A \mid I_n]\)

2. 对每一列 \(k = 0, \dots, n-1\)：

3. 选主元：在 \([A \mid I]\) 的第 \(k\) 列 \(k..n-1\) 行中找绝对值最大的元素

4. 换行：将该行交换到第 \(k\) 行
5. 归一化：将第 \(k\) 行除以主元值，使 \(A_{kk} = 1\)
6. 消元：对所有其他行 \(i \neq k\)，用第 \(k\) 行消去第 \(i\) 行的第 \(k\) 列
7. 左半部分变为 \(I\)，右半部分即为 \(A^{-1}\)

复杂度：\(O(n^3)\)。比伴随矩阵法（\(O(n^4)\)）快得多。

注意：

• 若主元接近零，归一化会放大舍入误差
• 对多个右端项，优先用 LU 分解 + 前向/回代，而非显式求逆

带状矩阵求解（Thomas 算法）¶

目标：高效求解三对角系统 \(A \cdot x = b\)，其中 \(A\) 仅有主对角线及相邻两条对角线非零。

算法（Thomas 算法，三对角特例）：

1. 前向消元（\(i = 1, \dots, n-1\)）：

2. \(w = A_{i,i-1} / A_{i-1,i-1}\)

3. \(A_{i,i} \;\;{-}{=}\;\; w \cdot A_{i-1,i}\)
4. \(b_i \;\;{-}{=}\;\; w \cdot b_{i-1}\)

5. 回代（\(i = n-1, \dots, 0\)）：

6. \(x_i = (b_i - A_{i,i+1} \cdot x_{i+1}) / A_{i,i}\)

复杂度：\(O(n)\) vs 普通高斯消元 \(O(n^3)\)。

适用条件：矩阵必须对角占优或对称正定以保证数值稳定。

线性代数¶

转置矩阵¶

Mat Mat::transpose();

描述:

计算矩阵的转置，返回新矩阵。矩阵A的转置A^T 通过交换行和列得到：(A^T)ᵢⱼ = Aⱼᵢ。

数学原理:

• 对于任何矩阵A，(A^T) ^T = A

• (A + B)^T = A^T + B^T

• (AB)^T = B^T * A^T

• 对于方阵，det(A) = det(A^T)

参数:

无。

返回值:

Mat - 转置后的矩阵 (col × row)。

使用建议:

• 内存布局: 创建新矩阵，因此内存使用量暂时翻倍。对于大矩阵，考虑内存限制。

• 对称矩阵: 如果A = A^T，矩阵是对称的。使用is_symmetric()检查。

• 应用:

• 内积: u^T * v

• 二次型: x^T * A * x

• 矩阵方程: A^T * A（正规方程）

余子式矩阵¶

Mat Mat::minor(int target_row, int target_col);

描述:

通过移除指定的行和列来计算余子式矩阵。余子式是移除一行一列后得到的子矩阵。

参数:

• int target_row: 要移除的行索引

• int target_col: 要移除的列索引

返回值:

Mat - (n-1)x(n-1) 的余子式矩阵

代数余子式矩阵¶

Mat Mat::cofactor(int target_row, int target_col);

描述:

计算代数余子式矩阵（与余子式矩阵相同）。代数余子式矩阵与余子式矩阵相同。符号 (-1)^(i+j) 在计算代数余子式值时应用，而不是应用到矩阵元素本身。

参数:

• int target_row: 要移除的行索引

• int target_col: 要移除的列索引

返回值:

Mat - (n-1)x(n-1) 的代数余子式矩阵（与余子式矩阵相同）

行列式（自动选择方法）¶

float Mat::determinant();

描述:

计算方阵的行列式，根据矩阵大小自动选择最优方法。对于小矩阵（n ≤ 4），使用拉普拉斯展开法；对于大矩阵（n > 4），使用LU分解法以提高效率。

数学原理:

行列式是方阵的一个重要数值特征，具有以下性质：

• det(AB) = det(A) * det(B)

• det(A^T) = det(A)

• det(A^(-1)) = 1 / det(A)

• 如果矩阵是奇异的，det(A) = 0

方法选择:

• 小矩阵 (n ≤ 4): 使用 determinant_laplace() - 拉普拉斯展开法，时间复杂度 O(n!)，对小矩阵更准确

• 大矩阵 (n > 4): 使用 determinant_lu() - LU分解法，时间复杂度 O(n³)，效率更高

参数:

void

返回值:

float - 行列式的值

使用建议:

• 自动选择: 对于大多数应用，直接使用 determinant() 即可，函数会自动选择最优方法

• 性能优化: 如果需要多次计算相同大小的矩阵行列式，可以考虑直接调用 determinant_lu() 或 determinant_gaussian()

• 精度要求: 对于小矩阵，determinant_laplace() 可能提供更好的数值精度

行列式 - 拉普拉斯展开法¶

float Mat::determinant_laplace();

描述:

使用拉普拉斯展开（余子式展开）计算方阵的行列式。时间复杂度为 O(n!)，仅适用于小矩阵（n ≤ 4）。

数学原理:

拉普拉斯展开是行列式的递归定义方法：

• 对于 1×1 矩阵: det([a]) = a

• 对于 2×2 矩阵: det([[a,b],[c,d]]) = ad - bc

• 对于 n×n 矩阵: det(A) = Σⱼ₌₁ⁿ (-1)ⁱ⁺ʲ aᵢⱼ * det(Mᵢⱼ)，其中 Mᵢⱼ 是余子式矩阵

本实现使用第一行展开，递归计算余子式的行列式。

参数:

void

返回值:

float - 行列式的值

行列式 - LU分解法¶

float Mat::determinant_lu();

描述:

使用LU分解计算方阵的行列式。时间复杂度为 O(n³)，适用于大矩阵。

数学原理:

LU分解将矩阵分解为 A = P * L * U，其中：

• P 是置换矩阵（如果使用主元）

• L 是单位对角线的下三角矩阵

• U 是上三角矩阵

行列式计算公式：det(A) = det(P) * det(L) * det(U)

其中：

• det(P) = (-1)^(置换的符号)，由行交换次数决定

• det(L) = 1（因为L是单位对角线的下三角矩阵）

• det(U) = ∏ᵢ Uᵢᵢ（U的对角线元素的乘积）

算法步骤:

1. 执行LU分解（带主元以提高数值稳定性）
2. 计算置换矩阵的行列式 det(P)
3. 计算上三角矩阵U的对角线元素乘积 det(U)
4. 返回 det(P) * det(U)

参数:

void

返回值:

float - 行列式的值。如果矩阵是奇异的或接近奇异的，返回 0.0

使用建议:

• 效率: 对于 n > 4 的矩阵，比拉普拉斯展开法快得多

• 数值稳定性: 使用主元（pivoting）提高数值稳定性

• 奇异矩阵: 如果矩阵是奇异的，LU分解会失败，函数返回 0.0

行列式 - 高斯消元法¶

float Mat::determinant_gaussian();

描述:

使用高斯消元法计算方阵的行列式。时间复杂度为 O(n³)，适用于大矩阵。

数学原理:

高斯消元法将矩阵转换为上三角形式，然后计算对角线元素的乘积。行列式的值等于上三角矩阵对角线元素的乘积，并根据行交换次数调整符号。

算法步骤:

1. 使用部分主元法进行高斯消元，将矩阵转换为上三角形式
2. 跟踪行交换次数
3. 计算上三角矩阵对角线元素的乘积
4. 根据行交换次数调整符号：每次行交换使行列式乘以 -1

参数:

void

返回值:

float - 行列式的值。如果矩阵是奇异的，返回 0.0

使用建议:

• 效率: 对于大矩阵，时间复杂度为 O(n³)，与LU分解法相当

• 数值稳定性: 使用部分主元法提高数值稳定性

• 实现简单: 相比LU分解，实现更直观，但功能较少（不能用于求解线性系统）

伴随矩阵¶

Mat Mat::adjoint();

描述:

计算方阵的伴随（或伴随转置）矩阵。

参数:

void

返回值:

Mat - 伴随矩阵对象

归一化¶

void Mat::normalize();

描述:

使用L2范数（Frobenius范数）归一化矩阵。归一化后，||Matrix|| = 1。

参数:

void

返回值:

void

范数¶

float Mat::norm() const;

描述:

计算矩阵的Frobenius范数（L2范数）。

参数:

void

返回值:

float - 计算得到的矩阵范数

矩阵求逆 -- 基于伴随矩阵¶

Mat Mat::inverse_adjoint();

描述:

使用伴随矩阵法计算方阵的逆矩阵。如果矩阵是奇异的，返回零矩阵。

参数:

void

返回值:

Mat - 逆矩阵对象。如果矩阵是奇异的，返回零矩阵

单位矩阵¶

static Mat Mat::eye(int size);

描述:

生成指定大小的单位矩阵。

参数:

• int size: 方阵的维度

返回值:

Mat - 单位矩阵 (size x size)

增广矩阵（水平连接）¶

static Mat Mat::augment(const Mat &A, const Mat &B);

描述:

通过水平连接两个矩阵创建增广矩阵 [A | B]。A和B的行数必须匹配。

参数:

• const Mat &A: 左侧矩阵

• const Mat &B: 右侧矩阵

返回值:

Mat - 增广矩阵 [A B]

垂直堆叠¶

static Mat Mat::vstack(const Mat &A, const Mat &B);

描述:

垂直堆叠两个矩阵 [A; B]。A和B的列数必须匹配。

参数:

• const Mat &A: 顶部矩阵

• const Mat &B: 底部矩阵

返回值:

Mat - 垂直堆叠的矩阵 [A; B]

Gram-Schmidt正交化¶

static bool Mat::gram_schmidt_orthogonalize(const Mat &vectors, Mat &orthogonal_vectors, 
                                            Mat &coefficients, float tolerance = 1e-6f);

描述:

使用Gram-Schmidt过程对一组向量进行正交化。这是一个通用的正交化函数，可重复用于QR分解和其他需要正交基的应用。使用改进的Gram-Schmidt算法和重新正交化以提高数值稳定性。

数学原理:

给定一组向量 {v₁, v₂, ..., vₙ}，Gram-Schmidt过程产生正交集合 {q₁, q₂, ..., qₙ}，其中：

• q₁ = v₁ / ||v₁||

• qⱼ = (vⱼ - Σᵢ₌₁ʲ⁻¹⟨vⱼ, qᵢ⟩qᵢ) / ||vⱼ - Σᵢ₌₁ʲ⁻¹⟨vⱼ, qᵢ⟩qᵢ||

改进版本立即减去投影，这提高了数值稳定性。

参数:

• const Mat &vectors : 输入矩阵，其中每列是要正交化的向量 (m × n)。

• Mat &orthogonal_vectors : 输出矩阵，用于正交化向量 (m × n)，每列都是正交且归一化的。

• Mat &coefficients : 输出矩阵，用于投影系数 (n × n，上三角)，类似于QR分解中的R矩阵。

• float tolerance : 线性独立性检查的最小范数阈值（默认：1e-6）。

返回值:

bool - 成功返回true，输入无效返回false。

使用建议:

• 数值稳定性: 实现使用改进的Gram-Schmidt和重新正交化，这显著提高了近线性相关向量的稳定性。

• QR分解: 此函数由qr_decompose()内部使用。对于QR分解，系数矩阵对应于R矩阵。

• 基构造: 用于从一组向量构造正交基，这在许多线性代数应用中都是基础。

• 性能: 对于大矩阵，考虑计算成本。对于m维向量和n个向量，复杂度为O(mn²)。

全1矩阵（矩形）¶

static Mat Mat::ones(int rows, int cols);

描述:

创建指定大小的全1矩阵。

参数:

• int rows: 行数

• int cols: 列数

返回值:

Mat - 矩阵 [rows x cols]，所有元素 = 1

全1矩阵（方阵）¶

static Mat Mat::ones(int size);

描述:

创建指定大小的方阵，所有元素为1。

参数:

• int size: 方阵的大小（行数 = 列数）

返回值:

Mat - 方阵 [size x size]，所有元素 = 1

高斯消元法¶

Mat Mat::gaussian_eliminate() const;

描述:

执行高斯消元法，将矩阵转换为行阶梯形式（REF）。这是求解线性系统和计算矩阵秩的第一步。

数学原理:

高斯消元通过初等行变换将矩阵转换为行阶梯形式：

1. 行交换: 交换两行

2. 行缩放: 将一行乘以非零标量

3. 行加法: 将一行的倍数加到另一行

行阶梯形式（REF）的性质:

• 所有零行在底部

• 每个非零行的前导系数（主元）位于上一行主元的右侧

• 主元下方的所有条目为零

参数:

无。

返回值:

Mat - 上三角矩阵（REF形式）。

使用建议:

• 线性系统求解: 求解Ax = b的第一步。REF之后，使用回代。

• 秩计算: 秩等于REF中非零行的数量。

• 行列式: 可以从REF计算行列式（对角元素的乘积，根据行交换进行调整）。

• 数值稳定性: 实现使用部分主元以提高数值稳定性。

• 性能: 对于n×n矩阵为O(n³)。对于多个系统，优先使用LU分解。

坑点:

• 接近零的主元: 如果主元非常小，消元会放大舍入误差。部分主元只能缓解，不能彻底消除。

• 奇异或秩亏矩阵: 你可能仍然得到一个“看起来正常”的REF，但这不代表方程组有唯一解。

• 行列式符号: 每做一次行交换，行列式符号就会翻转。若用REF估算det(A)，必须把交换次数算进去。

• 工作副本视角: 算法改变的是等价的工作副本，不是数学意义上的原矩阵对象。不要假设原始行顺序会保留。

从高斯消元到行最简形式¶

Mat Mat::row_reduce_from_gaussian() const;

描述:

将矩阵（假设已为行阶梯形式）转换为简化行阶梯形式（RREF）。

数学原理:

这本质上是高斯消元的回代阶段再加一次规范化：先把每个主元行缩放到主元为1，再把主元上方的元素全部消去，使每一列的主元成为唯一非零项。

使用建议:

• 最好先有REF: 先调用 gaussian_eliminate() 再调用这个函数最可靠。直接对任意矩阵调用有时也能工作，但结果不一定有明确的数值意义。

• 不是独立求解器: RREF更适合暴露解空间结构、自由变量和不相容性，而不是单独替代求解过程。

• 浮点噪声清理: 非常小的值可能以浮点噪声形式残留，应该使用容差判断，不要期待绝对零。

坑点:

• 主元归一化: 如果主元接近0，除法会放大噪声。这也是为什么最好从带主元的REF开始。

• 精确代数 vs 浮点数: 不要把符号推导中的“严格为零”直接套到数值结果上。

参数:

void

返回值:

Mat - RREF形式的矩阵

高斯-约旦消元法求逆¶

Mat Mat::inverse_gje() const;

描述:

使用高斯-约旦消元法计算方阵的逆矩阵。

数学原理:

该方法构造增广矩阵 [A | I]，然后做行变换直到左半边变成 I。此时右半边就是 A⁻¹。这个过程直观，但数值稳定性通常不是最好的。

使用建议:

• 适合教学和小矩阵: 它能清晰展示逆矩阵是如何构造出来的。

• 不适合重复求解: 如果只是要解 Ax = b，直接求解方程组通常比显式求逆更稳、更快。

• 要关注条件数: 矩阵可逆不代表数值上安全；病态矩阵会让逆矩阵非常不可靠。

坑点:

• 必须是方阵: 非方阵应视为无效输入。

• 奇异或近奇异矩阵: 严格奇异时会失败；近奇异时也可能得到非常不稳定的结果。

• 优先考虑分解法: 多右端项场景下，LU/Cholesky 通常是更好的工程选择。

参数:

void

返回值:

Mat - 如果矩阵可逆则返回逆矩阵，否则返回空矩阵

点积¶

float Mat::dotprod(const Mat &A, const Mat &B);

描述:

计算两个向量（Nx1）的点积。

参数:

• const Mat &A: 输入向量 A (Nx1)

• const Mat &B: 输入向量 B (Nx1)

返回值:

float - 计算得到的点积值

解线性方程组¶

Mat Mat::solve(const Mat &A, const Mat &b) const;

描述:

使用高斯消元法和回代求解线性系统Ax = b。这是适用于良条件系统的直接方法。

数学原理:

该方法包括两个阶段：

1. 前向消元: 将增广矩阵[A|b]转换为上三角形式

2. 回代: 从下到上求解Ux = y

算法:

• 创建增广矩阵 [A | b]

• 应用高斯消元得到 [U | y]，其中U是上三角矩阵

• 使用回代求解Ux = y：xᵢ = (yᵢ - Σⱼ₌ᵢ₊₁ⁿ Uᵢⱼxⱼ) / Uᵢᵢ

参数:

• const Mat &A : 系数矩阵 (N×N)，必须是方阵且非奇异。

• const Mat &b : 右端项向量 (N×1)。

返回值:

Mat - 解向量 (N×1)，包含方程Ax = b的根。如果系统是奇异的或不兼容的，返回空矩阵。

使用建议:

• 单个系统: 对于求解一个系统很高效。对于具有相同A的多个系统，使用LU分解 + solve_lu()。

• 条件数: 对于病态矩阵，性能会下降。如果结果不准确，检查条件数。

• 奇异系统: 如果A是奇异的（det(A) = 0），返回空矩阵。对于秩亏系统，使用SVD + 伪逆。

• 性能: 消元为O(n³)，回代为O(n²)。总计O(n³)。

• 替代方法:

• 对于SPD矩阵：使用Cholesky分解 + solve_cholesky()（更快）

• 对于多个右端项：使用LU分解 + solve_lu()（更高效）

• 对于超定系统：使用QR分解 + solve_qr()（最小二乘）

坑点:

• 维度必须匹配: A 必须是方阵，b 的行数必须一致。维度错了通常不是“输入小问题”，而是建模本身有问题。

• 隐含不相容性: 增广矩阵里如果存在矛盾行，消元后才会暴露出来。不要默认每个线性系统都有解。

• 不要先求逆: 显式求 A⁻¹ 再乘 b 通常更慢，也更不稳定。

带状矩阵求解¶

Mat Mat::band_solve(Mat A, Mat b, int k) const;

描述:

使用优化的高斯消元法求解带状矩阵方程组 Ax = b。

数学原理:

这仍然是高斯消元，只是把消元限制在对角线附近的非零带内。矩阵确实是带状的时候，算法就能少碰大量无意义的零元素，节省时间和内存访问。

使用建议:

• 带状结构必须真实存在: 如果矩阵只是“看起来有点带状”，优化收益就会很有限。

• 带宽决定收益: k 越小，收益越大；如果 k 大到覆盖了大半个矩阵，就接近普通消元了。

• 适合结构化问题: 常见于离散微分方程、样条系统、局部耦合模型。

坑点:

• 带宽填小会出错: 低估 k 可能会把重要耦合忽略掉，得到错误答案。

• 仍然可能奇异: 带状矩阵也一样可能奇异或病态。

• 不是稳定性魔法: 这是对消元的结构优化，不是数值稳定性的万能修复。

参数:

• Mat A: 系数矩阵 (NxN) - 带状矩阵

• Mat b: 结果向量 (Nx1)

• int k: 矩阵的带宽（非零带的宽度）

返回值:

Mat - 解向量 (Nx1)，包含方程 Ax = b 的根

线性系统求根¶

Mat Mat::roots(Mat A, Mat y) const;

描述:

使用不同方法求解矩阵。这是 'solve' 函数的另一种实现，原理上没有区别。此方法使用高斯消元法求解线性系统 A * x = y。

说明:

roots() 本质上是 solve() 的语义别名，便于在“求根”这种表达里直接使用。它的代数核心和数值行为与 solve() 一致。

坑点:

• 不会比 solve 更稳: 如果 solve() 对某个输入不稳定，roots() 也不会神奇地变好。

• 预条件相同: A 仍然必须是方阵且非奇异，y 的行数也必须匹配。

参数:

• Mat A: 矩阵 [N]x[N]，包含输入系数

• Mat y: 向量 [N]x[1]，包含结果值

返回值:

Mat - 矩阵 [N]x[1]，包含根