【python】【原创】常用线性代数矩阵运算函数的函数库
本帖最后由 白冥 于 2025-1-23 22:08 编辑本贴将详细介绍一组用于矩阵计算和线性代数运算的Python函数。函数库的函数都是自己尝试实现的，包含了一些经典的线性代数算法，包括QR分解、特征值求解、线性方程组求解、行列式计算、矩阵求逆、基变换、矩阵变换和矩阵收缩。所有函数均采用高效的算法，并且通过模块化设计，每个算法的实现和使用都尽可能简洁明了。
static/image/hrline/4.gif
### 目录
1. **QR 分解** (`QR` 函数)
2. **特征值和特征向量计算** (`eig` 函数)
3. **线性方程组求解** (`solve_linear_system` 函数)
4. **行列式计算** (`det` 函数)
5. **矩阵求逆** (`inverse_matrix` 函数)
6. **基变换** (`basis_transform` 函数)
7. **矩阵变换** (`transform` 函数)
8. **矩阵收缩** (`contract_matrix` 函数)
9. **代码实现和数学背景**
---
static/image/hrline/5.gif
### 1. **QR 分解** (QR函数)
QR分解是将一个矩阵 ( A ) 分解为正交矩阵 ( Q ) 和上三角矩阵 ( R ) 的过程。此函数实现了经典的Gram-Schmidt正交化方法。
def QR(A):
# 内部函数：转置矩阵
def trans(matrix):
return list(zip(*matrix))
# 内部函数：计算向量的模
def mod(v):
return sum() ** 0.5
# 内部函数：归一化向量
def unit(v):
return
# 内部函数：向量点积
def dot(v_0, v_1):
return sum()
# 内部函数：计算向量的投影
def proj(v_0, v_1):
return [(x * dot(v_0, v_1)) / mod(v_1) ** 2 for x in v_1]
# 内部函数：正交化
def orthogonalize(v, Q_T):
return )]
A_T = trans(A)# 获取矩阵A的转置
if any(mod(v) == 0 for v in A_T):# 如果A的某一列全为零，则返回None
(放过来只是为了吐槽本帖用不了用不了markdown)
(然后吐槽本帖上面不能用样式，排版还经常会乱)
_(:з」∠)__(:з」∠)_手动分割线_(:з」∠)__(:з」∠)_
QR分解是一种将一个矩阵分解为两个矩阵的方法，通常用于线性代数中的数值计算。它将一个给定的矩阵 A 分解为两个矩阵的乘积：A = QR
其中， Q 是一个正交矩阵，即 Q^T Q = I，其中I是单位矩阵），而R 是一个上三角矩阵。
QR分解的主要应用包括：
1. **求解线性方程组**：通过分解矩阵，可以更有效地求解线性方程组。
2. **最小二乘问题**：当需要拟合数据时，QR分解可以用来解最小二乘问题。
3. **特征值计算**：QR分解是计算矩阵特征值的一种重要方法。
QR分解的常见方法包括**格拉姆-施密特过程**、**豪斯霍尔德变换**和**吉文斯旋转**等。
这里用的Gram-Schmidt就是工程线代课程里教的格拉姆-施密特正交化。
**函数解释：**
- **输入：** 矩阵 \( A \)（任意形状的二维数组）
- **输出：** 矩阵 \( Q \) 和矩阵 \( R \)，其中 \( A = QR \)
- **功能：** 该函数使用Gram-Schmidt方法将矩阵 \( A \) 分解为正交矩阵 \( Q \) 和上三角矩阵 \( R \)。正交矩阵的列是单位向量，且互相正交。上三角矩阵 \( R \) 存储了矩阵 \( A \) 中列之间的线性关系。
static/image/hrline/4.gif
### 2. **特征值和特征向量计算** (`eig` 函数)
def eig(A, max_iter, tol):
def trans(matrix):
return list(zip(*matrix))
def dot(v_0, v_1):
return sum()
def multiplication(M, N):
return [ for row in M]
def converges(next_A):
squares = sum( ** 2 for i in range(len(A)) for j in range(len(A)) if i != j])
return squares < tol
next_A = trans(trans(A))# 初始化矩阵A
for _ in range(max_iter):
Q, R = QR(next_A)# QR分解
next_A = multiplication(R, Q)# 计算新的矩阵A
if converges(next_A):# 判断是否收敛
break
eigenvalue = for i in range(len(A))]# 提取对角线上的特征值
return eigenvalue, trans(Q)# 返回特征值和特征向量
_(:з」∠)__(:з」∠)_手动分割线_(:з」∠)__(:з」∠)_
QR算法是一种用于计算矩阵特征值的数值算法。它通过迭代方式将一个方阵逐渐转化为一种更容易分析的形式，通常是上三角矩阵，从而使得其特征值容易计算。
### QR算法的基本步骤：
1. **分解**：给定一个方阵 \( A \)，首先进行QR分解，即将矩阵 \( A \) 分解为 \( A = QR \)，其中 \( Q \) 是正交矩阵，\( R \) 是上三角矩阵。
2. **重构**：然后构造新的矩阵 \( A' = RQ \)。
3. **迭代**：重复执行QR分解和重构的过程，得到一系列矩阵 \( A^{(k)} \)，每一次的分解结果都接近于一个上三角矩阵。
随着迭代的进行，矩阵 \( A^{(k)} \) 会逐渐接近一个上三角矩阵，其对角线上的元素就是原矩阵 \( A \) 的特征值。
### QR算法的特点：
- **收敛性**：QR算法通常是收敛的，尤其是在没有特殊结构的矩阵下。它的收敛速度较快，通常几次迭代后就能得到良好的近似特征值。
- **稳定性**：QR算法在数值计算中较为稳定，适用于大多数实对称矩阵和一般矩阵的特征值计算。
### 应用：
QR算法广泛用于特征值问题的求解，尤其是在需要计算大规模矩阵的特征值时。例如，它被用于谱分解、主成分分析（PCA）、控制系统分析等领域。
说这些是因为大家都知道特征值和特征向量，但是不知道QR算法。QR算法这东西用一句话说就是，近似。就像五次方程没有根式解一样，因为求特征值本质上就是要解多项式，QR算法近似特征值算是很精准的了。
static/image/hrline/4.gif
### 3. **线性方程组求解**
def solve_linear_system(A, b):
def trans(matrix):
return list(zip(*matrix))
A_T = trans(A)# 转置矩阵A
solution = * len(A_T)# 初始化解向量
Q, R = QR(A)# 进行QR分解
if not Q:
return None
b = ) for row, x in zip(trans(Q), b)]# 转换b向量
for i in range(len(A) - 1, -1, -1):
terms = []
for j in range(i + 1, len(A)):
terms.append(R * solution)
solution = (b - sum(terms)) / R# 使用回代法解线性方程
return solution# 返回解向量
希望论坛不要吞代码（
**函数解释：**
- **输入：** 矩阵 \( A \) 和向量 \( b \)（线性方程组 \( A \times x = b \)）
- **输出：** 向量 \( x \)，即线性方程组的解
- **功能：** 通过QR分解将线性方程组转化为上三角形式，并使用回代法求解。首先对矩阵 \( A \) 进行QR分解，接着利用正交矩阵 \( Q \) 将右侧向量 \( b \) 转换，然后通过回代法求解得到解向量。
这个其实没什么好说的了（
解方程嘛（
又不是高次的（
其实本来想写一个处理最小二乘问题的函数，不过这样就显得QR分解疑似有点太常用了（
不过雀食常用（
有人问为什么解线性方程为什么提最小二乘，这也没有啊？
最小二乘问题通常的形式是，给定一个矩阵 A 和一个向量 b，我们希望找到一个向量 x，使得 Ax ≈ b 的误差最小。换句话说，我们要最小化目标函数：||Ax ≈ b||²
在这种情境下，QR分解被用来将矩阵 A分解为一个正交矩阵 Q和一个上三角矩阵 R，即：A = QR
通过将最小二乘问题代入QR分解形式，可以将原问题转化为一个更容易求解的形式。具体来说，最小化 ||Ax - b||² 就变成了最小化 ||QRx - b||²，由于 Q是正交矩阵，我们可以通过左乘 Q^T来简化为：||Rx - Q^Tb||²
接着，由于 R是上三角矩阵，问题变得易于通过回代法直接求解。
裆燃啦，这里就没有写这个函数，只是解普通的线性方程组。∠( ᐛ 」∠)_
static/image/hrline/4.gif
### 4. **行列式计算** (`det` 函数)
def det(A):
_, R = QR(A)# 获取QR分解
return math.prod( for i in range(len(R))])# 计算行列式
嘿咻咻，这里就不写解释了，应该没有人不知道行列式的值是上三角矩阵主对角线元素之积吧（莫名狗头
_(:з」∠)__(:з」∠)_哎，就是分割_(:з」∠)__(:з」∠)_
哎，分割个寂寞~
static/image/hrline/4.gif
### 5. **矩阵求逆** (`inverse_matrix` 函数)
def inverse_matrix(A):
def multiply(num, M):
return [ for r in M]
def cof(M, i, j):
return ((-1) ** (i + j + 2)) * det((M[:i] + M)[:i] + (M[:i] + M))
def adj(M):
return [ for i in range(len(M))]
d = det(A)
if d == 0:
return None
adj_A = adj(A)# 计算伴随矩阵
return multiply(1 / d, adj_A)# 通过伴随矩阵计算逆矩阵
**函数解释：**
- **输入：** 矩阵 \( A \)
- **输出：** 矩阵 \( A \) 的逆矩阵
- **功能：** 通过伴随矩阵法计算矩阵的逆。如果矩阵的行列式为零，则返回 `None`，表示矩阵不可逆。
嗯，是的，不用增广矩阵求，主要是太麻烦了∠( ᐛ 」∠)_
static/image/hrline/4.gif
### 6. **基变换** (`basis_transform` 函数)
（这个其实没什么好说的）
def basis_transform(bvgroup_0, bvgroup_1):
def trans(M):
return list(zip(*M))
A = trans(bvgroup_0)# 转置矩阵
B = trans(bvgroup_1)# 转置矩阵
B_r = inverse_matrix(B)# 计算B的逆
return multiplication(B_r, A)# 返回基变换矩阵
看吧～
怎么说呢？
保险起见，把7安排上
### 7. **矩阵变换** (`transform` 函数)
def transform(A, bvgroup_0, bvgroup_1):
def trans(M):
return list(zip(*M))
P = basis_transform(bvgroup_0, bvgroup_1)# 获取基变换矩阵
P_r = inverse_matrix(P)# 获取基变换矩阵的逆
return multiplication(multiplication(P_r, A), P)# 执行矩阵变换好了，其实放代码你不一定看得懂，数学原理：
若A是一个n阶方阵，y₁是一个n维向量，那么A对y的映射结果为Ay₁，设这个映射结果为y'₁，y'₁=Ay₁
这个结果默认是以(α₁,α₂,…,αₙ)作为线性空间的一组基
由于向量在线性空间中是客观存在的，则以任何基的同一线性空间中，都有：
y₀=(α₁,α₂,…,αₙ)y₁,y'₀=(α₁,α₂,…,αₙ)y'₁
把y'₀用已知量y₁表示：
y'₀=(α₁,α₂,…,αₙ)y'₁=(α₁,α₂,…,αₙ)Ay₁
假设以{β₁,β₂,…,βₙ}作为线性空间的一组基，也有y₀=(β₁,β₂,…,βₙ)y₂，已知对于的从(β₁,β₂,…,βₙ)到(α₁,α₂,…,αₙ)的基变换矩阵K=(β₁,β₂,…,βₙ)⁻¹(α₁,α₂,…,αₙ)
由于以{β₁,β₂,…,βₙ}作为线性空间的一组基时，y₀在线性空间中的坐标不再是y₁，因此对应的变换矩阵不再是A，设这个未知变换矩阵为B，即y'₀=(β₁,β₂,…,βₙ)By₂，有y₂=Ky₁
通过基变换，有(α₁,α₂,…,αₙ)=(β₁,β₂,…,βₙ)K，(β₁,β₂,…,βₙ)=(α₁,α₂,…,αₙ)K⁻¹
y'₀=(β₁,β₂,…,βₙ)By₂=(α₁,α₂,…,αₙ)K⁻¹BKy₁，又y'₀=(α₁,α₂,…,αₙ)Ay₁，则A=K⁻¹BK
再经过运算，有：
A=K⁻¹BK
KA=(KK⁻¹)BK=BK
KAK⁻¹=B
综上所述，若已知在以(α₁,α₂,…,αₙ)为基时，变换矩阵为A，要求同一线性空间下，以(β₁,β₂,…,βₙ)为基时，新的未知变换矩阵B，引入基变换矩阵K，有：
B=KAK⁻¹，K=(β₁,β₂,…,βₙ)⁻¹(α₁,α₂,…,αₙ)
哎！ 啊啊啊，我瞎了，向量什么的，矩阵什么的，相乘什么的，分解什么的，放过我吧，本可已经被你们轰到外翻了 靴靴有点像在看天书惹
知识轰炸了我的大脑
炸翻天惹
https://img.gamemale.com/album/202412/11/164541bmmaavuvijmxif8v.jpg.thumb.jpg 我天，我居然在GM论坛学PYTHON！！！！！ 内容全面，原理清晰，实用性强。楼楼真厉害
{:6_188:} 本坛就是需要更多楼主这样的人才 以前在学校的时候我们专业只需要学一部分线性代数，考试还得手算，学得头疼 以前在项目写算法都是从网上直接抄一个，不研究具体的写法了 天哪，本来还想着要不要重新学一下高中数学拓展业务的，看了两眼楼主帖子，决定去复习物理了。。 这个也可以吗？感觉我也可以把我上课的课件上传，保证原创 惹，看到线性代数就想到我大学时候的受苦受难惹... 都忘掉了谢谢，线代学了跟没学一样:'(全都还给老师惹 所以我为啥不直接用包装好的函数而是学底层？ 看起来好专业惹，梦回被高数支配的日子噜(˃ ⌑ ˂ഃ ) https://img.gamemale.com/album/202408/03/102121ubl2b4leibsmzs7s.gif感觉不太想研究算法原理，要用到就用呢，调用一下matlab的库 感觉又回忆起在电脑面前不停敲代码，不停调错的暴躁日子了 啊？不是，虽然博主很厉害。但我是真的想不到能在这看到如此高知识含量的帖子{:4_87:} 太强了惹，竟然还有第二个python数学帖{:6_172:} 强烈要求线性代数考试的时候可以带电脑 {:4_114:} 又来学习了U•ェ•*U应该先这样然后再这样然后结论就出来（胡言乱语ing）
大佬太厉害了，继续膜拜大佬٩(๑•̀ω•́๑)۶
页:
[1]
2