⚠️ Alpha内测版本警告:此为早期内部构建版本,尚不完整且可能存在错误,欢迎大家提Issue反馈问题或建议。
Skip to content

第15章 奇异值分解

习题15.1

  试求矩阵

A=[120202]

的奇异值分解。

解答:

解答思路:

  1. 给出奇异值分解的定义
  2. 给出奇异值分解的步骤
  3. 使用numpy实现奇异值分解
  4. 自编程实现奇异值分解
  5. 推导计算奇异值分解

解答步骤:

第1步:奇异值分解

  根据书中第15.1节的奇异值分解的定义:

定义15.1(奇异值分解) 矩阵的奇异值分解是指,将一个非零的m×n实矩阵AARm×n,表示为以下三个实矩阵乘积形式的运算,即进行矩阵的因子分解:

A=UΣVT

其中Um阶正交矩阵,Vn阶正交矩阵,Σ是由降序排列的非负的对角线元素组成的m×n矩形对角矩阵,

UUT=IVVT=IΣ=diag(σ1,σ2,,σp)σ1σ2σp0p=min(m,n)

UΣVT称为矩阵A的奇异值分解,σi称为矩阵A的奇异值,U的列向量称为左奇异向量,V的列向量称为右奇异向量。

第2步:奇异值分解的一般步骤

  根据书中第15.2节的奇异值分解的计算:

给定m×n矩阵A,可以写出矩阵奇异值分解的计算过程。
(1)求ATA的特征值和特征向量。
计算对称矩阵W=ATA
求解特征方程

(WλI)x=0

得到特征值λi,并将特征值由大到小排列

λ1λ2λn0

将特征值λi (i=1,2,,n)代入特征方程求得对应的特征向量。
(2)求n阶正交矩阵V
将特征向量单位化,得到单位特征向量v1,v2,,vn,构成n阶正交矩阵V

V=[v1 v2  vn]

(3)求m×n对角矩阵Σ
计算A的奇异值

σi=λi,i=1,2,,n

构造m×n矩形对角矩阵Σ,主对角线元素是奇异值,其余元素是零

Σ=diag(σ1,σ2,,σn)

(4)求m阶正交矩阵U
A的前r个正奇异值,令

uj=1σjAvj,j=1,2,,r

得到

U1=[u1 u2  ur]

AT的零空间的一组标准正交基{ur+1,ur+2,,um},令

U2=[ur+1 ur+2  um]

并令

U=[U1 U2]

(5)得到奇异值分解

A=UΣVT

第3步:使用numpy实现奇异值分解

python
import numpy as np

A = np.array([[1, 2, 0],
              [2, 0, 2]])

# 调用numpy的svd方法
U, S, V = np.linalg.svd(A)

# 设置精度为2
np.set_printoptions(precision=2, suppress=True)

print("U=", U)
print("S=", S)
print("V=", V.T)
Sigma = np.zeros_like(A, float)
np.fill_diagonal(Sigma, S)
calc = np.dot(np.dot(U, Sigma), V)
print("A=", calc)
U= [[ 0.45  0.89]
 [ 0.89 -0.45]]
S= [3. 2.]
V= [[ 0.75 -0.   -0.67]
 [ 0.3   0.89  0.33]
 [ 0.6  -0.45  0.67]]
A= [[ 1.  2.  0.]
 [ 2. -0.  2.]]

第4步:自编程实现奇异值分解

python
import numpy as np
from scipy.linalg import null_space


def my_svd(A):
    m = A.shape[0]

    # (1) 计算对称矩阵 A^T A 的特征值与特征向量,
    W = np.dot(A.T, A)
    # 返回的特征值lambda_value是升序的,特征向量V是单位化的特征向量
    lambda_value, V = np.linalg.eigh(W)
    # 并按特征值从大到小排列
    lambda_value = lambda_value[::-1]
    lambda_value = lambda_value[lambda_value > 0]
    # (2)计算n阶正价矩阵V
    V = V[:, -1::-1]

    # (3) 求 m * n 对角矩阵
    sigma = np.sqrt(lambda_value)
    S = np.diag(sigma) @ np.eye(*A.shape)

    # (4.1) 求A的前r个正奇异值
    r = np.linalg.matrix_rank(A)
    U1 = np.hstack([(np.dot(A, V[:, i]) / sigma[i])[:, np.newaxis] for i in range(r)])
    # (4.2) 求A^T的零空间的一组标准正交基
    U = U1
    if r < m:
        U2 = null_space(A.T)
        U2 = U2[:, r:]
        U = np.hstack([U, U2])

    return U, S, V
python
A = np.array([[1, 2, 0],
              [2, 0, 2]])

np.set_printoptions(precision=2, suppress=True)

U, S, V = my_svd(A)
print("U=", U)
print("S=", S)
print("V=", V)
calc = np.dot(np.dot(U, S), V.T)
print("A=", calc)
U= [[-0.45 -0.89]
 [-0.89  0.45]]
S= [[3. 0. 0.]
 [0. 2. 0.]]
V= [[-0.75 -0.   -0.67]
 [-0.3  -0.89  0.33]
 [-0.6   0.45  0.67]]
A= [[ 1.  2. -0.]
 [ 2. -0.  2.]]

第5步:推导计算奇异值分解

  1. 求解ATA的特征值和特征向量

  计算对称矩阵W=ATA

ATA=[122002][120202]=[524240404]

  特征值和特征向量满足方程

(ATAλI)x=0

  得到齐次方程组

[5λ2424λ0404λ]x=0

  该方程组有非零解的充要条件是

|5λ2424λ0404λ|=0

  解得矩阵 W=ATA 的特征值 λ1=9λ2=4λ3=0 (按从大到小排列)。

  特征值代入对应的齐次方程组,得到对应的单位特征向量 (特征向量的选取不唯一)

v1=[535235435],v2=[02515],v3=[231323]
  1. 求正交矩阵V

  构造正交矩阵V

V=[v1v2v3]=[53502323525134351523]
  1. 求解对角矩阵Σ

  构造对角矩阵Σ

Σ=[300020]
  1. 求正交矩阵U

  基于A的正奇异值 σ1=3σ2=2 计算得到列向量u1,u2

u1=1σ1Av1=13[120202][535235435]=[1525]u2=1σ2Av2=12[120202][02515]=[2515]

  构造正交矩阵U

U=[15252515]
  1. 矩阵A的奇异值分解
A=UΣVT=[15252515][300020][53502323525134351523]T

习题15.2

  试求矩阵

A=[24130000]

的奇异值分解并写出其外积展开式。

解答:

解答思路:

  1. 给出矩阵的外积展开式定义
  2. 使用numpy计算矩阵A的奇异值分解
  3. 根据奇异值分解的结果,写出外积展开式

解答步骤:

第1步:矩阵 A 的外积展开式;

  根据书中第15.3.3节的矩阵的外积展开式:

  矩阵A的奇异值分解UΣVT也可以由外积形式表示。事实上,若将A的奇异值分解看成矩阵UΣVT的乘积,将UΣ按列向量分块,将VT按行向量分块,即得

UΣ=[σ1u1 σ2u2  σnun]VT=[v1Tv2TvnT]

A=σ1u1v1T+σ2u2v2T++σnunvnT

称为矩阵A的外积展开式,其中ukvkTm×n矩阵,是列向量uk和行向量vkT的外积,其第i行第j列元素为uk的第i个元素与vkT的第j个元素的乘积。即

uivjT=[u1iu2iumi][v1j v2j  vnj]=[u1iv1ju1iv2ju1ivnju2iv1ju2iv2ju2ivnjumiv1jumiv2jumivnj]

A的外积展开式也可以写成下面的形式

A=k=1nAk=k=1nσkukvkT

其中Ak=σkukvkTm×n矩阵。上式将矩阵A分解为矩阵的有序加权和。

  根据书中第14章的本章概要中的外积展开式:

任意一个实矩阵A可以由其外积展开式表示

A=σ1u1v1T+σ2u2v2T++σnunvnT

其中ukvkTm×n矩阵,是列向量uk和行向量vkT的外积,σk为奇异值,uk,vkT,σk通过矩阵A的奇异值分解得到。

第2步:计算矩阵A的奇异值分解

python
import numpy as np

A = np.array([[2, 4],
              [1, 3],
              [0, 0],
              [0, 0]])

# 调用numpy的svd方法
U, S, V = np.linalg.svd(A)
np.set_printoptions()

print("U=", U)
print("S=", S)
print("V=", V.T)
U= [[-0.82 -0.58  0.    0.  ]
 [-0.58  0.82  0.    0.  ]
 [ 0.    0.    1.    0.  ]
 [ 0.    0.    0.    1.  ]]
S= [5.46 0.37]
V= [[-0.4  -0.91]
 [-0.91  0.4 ]]

第3步:根据奇异值分解的结果,写出外积展开式

  矩阵A的外积展开式为

A=σ1u1v1T+σ2u2v2T

其中

σ1=5.4649857,σ2=0.36596619u1=[0.817415560.5760484400],u2=[0.576048440.8174155600]v1T=[0.40455358,0.9145143],v2T=[0.9145143,0.40455358]
python
calc = S[0] * np.outer(U[:, 0], V[:, 0]) + S[1] * np.outer(U[:, 1], V[:, 1])
print("A=", calc)
A= [[ 2.  4.]
 [ 1.  3.]
 [-0.  0.]
 [-0.  0.]]

习题15.3

  比较矩阵的奇异值分解与对称矩阵的对角化的异同。

解答:

解答思路:

  1. 给出矩阵的奇异值分解的定义
  2. 给出对称矩阵和方阵的对角化定义
  3. 比较两者的相同点
  4. 比较两者的不同点

解答步骤:

注意: 与书中一致,不考虑复数矩阵。

第1步:矩阵的奇异值分解的定义:

  矩阵的奇异值分解定义,详见书中第271页的定义15.1(奇异值分解)

第2步:对称矩阵和方阵对角化定义;

  根据《实对称矩阵对角化教学的应用案例》(张丽静、刘白羽、申亚男,2019)对称矩阵的对角化:

An阶实对称矩阵,则存在一个正交矩阵P=(p1,p2,,pn),使得

A=PΛPT

其中,Λ为对角元素是A的特征值λ1,λ2,,λn所构成的对角矩阵;p1,p2,,pn为特征值λ1,λ2,,λn对应的单位特征向量。

  根据参考资料:https://www.matongxue.com/parts/4628 中关于方阵的对角化:

如果n阶方阵An个线性无关的特征向量p1,p2,,pn,那么如下矩阵:

P=(p1,p2,,pn)

可以使得:$ A = P \Lambda P^{-1}\Lambda\Lambda = \text{diag} { \lambda_1, \lambda_2, \cdots, \lambda_n }$,特征值 λ1,λ2,,λn为特征向量p1,p2,,pn对应的特征值,该过程为对角化(Diagonalizable)。

第3步:两者的相同点

  1. 两者都可以实现特征分解,都需要求解特征值与特征向量
  2. 实对称矩阵的对角化是方阵对角化的特例,矩阵奇异值分解又是方阵对角化的推广
  3. 实对称矩阵一定可以对角化,矩阵的奇异值分解也一定存在
  4. 都可以被正交矩阵对角化
  5. 本质上都是用不同的基来描述线性变化

第4步:两者的不同点

  1. 矩阵对角化要求矩阵A是一个方阵,矩阵奇异值分解不需要对矩阵A有特殊要求
  2. 奇异值分解中要求对角元素降序排列,对称矩阵的对角化没有这一要求

习题15.4

  证明任何一个秩为1的矩阵可以写成两个向量的外积形式,并给出实例。

解答:

解答思路:

  1. 假设矩阵A的秩为1,写出其奇异值分解形式
  2. 将奇异值分解的矩阵Σ写成外积形式
  3. 将矩阵A的奇异值分解写成外积形式
  4. 给出实例

解答步骤:

第1步:假设矩阵 A 的秩为1,写出其奇异值分解形式

  设矩阵Am×n的秩为1。

  根据书中第15章的定理15.1的奇异值分解基本定理:

  定理15.1(奇异值分解基本定理)A为一m×n实矩阵,ARm×n,则A的奇异值分解存在

A=UΣVT

其中Um阶正交矩阵,Vn阶正交矩阵,Σm×n矩形对角矩阵,其对角线元素非负,且按降序排列。

  则存在矩阵A的奇异值分解

(1)A=UΣVT

第2步:将奇异值分解的矩阵Σ写成外积形式;

  根据矩阵A的秩与对角矩阵Σ的秩之间的关系:

rank(Σ)=rank(A)=1

  根据奇异值分解基本定理,对角矩阵Σ的奇异值是降序排列的,所以其中的非零元素一定位于第一行第一列,可设为

[σ1O(m1)×(n1)]m×n

其中,O(m1)×(n1)表示m1n1列的零矩阵。

  设两向量 a=[1,0,,0]m×1T,b=[σ1,0,,0]n×1T,可得:

(2)Σ=abT

第3步:将矩阵A 的奇异值分解写成外积形式

  将式(2)代入式(1)中可得:

A=UabTVT=(Ua)(Vb)T

其中,Uam×1阶的列向量,Vbn×1阶的列向量。

  所以,矩阵A可以写成向量Ua和向量Vb的外积形式,命题得证。

第4步:给出实例

  使用书中第283页例15.5的矩阵A=[112200]的秩为1,矩阵A的奇异值分解为

A=UΣVT=[1525025150001][1000000][12121212]

  可知:

Σ=[1000000]=[100][10, 0]=abTUa=[15250],Vb=[55]

  所以,(Ua)(Vb)T=[112200]=A

习题15.5

  搜索中的点击数据记录用户搜索时提交的查询语句,点击的网页URL以及点击的次数构成一个二部图,其中一个结点集合{qi}表示查询,另一个结点集合{uj}表示URL,边表示点击关系,边上的权重表示点击次数。图15.2是一个简化的点击数据例。点击数据可以由矩阵表示,试对该矩阵进行奇异值分解,并解释得到的三个矩阵所表示的内容。

解答:

解答思路:

  1. 根据二部图写出矩阵A
  2. 计算矩阵A的奇异值分解
  3. 解释奇异值分解得到的矩阵UΣV代表的内容

解答步骤:

第1步:根据二部图写出矩阵A

A=[0205001000300000100100]

其中,行向量分别对应(q1,q2,q3,q4),列向量分别对应(u1,u2,u3,u4,u5)

第2步:对矩阵 A 进行奇异值分解;

python
import numpy as np

A = np.array([[0, 20, 5, 0, 0],
              [10, 0, 0, 3, 0],
              [0, 0, 0, 0, 1],
              [0, 0, 1, 0, 0]]) 

# 调用numpy的svd方法
U, S, V = np.linalg.svd(A)

print("U=", U)
print("S=", S)
print("V=", V.T)
U= [[ 1.   -0.    0.   -0.01]
 [ 0.    1.    0.   -0.  ]
 [ 0.    0.    1.    0.  ]
 [ 0.01  0.    0.    1.  ]]
S= [20.62 10.44  1.    0.97]
V= [[ 0.    0.96 -0.   -0.    0.29]
 [ 0.97 -0.   -0.   -0.24 -0.  ]
 [ 0.24  0.    0.    0.97  0.  ]
 [ 0.    0.29  0.    0.   -0.96]
 [ 0.    0.    1.    0.    0.  ]]

第3步:解释奇异值分解得到的矩阵UΣV代表的内容

  根据题意,可知搜索中的点击数据记录用户搜索时提交的查询语句为{qi},点击的网页URL为{uj},点击的次数为矩阵的值。根据奇异值分解定理,可知:

  • U表示用户输入的查询语句与网页特征之间的关系矩阵
  • V表示网页与网页特征之间的相似性
  • Σ表示用户输入的查询语句映射到网页的权重

  举个例子:从矩阵U的第一列可以表示查询语句q1与网页特征1的关系密切,通过矩阵V可以表示网页u2与特征1的相似度最高,所以查询 q1的时候点击网页u2的次数是最多的。

  奇异值分解的另一个作用就是提取网页的特征来将用户输入与网页映射到一个低纬度空间中。通过第2步的计算,可以发现σ1,σ2的值比较大,计算截断的矩阵外积展开A=σ1u1v1T+σ2u2v2T,可以发现其与矩阵A大体上是相等的。奇异值下降的速度越快,那么矩阵包含的更多信息就越集中分布在前面几个值比较大的特征上面。

参考文献

【1】张丽静、刘白羽、申亚男. 实对称矩阵对角化教学的应用案例[C]. 大学数学,2019,35(2):116-121
【2】方阵的对角化:https://www.matongxue.com/parts/4628