python奇异值分解法代码

以下是Python中使用奇异值分解法（SVD）的示例代码：

import numpy as np
from scipy.linalg import svd

# 生成一个矩阵
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

# 使用SVD进行分解
U, S, V = svd(A)

# 打印结果
print("U = ")
print(U)

print("S = ")
print(np.diag(S))

print("V = ")
print(V)

# 重构原始矩阵
reconstructed_A = np.dot(U, np.dot(np.diag(S), V))

# 打印重构后的矩阵
print("Reconstructed A = ")
print(reconstructed_A)

# 计算矩阵的秩
rank_A = len(np.where(S > 1e-10)[0])
print("Rank of A = ", rank_A)

相关问题

矩阵的奇异值分解证明

### 矩阵奇异值分解的证明

#### 定义与预备知识
奇异值分解(SVD)是一种重要的矩阵分解方法，在数值分析中有广泛应用。对于任意 \( m \times n \) 的实数矩阵 \( A \)，存在两个正交矩阵 \( U_{m\times m} \) 和 \( V_{n\times n} \)，使得：

\[ A = U\Sigma V^T \]

其中，\( \Sigma \) 是一个 \( m \times n \) 对角矩阵，其对角线上元素称为奇异值，并按降序排列。

#### 推导过程

考虑给定矩阵 \( A \in R^{m×n} \)[^1]，定义如下优化问题来寻找最大奇异值及其对应的左、右奇异向量:

\[
\begin{aligned}
& \max _{\|u\|=1,\|v\|=1}\left(u^{T}Av\right)^2 \\
& s.t.\quad u^Tu=1, v^Tv=1.
\end{aligned}
\]

通过拉格朗日乘子法可以得到最优解条件为：

\[ Av=\sigma_1u \]
\[ A^Tu=\sigma_1v \]

这里 \( \sigma_1=(u^TAv)\geqslant0 \) 即为最大的奇异值，而 \( (u,v) \) 则是最优单位长度的向量对[^3]。

进一步地，假设已经找到了前 k-1 个相互垂直的最大奇异值所对应的向量，则第 k 大奇异值可通过下面的方式获得：

\[
\begin{array}{l}
\max _{\substack{
\|u_k\|=1 \\
u_i^Tu_k=0,i<k\\
}}
\max _{\substack{
\|v_k\|=1 \\
v_j^Tv_k=0,j<k\\
}}(u_k^TAv_k)^2,
\end{array}
\]

重复此操作直到找到所有的非零奇异值为止。此时便完成了整个 SVD 过程。

最后注意到当且仅当输入矩阵是对称半正定时，SVD 才会退化成特征值分解的形式；而在更普遍的情况下两者并不相同[^4]。

import numpy as np

def svd(A):
    """
    计算矩阵A的奇异值分解
    
    参数:
        A : 输入矩阵
        
    返回:
        U : 左奇异矩阵
        sigma : 奇异值数组
        VT : 右奇异矩阵转置后的形式
    """
    
    # 计算AA'和A'A的特征值和特征向量
    ATA_eigvals, V = np.linalg.eigh(np.dot(A.T,A))
    AAT_eigvals, UT = np.linalg.eigh(np.dot(A,A.T))

    # 获取奇异值并排序
    sigmas_squared = sorted([val for val in ATA_eigvals if abs(val)>1e-8], reverse=True)
    sigmas = [np.sqrt(x) for x in sigmas_squared]

    # 构建Σ矩阵
    Sigma = np.zeros((len(sigmas), len(sigmas)))
    np.fill_diagonal(Sigma, sigmas)

    return UT[:,:len(sigmas)], Sigma, V[:len(sigmas)].T


if __name__ == "__main__":
    test_matrix = [[1,-2],[4,5]]
    U,s,Vt = svd(test_matrix)
    print("U:\n",U,"\ns:\n",s,"V transpose:\n",Vt)