当前位置:   article > 正文

Python奇异值分解_奇异值分解python

奇异值分解python

对于一个矩阵A, 它的奇异值分解

A = U Σ V − 1 A=U\Sigma V^{-1} A=UΣV1
U U U为左奇异矩阵, V V V是右奇异矩阵且都是正交阵, 即 U U ⊤ = I UU^{\top}=I UU=I. Σ \Sigma Σ为奇异值构成的对角阵
example I
A = ( 0 1 1 1 1 0 ) 3 × 2 A=

(011110)
_{3\times2} A=0111103×2
U U U的计算方法
左乘矩阵 A A ⊤ AA^{\top} AA做特征值分解
A ∗ = A A ⊤ ( 0 1 1 1 1 0 ) 3 × 2 ( 0 1 1 1 1 0 ) 2 × 3 = ( 1 1 0 1 2 1 0 1 1 ) 3 × 3 A^*= AA^{\top}
(011110)
_{3\times2}
(011110)
_{2\times3}=
(110121011)
_{3\times3}
A=AA0111103×2(011110)2×3=1101210113×3

求解左乘矩阵 A A ⊤ AA^{\top} AA的特征值和特征向量. 由 A x = λ x Ax=\lambda x Ax=λx

A x − λ x = 0 → ( A − λ I ) x = 0 Ax-\lambda x=0 \to (A-\lambda I)x=0 Axλx=0(AλI)x=0

其中 I I I为单位阵, 假设 x x x不为零, 则有 det ⁡ ( A − λ I ) = 0 \det(A-\lambda I)=0 det(AλI)=0 ∣ A − λ I ∣ = 0 |A-\lambda I|=0 AλI=0, ∣ ⋅ ∣ |\cdot| 为求行列式. 此时有

det ⁡ ( A ∗ − λ I ) = ∣ ( 1 1 0 1 2 1 0 1 1 ) − λ ( 1 0 0 0 1 0 0 0 1 ) ∣ = 0 = ∣ ( 1 − λ 1 0 1 2 − λ 1 0 1 1 − λ ) ∣ = 0 = ∣ 1 − λ 1 0 1 2 − λ 1 0 1 1 − λ ∣ = 0 → ∣ a 1 b 1 c 1 a 2 b 2 c 2 a 3 b 3 c 3 ∣

det(AλI)=|(110121011)λ(100010001)|=0=|(1λ1012λ1011λ)|=0=|1λ1012λ1011λ|=0|a1b1c1a2b2c2a3b3c3|
det(AλI)=110121011λ100010001=0=1λ1012λ1011λ=0=1λ1012λ1011λ=0a1a2a3b1b2b3c1c2c3

则有三阶行列式计算方法:

a 1 ⋅ b 2 ⋅ c 3 + b 1 ⋅ c 2 ⋅ a 3 + c 1 ⋅ a 2 ⋅ b 3 − a 3 ⋅ b 2 ⋅ c 1 − b 3 ⋅ c 2 ⋅ a 1 − c 3 ⋅ a 2 ⋅ b 1 a_1\cdot b_2\cdot c_3+b_1\cdot c_2 \cdot a_3+c_1 \cdot a_2\cdot b_3-a_3\cdot b_2\cdot c_1-b_3\cdot c_2\cdot a_1-c_3\cdot a_2\cdot b_1 a1b2c3+b1c2a3+c1a2b3a3b2c1b3c2a1c3a2b1

从而有特征值 λ 1 = 3 , λ 2 = 1 , λ 3 = 0 \lambda_1=3,\lambda_2=1,\lambda_3=0 λ1=3,λ2=1,λ3=0, 对应特征向量为
λ 1 : A ∗ x = 3 x \lambda_1:A^*x=3x λ1:Ax=3x
→ { x 1 + x 2 + 0 = 3 x 1 x 1 + 2 x 2 + x 3 = 3 x 2 0 + x 2 + x 3 = 3 x 3 → 基 础 解 { x 1 = 1 x 2 = 2 x 3 = 1 → 特 征 向 量 ( 1 2 1 ) → 单 位 化 ( 1 / 6 2 / 6 1 / 6 )

{x1+x2+0=3x1x1+2x2+x3=3x20+x2+x3=3x3{x1=1x2=2x3=1(121)(1/62/61/6)
\\ x1+x2+0=3x1x1+2x2+x3=3x20+x2+x3=3x3x1=1x2=2x3=1 121 1/6 2/6 1/6
λ 2 : A ∗ x = x \lambda_2:A^*x=x λ2:Ax=x
→ { x 1 + x 2 + 0 = x 1 x 1 + 2 x 2 + x 3 = x 2 0 + x 2 + x 3 = x 3 → 基 础 解 { x 1 = 1 x 2 = 0 x 3 = − 1 → 特 征 向 量 ( 1 0 − 1 ) → 单 位 化 ( 1 / 2 0 − 1 / 2 )
{x1+x2+0=x1x1+2x2+x3=x20+x2+x3=x3{x1=1x2=0x3=1(101)(1/201/2)
\\
x1+x2+0=x1x1+2x2+x3=x20+x2+x3=x3x1=1x2=0x3=1 101 1/2 01/2

λ 3 : A ∗ x = 0 x \lambda_3:A^*x=0x λ3:Ax=0x
→ { x 1 + x 2 + 0 = 0 x 1 + 2 x 2 + x 3 = 0 0 + x 2 + x 3 = 0 → 基 础 解 { x 1 = 1 x 2 = − 1 x 3 = 1 → 特 征 向 量 ( 1 − 1 1 ) → 单 位 化 ( 1 / 3 − 1 / 3 1 / 3 )
{x1+x2+0=0x1+2x2+x3=00+x2+x3=0{x1=1x2=1x3=1(111)(1/31/31/3)
\\
x1+x2+0=0x1+2x2+x3=00+x2+x3=0x1=1x2=1x3=1 111 1/3 1/3 1/3

计算非 0 0 0奇异值 σ i = λ i \sigma_i=\sqrt{\lambda_i} σi=λi σ 1 = 3 , σ 2 = 1 \sigma_1=\sqrt{3},\sigma_2=1 σ1=3 ,σ2=1
组合特征向量(stack)
U = ( v 1 v 2 v 2 ) = ( 1 / 6 1 / 2 1 / 3 2 / 6 0 − 1 / 3 1 / 6 − 1 / 2 1 / 3 ) U=
(v1v2v2)
=
(1/61/21/32/601/31/61/21/3)
U=(v1v2v2)=1/6 2/6 1/6 1/2 01/2 1/3 1/3 1/3

V V V的计算方法:
右乘矩阵 A ⊤ A A^{\top}A AA做特征值分解
A ∗ = A ⊤ A = ( 0 1 1 1 1 0 ) 2 × 3 ( 0 1 1 1 1 0 ) 3 × 2 = ( 2 1 1 2 ) 2 × 2 A^*=A^{\top}A=
(011110)
_{2\times3}
(011110)
_{3\times2}=
(2112)
_{2\times2}
A=AA=(011110)2×30111103×2=(2112)2×2

求解左乘矩阵 A A ⊤ AA^{\top} AA的特征值和特征向量. 由 A x = λ x Ax=\lambda x Ax=λx
A x − λ x = 0 → ( A − λ I ) x = 0 Ax-\lambda x=0 \to (A-\lambda I)x=0 Axλx=0(AλI)x=0
其中 I I I为单位阵, 假设 x x x不为零, 则有 det ⁡ ( A − λ I ) = 0 \det(A-\lambda I)=0 det(AλI)=0 ∣ A − λ I ∣ = 0 |A-\lambda I|=0 AλI=0, ∣ ⋅ ∣ |\cdot| 为求行列式. 此时有
det ⁡ ( A ∗ − λ I ) = ∣ ( 2 1 1 2 ) − λ ( 1 0 0 1 ) ∣ = 0 = ∣ ( 2 − λ 1 1 2 − λ ) ∣ = 0 = ∣ 2 − λ 1 1 2 − λ ∣ = 0 → ( 2 − λ ) 2 − 1 = 0
det(AλI)=|(2112)λ(1001)|=0=|(2λ112λ)|=0=|2λ112λ|=0(2λ)21=0
det(AλI)=(2112)λ(1001)=0=(2λ112λ)=0=2λ112λ=0(2λ)21=0

有特征值 λ 1 = 3 , λ 2 = 1 \lambda_1=3,\lambda_2=1 λ1=3,λ2=1, 对应特征向量为
λ 1 : A ∗ x = 3 x \lambda_1:A^*x=3x λ1:Ax=3x
→ { 2 x 1 + x 2 = 3 x 1 x 1 + 2 x 2 = 3 x 2 → 基 础 解 { x 1 = 1 x 2 = 1 → 特 征 向 量 ( 1 1 ) → 单 位 化 ( 1 / 2 1 / 2 )
{2x1+x2=3x1x1+2x2=3x2{x1=1x2=1(11)(1/21/2)
\\
{2x1+x2=3x1x1+2x2=3x2{x1=1x2=1 (11) (1/2 1/2 )

λ 2 : A ∗ x = 1 x \lambda_2:A^*x=1x λ2:Ax=1x
→ { 2 x 1 + x 2 = x 1 x 1 + 2 x 2 = x 2 → 基 础 解 { x 1 = − 1 x 2 = 1 → 特 征 向 量 ( − 1 1 ) → 单 位 化 ( − 1 / 2 1 / 2 )
{2x1+x2=x1x1+2x2=x2{x1=1x2=1(11)(1/21/2)
\\
{2x1+x2=x1x1+2x2=x2{x1=1x2=1 (11) (1/2 1/2 )

计算非 0 0 0奇异值 σ i = λ i \sigma_i=\sqrt{\lambda_i} σi=λi σ 1 = 3 , σ 2 = 1 \sigma_1=\sqrt{3},\sigma_2=1 σ1=3 ,σ2=1
组合特征向量
V = ( 1 / 2 − 1 / 2 1 / 2 1 / 2 ) V=
(1/21/21/21/2)
V=(1/2 1/2 1/2 1/2 )

将奇异值组合成为对角矩阵
Σ = ( 3 0 0 1 ) \Sigma=
(3001)
Σ=(3 001)

finally
A = U Σ V ⊤ = ( 1 / 6 1 / 2 1 / 3 2 / 6 0 − 1 / 3 1 / 6 − 1 / 2 1 / 3 ) 3 × 3 ( 3 0 0 1 0 0 ) 3 × 2 ( 1 / 2 1 / 2 − 1 / 2 1 / 2 ) 2 × 2 A = U\Sigma V^{\top}=
(1/61/21/32/601/31/61/21/3)
_{3\times3}
(300100)
_{3\times2}
(1/21/21/21/2)
_{2\times2}
A=UΣV=1/6 2/6 1/6 1/2 01/2 1/3 1/3 1/3 3×33 000103×2(1/2 1/2 1/2 1/2 )2×2

Python验证

U = np.mat([[1 / np.sqrt(6), 1 / np.sqrt(2), 1 / np.sqrt(3)],
            [2 / np.sqrt(6), 0, -1 / np.sqrt(3)],
            [1 / np.sqrt(6), -1 / np.sqrt(2), 1 / np.sqrt(3)]])
Sigma = np.mat([[np.sqrt(3), 0], [0, 1], [0, 0]])
V = np.mat([[1 / np.sqrt(2), 1 / np.sqrt(2)],
            [-1 / np.sqrt(2), 1 / np.sqrt(2)]])
A = np.mat([[0, 1], [1, 1], [1, 0]])

U.dot(Sigma).dot(np.linalg.inv(V))
>>>
matrix([[ 1., -0.],
        [ 1., -1.],
        [ 0., -1.]])
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13

计算结果显然与真实的 A A A不符合, why?. 这是因为在计算左乘、右乘矩阵的特征向量时是不只有一组解的. 例如左乘矩阵 A A ⊤ AA^{\top} AA
λ 1 : A ∗ x = 3 x \lambda_1:A^*x=3x λ1:Ax=3x
→ { x 1 + x 2 + 0 = 3 x 1 x 1 + 2 x 2 + x 3 = 3 x 2 0 + x 2 + x 3 = 3 x 3 → 基 础 解 { x 1 = 1 x 2 = 2 x 3 = 1 → 特 征 向 量 ( 1 2 1 ) → 单 位 化 ( 1 / 6 2 / 6 1 / 6 ) → 另 一 组 解 { x 1 = − 1 x 2 = − 2 x 3 = − 1 → 特 征 向 量 ( − 1 − 2 − 1 ) → 单 位 化 ( − 1 / 6 − 2 / 6 − 1 / 6 )

{x1+x2+0=3x1x1+2x2+x3=3x20+x2+x3=3x3{x1=1x2=2x3=1(121)(1/62/61/6){x1=1x2=2x3=1(121)(1/62/61/6)
\\ x1+x2+0=3x1x1+2x2+x3=3x20+x2+x3=3x3x1=1x2=2x3=1 121 1/6 2/6 1/6 x1=1x2=2x3=1 121 1/6 2/6 1/6
右乘矩阵 A ⊤ A A^{\top}A AA
λ 1 : A ∗ x = 3 x \lambda_1:A^*x=3x λ1:Ax=3x
→ { 2 x 1 + x 2 = 3 x 1 x 1 + 2 x 2 = 3 x 2 → 基 础 解 { x 1 = 1 x 2 = 1 → 特 征 向 量 ( 1 1 ) → 单 位 化 ( 1 / 2 1 / 2 ) → 另 一 组 解 { x 1 = − 1 x 2 = − 1 → 特 征 向 量 ( − 1 − 1 ) → 单 位 化 ( − 1 / 2 − 1 / 2 )
{2x1+x2=3x1x1+2x2=3x2{x1=1x2=1(11)(1/21/2){x1=1x2=1(11)(1/21/2)
\\
{2x1+x2=3x1x1+2x2=3x2{x1=1x2=1 (11) (1/2 1/2 ){x1=1x2=1 (11) (1/2 1/2 )

这样得到的 U , V U,V U,V就与之前的有所不同, 再次代入python进行验证

U = np.mat([[-1 / np.sqrt(6), 1 / np.sqrt(2), 1 / np.sqrt(3)],
            [-2 / np.sqrt(6), 0, -1 / np.sqrt(3)],
            [-1 / np.sqrt(6), -1 / np.sqrt(2), 1 / np.sqrt(3)]]) 
                                 #因解的不同导致U的第一列发生了符号上的变化
Sigma = np.mat([[np.sqrt(3), 0], [0, 1], [0, 0]])
V = np.mat([[-1 / np.sqrt(2), -1 / np.sqrt(2)],
            [-1 / np.sqrt(2), 1 / np.sqrt(2)]])  
                                  #因解的不同导致V转置的符号发生了变化
A = np.mat([[0, 1], [1, 1], [1, 0]])

U.dot(Sigma).dot(np.linalg.inv(V))
>>>
matrix([[0., 1.],
        [1., 1.],
        [1., 0.]])
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15

此时与示例中的 A A A保持了一致. 所以 A A A的奇异值分解应当为
A = U Σ V ⊤ = ( − 1 / 6 1 / 2 1 / 3 − 2 / 6 0 − 1 / 3 − 1 / 6 − 1 / 2 1 / 3 ) 3 × 3 ( 3 0 0 1 0 0 ) 3 × 2 ( − 1 / 2 − 1 / 2 − 1 / 2 1 / 2 ) 2 × 2 A = U\Sigma V^{\top}=

(1/61/21/32/601/31/61/21/3)
_{3\times3}
(300100)
_{3\times2}
(1/21/21/21/2)
_{2\times2} A=UΣV=1/6 2/6 1/6 1/2 01/2 1/3 1/3 1/3 3×33 000103×2(1/2 1/2 1/2 1/2 )2×2
直接Python实现

np.set_printoptions(suppress=True)
C = np.mat([[0, 1], [1, 1], [1, 0]])
c1, c2 = np.linalg.eig(C.dot(C.T))
d1, d2 = np.linalg.eig(C.T.dot(C))
print("c1:", c1)
print("c2:", c2)
print("d1:", d1)
print("d2:", d2)

U, Sigma, V = np.linalg.svd(C,full_matrices=False)
print("U:", U)
print("Sigma:", Sigma)
print("V:", V)
print(U*np.diag(Sigma)*V)

>>>
左乘矩阵特征根: [ 3.  1. -0.]
左乘矩阵特征向量: [[-0.40824829  0.70710678  0.57735027]
                [-0.81649658  0.         -0.57735027]
                [-0.40824829 -0.70710678  0.57735027]]
    
右乘矩阵特征根: [3. 1.]
右乘矩阵特征向量: [[ 0.70710678 -0.70710678]
                [ 0.70710678  0.70710678]]
    
U: [[-0.40824829  0.70710678]
    [-0.81649658  0.        ]
    [-0.40824829 -0.70710678]]
    
Sigma: [1.73205081 1.        ]
    
V: [[-0.70710678 -0.70710678]
    [-0.70710678  0.70710678]]
    
A: [[ 0.  1.]
    [ 1.  1.]
    [ 1. -0.]]
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37

计算结果与 A A A一致
但反过来回想一下, 为什么说奇异值分解首先是对左、右乘矩阵做特征值分解?因为特征值分解只适用于方阵, 左乘矩阵 A A ⊤ AA^{\top} AA、右乘矩阵 A ⊤ A A^{\top}A AA恰都是方阵, 所以适用特征值分解. 现在给出特征值分解式
A = Q Λ Q − 1 A=Q\Lambda Q^{-1} A=QΛQ1
上式中的 Q Q Q是由 A A A的特征向量构成的, that is to say Q = U Q=U Q=U, 实际上我们并不需要计算 Q − 1 Q^{-1} Q1, 因为在奇异值分解的定义中规定了 U U U是正交阵, 所以 Q = U Q=U Q=U也是正交阵, 所以 Q − 1 = Q ⊤ Q^{-1}=Q^{\top} Q1=Q, 这既是正交阵的性质, we also can improve this in python:

#借用上述代码中的变量c2, 因为c2即左乘矩阵的特征向量阵=Q的定义
np.linalg.inv(c2)
>>>
matrix([[-0.40824829, -0.81649658, -0.40824829],
        [ 0.70710678,  0.        , -0.70710678],
        [ 0.57735027, -0.57735027,  0.57735027]])

c2.T
>>>
matrix([[-0.40824829, -0.81649658, -0.40824829],
        [ 0.70710678,  0.        , -0.70710678],
        [ 0.57735027, -0.57735027,  0.57735027]])
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12

所对于左乘矩阵的特征值分解为

c2.dot(np.diag(c1)).dot(np.linalg.inv(c2))
>>>
matrix([[ 1.,  1., -0.],
        [ 1.,  2.,  1.],
        [-0.,  1.,  1.]])
  • 1
  • 2
  • 3
  • 4
  • 5

即左乘矩阵, 对于右乘矩阵同理可得

除此之外再考虑一个问题, 为什么在有些计算中考虑了 U U U的秩即 r a n k ( U ) \rm{rank}(U) rank(U), 此时奇异值分解为

A = U Σ V ⊤ = ( − 1 / 6 1 / 2 − 2 / 6 0 − 1 / 6 − 1 / 2 ) 3 × 2 ( 3 0 0 1 ) 2 × 2 ( − 1 / 2 − 1 / 2 − 1 / 2 1 / 2 ) 2 × 2 A = U\Sigma V^{\top}=

(1/61/22/601/61/2)
_{3\times2}
(3001)
_{2\times2}
(1/21/21/21/2)
_{2\times2} A=UΣV=1/6 2/6 1/6 1/2 01/2 3×2(3 001)2×2(1/2 1/2 1/2 1/2 )2×2
通过python进行验证

U = np.mat([[-1 / np.sqrt(6), 1 / np.sqrt(2)],
            [-2 / np.sqrt(6), 0],
            [-1 / np.sqrt(6), -1 / np.sqrt(2)]])
Sigma = np.mat([[np.sqrt(3), 0], [0, 1]])
V = np.mat([[-1 / np.sqrt(2), -1 / np.sqrt(2)],
            [-1 / np.sqrt(2), 1 / np.sqrt(2)]])
A = np.mat([[0, 1], [1, 1], [1, 0]])

U.dot(Sigma).dot(np.linalg.inv(V))
>>>
matrix([[0., 1.],
        [1., 1.],
        [1., 0.]])
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13

可以看到考虑 r a n k ( U ) = 2 \rm{rank}(U)=2 rank(U)=2 2 2 2阶奇异值分解与 3 3 3阶的计算结果是完全相同的. 这是由于秩反映了一个矩阵的极大线性无关组数量, 由 r a n k ( U ) = 2 \rm{rank}(U)=2 rank(U)=2可得对于三级矩阵 U U U仅有两个极大线性无关组, 从本质上来说这意味着 U U U中的第三个列向量并没有提供任何有用的信息, 且其对应特征值也为 0 0 0, 所以删除这一列向量不会导致信息损失, 保留也不会提供任何信息, 从奇异值矩阵
Σ = ( 3 0 0 1 0 0 ) 3 × 2 \Sigma=

(300100)
_{3\times2} Σ=3 000103×2
也能看出新加入的第三行全为 0 0 0, 同样不提供任何信息. 但考虑秩的情况能有一个好处: 对高维矩阵进行奇异值分解时能极大地提高计算效率.

  • 思考

根据矩阵秩的定义, 矩阵的行秩 = = =列秩 = = =矩阵的秩, 我们可以将无用的列向量删除来简化计算, 那是否意味着我们也通过删除行向量的方式得到相同的奇异值分解结果?

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/小小林熬夜学编程/article/detail/713917
推荐阅读
相关标签
  

闽ICP备14008679号