|
本帖最后由 白冥 于 2025-1-23 22:08 编辑
本贴将详细介绍一组用于矩阵计算和线性代数运算的Python函数。函数库的函数都是自己尝试实现的,包含了一些经典的线性代数算法,包括QR分解、特征值求解、线性方程组求解、行列式计算、矩阵求逆、基变换、矩阵变换和矩阵收缩。所有函数均采用高效的算法,并且通过模块化设计,每个算法的实现和使用都尽可能简洁明了。

### 目录
1. **QR 分解** (`QR` 函数)
2. **特征值和特征向量计算** (`eig` 函数)
3. **线性方程组求解** (`solve_linear_system` 函数)
4. **行列式计算** (`det` 函数)
5. **矩阵求逆** (`inverse_matrix` 函数)
6. **基变换** (`basis_transform` 函数)
7. **矩阵变换** (`transform` 函数)
8. **矩阵收缩** (`contract_matrix` 函数)
9. **代码实现和数学背景**
---

### 1. **QR 分解** (QR函数)
QR分解是将一个矩阵 ( A ) 分解为正交矩阵 ( Q ) 和上三角矩阵 ( R ) 的过程。此函数实现了经典的Gram-Schmidt正交化方法。
- def QR(A):
- # 内部函数:转置矩阵
- def trans(matrix):
- return list(zip(*matrix))
-
- # 内部函数:计算向量的模
- def mod(v):
- return sum([x ** 2 for x in v]) ** 0.5
-
- # 内部函数:归一化向量
- def unit(v):
- return [x / mod(v) for x in v]
-
- # 内部函数:向量点积
- def dot(v_0, v_1):
- return sum([x * y for x, y in zip(v_0, v_1)])
-
- # 内部函数:计算向量的投影
- 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 [x - y for x, y in zip(v, *[proj(v, p) for p in Q_T])]
-
- 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 \) 中列之间的线性关系。

### 2. **特征值和特征向量计算** (`eig` 函数)
- def eig(A, max_iter, tol):
- def trans(matrix):
- return list(zip(*matrix))
-
- def dot(v_0, v_1):
- return sum([x * y for x, y in zip(v_0, v_1)])
-
- def multiplication(M, N):
- return [[dot(row, col) for col in trans(N)] for row in M]
-
- def converges(next_A):
- squares = sum([next_A[i][j] ** 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 = [next_A[i][i] 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算法近似特征值算是很精准的了。

### 3. **线性方程组求解**
- def solve_linear_system(A, b):
- def trans(matrix):
- return list(zip(*matrix))
-
- A_T = trans(A) # 转置矩阵A
- solution = [0] * len(A_T) # 初始化解向量
-
- Q, R = QR(A) # 进行QR分解
- if not Q:
- return None
-
- b = [sum([x * y for y in row]) 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[i][j] * solution[j])
- solution[i] = (b[i] - sum(terms)) / R[i][i] # 使用回代法解线性方程
- 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是上三角矩阵,问题变得易于通过回代法直接求解。
裆燃啦,这里就没有写这个函数,只是解普通的线性方程组。∠( ᐛ 」∠)_

### 4. **行列式计算** (`det` 函数)
- def det(A):
- _, R = QR(A) # 获取QR分解
- return math.prod([R[i][i] for i in range(len(R))]) # 计算行列式
复制代码
嘿咻咻,这里就不写解释了,应该没有人不知道行列式的值是上三角矩阵主对角线元素之积吧(莫名狗头
_(:з」∠)__(:з」∠)_哎,就是分割_(:з」∠)__(:з」∠)_
哎,分割个寂寞~

### 5. **矩阵求逆** (`inverse_matrix` 函数)
- def inverse_matrix(A):
- def multiply(num, M):
- return [[x * num for x in r] for r in M]
-
- def cof(M, i, j):
- return ((-1) ** (i + j + 2)) * det((M[:i] + M[i + 1:])[:i] + (M[:i] + M[i + 1:])[j + 1:])
-
- def adj(M):
- return [[cof(M, i, j) for j in range(len(M))] 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`,表示矩阵不可逆。
嗯,是的,不用增广矩阵求,主要是太麻烦了∠( ᐛ 」∠)_

### 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=(β₁,β₂,…,βₙ)⁻¹(α₁,α₂,…,αₙ)
哎! |
本帖子中包含更多资源
您需要 登录 才可以下载或查看,没有账号?立即注册
x
评分
-
查看全部评分
白冥
:啊!怎么帖子还有字数限制
白冥
:真吞字啦
白冥
:源代码:
import math
def QR(A):
def trans(matrix):
return list(zip(*matrix))
def mod(v):
return sum([x ** 2 for x in v]) ** 0.5
def unit(v):
return [x / mod(v) for x in v]
def dot(v_0,v_1):
return sum([x*y for x,y in zip(v_0,v_1)])
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 [x - y for x, y in zip(v, *[proj(v, p) for p in Q_T])]
A_T=trans(A)
if any(mod(v)==0 for v in A_T):
return None,None
Q_T = [unit(orthogonalize(v, Q_T)) for v in A_T]
Q = trans(Q_T)
R = [[dot(A_T[k], Q_T) for k in range(i, len(A))] for i in range(len(A))]
return Q, R
def eig(A,max_iter,tol):
def trans(matrix):
return list(zip(*matrix))
def dot(v_0,v_1):
return sum([x*y for x,y in zip(v_0,v_1)])
def multiplication(M,N):
return [[dot(row,col) for col in trans(N)] for row in M]
def converges(next_A):
squares=sum([next_A[j] ** 2 for i in range(len(A)) for j in range(len(A)) if i != j])
return squares < tol
next_A = trans(trans(A))
for _ in range(max_iter):
Q,R = QR(next_A)
next_A=multiplication(R,Q)
if converges(next_A):
break
eigenvalue=[next_A for i in range(len(A))]
return eigenvalue,trans(Q)
def solve_linear_system(A, b):
def trans(matrix):
return list(zip(*matrix))
A_T=trans(A)
solution=[0]*len(A_T)
Q, R=QR(A)
if not Q:
return None
b=[sum([x * y for y in row]) for row, x in zip(trans(Q), b)]
for i in range(len(A) - 1, -1, -1):
terms=[]
for j in range(i + 1, len(A)):
terms.append(R[j] * solution[j])
solution = (b - sum(terms)) / R
return solution
def det(A):
_, R = QR(A)
return math.prod([R for i in range(len(R))])
def inverse_matrix(A):
def multiply(num,M):
return [[x*num for x in r] for r in M]
def cof(M,i,j):
return ((-1) ** (i+j+2)) * det((M[:i] + M[i+1:])[:i] + (M[:i] + M[i+1:])[j+1:])
def adj(M):
return [[cof(M,i,j) for j in range(len(M))] for i in range(len(M))]
d=det(A)
if d==0:
return None
adj_A=adj(A)
return multiply(1/d,adj_A)
def basis_transform(bvgroup_0,bvgroup_1):
def trans(M):
return list(zip(*M))
def dot(v_0,v_1):
return sum([x*y for x,y in zip(v_0,v_1)])
def multiplication(M,N):
return [[dot(row,col) for col in trans(N)] for row in M]
A=trans(bvgroup_0)
B=trans(bvgroup_1)
B_r=inverse_matrix(B)
return multiplication(B_r,A)
def transform(A,bvgroup_0,bvgroup_1):
def trans(M):
return list(zip(*M))
def dot(v_0,v_1):
return sum([x*y for x,y in zip(v_0,v_1)])
def multiplication(M,N):
return [[dot(row,col) for col in trans(N)] for row in M]
P=basis_transform(bvgroup_0,bvgroup_1)
P_r=inverse_matrix(P)
return multiplication(multiplication(P_r,A),P)
def contract_matrix(A):
def trans(M):
return list(zip(*M))
eigvals,eigvecs=eig(A,100,1e-6)
P,_=QR(trans(eigvecs))
D=[[(eigvals if j==i else 0) for j in range(len(A))] for i in range(len(A))]
P_r=inverse_matrix(P)
return P,D,P_r
白冥
:求一求免费的追随
-
|