TheAlgorithms/Python 线性代数模块解析：从纯 Python 的 Vector/Matrix 库到 NumPy 数值算法实现

本篇以 linear_algebra/README.md 为核心，系统讲解 TheAlgorithms/Python 仓库中的线性代数模块：先吃透 src/lib.py 中纯 Python 实现的 Vector / Matrix 类与 axpy、unit_basis_vector 等工具函数的用法和源码细节，再结合仓库中基于 NumPy 的高斯消元、Jacobi 迭代、LU 分解、共轭梯度与幂迭代等算法文件，形成"基础数据结构 + 数值算法"两条完整学习路径。

模块定位与目录结构

该模块在仓库中的入口说明位于 linear_algebra/README.md，其核心定位是"包含用于执行线性代数的类和函数的库"。从目录结构看，模块分为两层：

• 库层：linear_algebra/src/lib.py 提供不依赖第三方库的 Vector、Matrix 类与若干函数，配套单元测试 linear_algebra/src/test_linear_algebra.py；
• 算法层：模块根目录与 src 子目录下各独立脚本实现具体的数值线性代数算法，例如 gaussian_elimination.py、jacobi_iteration_method.py、lu_decomposition.py、matrix_inversion.py、conjugate_gradient、power_iteration 等。

Vector 类：任意维向量的完整操作集

README 对 Vector 的概括是"表示任意大小的向量及相关操作的类"。其完整方法清单（继承自原文档并对照源码补全）如下：

成员	说明
constructor(components)	用分量列表初始化向量，components 缺省时得到空向量
set(components) / change_component(pos, value)	替换向量的分量（change_component 支持负下标，源码中使用 assert 做前置条件检查）
__str__()	字符串表示，如 (1,2,3)
component(i)	获取第 i 个分量（0 起），支持负索引
__len__()	向量规模（分量个数），使 len(v) 可用
euclidean_length()	欧几里得长度，空向量会抛出 Exception: Vector is empty
运算符 + / -	向量加/减，尺寸不一致时抛出异常
运算符 *	双重语义：与标量相乘得向量，与向量相乘得点积（内积）
copy()	拷贝并返回新向量

源码中的几个值得注意的实现细节：

1. * 运算符的双重重载。lib.py 的 __mul__ 方法 使用 @overload 标注：Vector * float 返回 Vector，Vector * Vector 返回 float（点积）。点积实现为逐分量相乘后 sum，这也是 angle 方法（两向量夹角，支持弧度/角度两种输出）的基础。
2. 内置 docstring 与 doctest。euclidean_length 的 docstring 直接内嵌可运行的测试用例（见 lib.py 第 155-173 行）：Vector([2, 3, 4]).euclidean_length() 返回 5.385164807134504，空向量抛出异常。这是 README "Documentation" 一节提到的用 Python 内建 help(...) 查看文档的基础。
3. README 未列出但源码中存在的方法：从源码结构看，Vector 还实现了 __eq__（逐分量相等判断，见 lib.py 第 98-106 行）与 angle(other, deg=False)（夹角，见 lib.py 第 175-193 行），实际 API 比 README 的概览更完整。

最小使用示例（与测试用例 test_linear_algebra.py 中的断言一致）：

from src.lib import Vector, unit_basis_vector, axpy

x = Vector([1, 2, 3])
print(len(x), x.component(0), x.component(-1))   # 3 1 3
print(x + unit_basis_vector(3, 2))               # (1,2,4)
print(x * 2, x * Vector([1, 1, 1]))              # (2,4,6) 6
print(axpy(2, x, unit_basis_vector(3, 0)))      # (3,4,6)  即 2x + y

模块级工具函数

README 列出的四个函数在源码 lib.py 中均有对应实现，且都通过 assert 做前置条件校验：

• zero_vector(dimension)（L196-L202）：返回指定维度的零向量，是构造 Matrix 矩阵向量乘结果时的起点；
• unit_basis_vector(dimension, pos)（L205-L215）：返回第 pos 位为 1 的基向量（0 起索引），是构造测试向量的高频工具；
• axpy(scalar, x, y)（L218-L228）：实现数值计算中最基础的 BLAS-1 运算 y + scalar*x，实现仅一行 return x * scalar + y，直观展示了运算符重载如何拼装成算法原语；
• random_vector(n, a, b)（L231-L240）：生成分量在 [a, b] 区间内的随机整数向量（调用 random.seed(None) 重新播种），常用于算法的随机化测试输入。

Matrix 类：从基本运算到行列式

Matrix 类（lib.py 第 243 行起）表示任意 W×H 矩阵，README 概括的方法与源码实际 API 对照如下：

README 方法	源码实现要点
__str__()	每行以 | 包围、逗号分隔输出
运算符 + / -	逐元素加减，尺寸不匹配抛 Exception("matrix must have the same dimension!")
运算符 *	双重重载：Matrix * Vector 执行矩阵向量乘（要求向量长度等于矩阵宽度），Matrix * 标量 返回放大后的新矩阵，见 lib.py 第 319-351 行
component(x, y) / change_component(x, y, value)	越界时抛 Exception("change_component: indices out of bounds")
width() / height()	返回列数/行数
determinant()	仅对方阵有效，采用拉普拉斯展开（沿第一行）递归求解

源码中超出 README 概览的部分（"从源码结构看"）：Matrix 额外提供 minor(x, y)（余子式）与 cofactor(x, y)（代数余子式，(-1)^(x+y) 乘余子式），二者正是 determinant 递归的支撑——determinant 实现 对 1×1、2×2 有直接公式，更大的方阵按第一行展开为 matrix[0][y] * cofactor(0, y) 之和。需要说明的是，拉普拉斯展开的时间复杂度较高，适合理解原理和小规模验证，不适合大规模数值计算。

与之配套的工厂函数：

• square_zero_matrix(n)（L427-L432）：返回 N×N 零矩阵；
• random_matrix(width, height, a, b)（L435-L444）：返回 W×H、分量为 [a, b] 内随机整数的矩阵。

矩阵向量乘的典型用法（对应 README "operator *" 条目）：

from src.lib import Matrix, Vector, random_matrix

m = Matrix([[1, 2], [3, 4]], 2, 2)
v = Vector([1, 1])
print(m * v)          # (3,7)
print(m * 2)          # |2,4|
                      # |6,8|
print(Matrix.random if False else random_matrix(2, 2, 1, 9).width())  # 2

使用与文档：import、help 与单元测试

README 的 Usage 一节给出了两种使用方式：将 src 目录下的 lib.py 导入项目，或直接使用编译后的字节码文件 lib.pyc。在当前仓库中，推荐的验证路径是：

# 以模块方式运行 lib 的配套单元测试
python3 -m unittest -v

README "Tests" 一节明确指出测试位于 src/tests.py（当前仓库中对应文件为 linear_algebra/src/test_linear_algebra.py，同时使用 unittest 与 pytest 框架，覆盖 component、__str__、__len__、euclidean_length 等行为）；而 lib.py 内部各方法的 docstring 本身也是 doctest 用例，可用 python3 -m doctest linear_algebra/src/lib.py -v 方式验证。文档方面，README 建议直接尝试 help(Vector)、help(unit_basis_vector) 以及 help(CLASSNAME.METHODNAME)。

算法层：基于 NumPy 的求解与分解实现

模块根目录的独立脚本则面向"用 NumPy 解决具体数值问题"，每个文件自带可运行的 doctest：

高斯消元（gaussian_elimination.py）

gaussian_elimination(coefficients, vector) 将系数矩阵与常数列拼成增广矩阵，前向消元化为上三角，再调用 retroactive_resolution 回代求解。docstring 中给出的完整例子：

>>> gaussian_elimination([[1, -4, -2], [5, 2, -2], [1, -1, 0]], [[-2], [-3], [4]])
array([[ 2.3 ],
       [-1.7 ],
       [ 5.55]])

注意：该实现未做主元选择（pivoting），非方阵输入会返回空数组。带主元选取的改进版本见 src/gaussian_elimination_pivoting.py 中的 solve_linear_system：它先按列最大值交换行，并对主元绝对值小于 1e-8 的情形抛出 ValueError: Matrix is singular，docstring 示例的 4×4 系统可直接运行验证。

Jacobi 迭代（jacobi_iteration_method.py）

针对严格对角占优的线性方程组的迭代求解。jacobi_iteration_method(coefficient_matrix, constant_matrix, init_val, iterations) 在进入迭代前会依次校验：系数矩阵必须为 n×n、常数列必须为 n×1、初值个数与行数一致、迭代次数至少为 1，并通过辅助函数 strictly_diagonally_dominant 检查对角占优性（不满足则抛 ValueError）。其 doctest 示例：

>>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
>>> constant = np.array([[2], [-6], [-4]])
>>> jacobi_iteration_method(coefficient, constant, [0.5, -0.5, -0.5], 3)
[0.909375, -1.14375, -0.7484375]

实现上用 NumPy 布尔掩码 ~np.eye(...) 一次性剔除对角元，迭代更新公式 new_val = (sum_product_rows + val_last) / denominator 向量化完成，是"教科书循环 + NumPy 向量化"的典型对照样本。

LU 分解（lu_decomposition.py）

lower_upper_decomposition(table) 在不选主元的前提下将方阵分解为下三角 L（对角线为 1）与上三角 U。文件头部 docstring 完整陈述了存在性条件：可逆矩阵当且仅当所有顺序主子式非零时存在 LU 分解；奇异矩阵秩为 k 时，前 k 个顺序主子式非零即可。doctest 覆盖了三类典型情形，包括 [[0, 1], [1, 0]]（可逆但第一个顺序主子式为 0，抛 ArithmeticError: No LU decomposition exists）与 [[1, 0], [1, 0]]（奇异但分解存在）。

矩阵求逆（matrix_inversion.py）

invert_matrix(matrix) 是对 np.linalg.inv 的薄封装：捕获 np.linalg.LinAlgError 并转为语义更清晰的 ValueError: Matrix is not invertible。doctest 示例 invert_matrix([[4.0, 7.0], [2.0, 6.0]]) 返回 [[0.6, -0.7], [-0.2, 0.4]]。

src/ 目录下的进阶算法

src 子目录除 lib.py 外，还收录了一批面向特征值与矩阵结构的算法实现，均带 doctest 自验证：

• 共轭梯度法 src/conjugate_gradient.py：求解对称正定（SPD）矩阵系统 Ax = b。conjugate_gradient(spd_matrix, load_vector, max_iterations=1000, tol=1e-8) 内部先用 np.linalg.eigh 的特征值符号判定 SPD 性，再以 alpha = r·r / (p·Ap) 更新解向量，收敛判据取残差变化与解变化的较大者。配套函数 test_conjugate_gradient 用随机生成的 SPD 系统将其结果与 np.linalg.solve 对照（误差阈值 1e-6），是验证迭代算法正确性的实用范式；
• 幂迭代法 src/power_iteration.py：求矩阵最大特征值及对应特征向量。power_iteration(input_matrix, vector, error_tol=1e-12, max_iterations=100) 每轮做 w = A @ v 并归一化，用瑞利商估计特征值，以相对变化量判断收敛；docstring 明确约束输入为实矩阵或酉（Hermitian）矩阵，test_power_iteration 同时验证了实数与复数（构造 A = A + i*triu(A,1) - i*triu(A,1).T 的酉矩阵）两种情形并与 np.linalg.eigh 结果比对；
• 矩阵秩 src/rank_of_matrix.py：纯列表输入求秩，doctest 中 [[1,2,3],[4,5,6],[7,8,9]] 的秩为 2；
• 瑞利商 src/rayleigh_quotient.py：含 is_hermitian 判定与瑞利商计算，服务于特征值估计；
• 舒尔补 src/schur_complement.py：对 2×2 块矩阵 [A B; C D] 计算 Schur 补，A 奇异时可通过 pseudo_inv 参数提供伪逆；
• 点拟合多项式 src/polynom_for_points.py：points_to_polynomial(coordinates) 给定若干 (x, y) 点求拟合多项式，空输入抛 ValueError；
• 2D 变换矩阵 src/transformations_2d.py：提供 scaling、rotation、projection、reflection 四个函数，文件头部 docstring 直接给出 45° 旋转/投影矩阵的数值结果，是线性代数与图形学交叉的轻量示例。

小结

linear_algebra 模块的价值在于两条互补路径：src/lib.py 以约 450 行纯 Python 代码讲透了向量/矩阵数据结构、运算符重载、余子式与拉普拉斯展开等"原理层"内容，并通过 docstring 内嵌的 doctest 与 test_linear_algebra.py 保证可验证性；而根目录与 src 中的 NumPy 脚本则覆盖高斯消元、选主元消元、Jacobi 迭代、LU 分解、共轭梯度、幂迭代等"算法层"实现。读者可先按 README 的用法导入 lib.py 建立直觉，再逐个运行各算法文件的 doctest（每个文件尾部均以 if __name__ == "__main__": doctest.testmod() 自测），以最小成本完成从数据结构到数值算法的完整练习闭环。