趁旧智秃1. 引言
复习上一篇文章《最小二乘问题详解1：线性最小二乘》中的知识，对于一个线性问题模型：
f
(
x
;
θ
)
=
A
θ
那么线性最小二乘问题可以表达为求一组待定值
θ
，使得残差的平方和最小：
min
θ
∥
A
θ
?
b
∥
2
本质上是求解超定线性方程组：
A
θ
=
b
具体的线性最小二乘解是：
θ
?
=
(
A
T
A
)
?
1
A
T
b
(1)
2. 求解
2.1 问题
虽然线性最小二乘解已经给出，但是并不意味着在实际的数值计算中就能按照式(1)来进行求解。一个典型的问题就是求逆矩阵：在工程实践和数值计算中，直接求解逆矩阵通常是一个性能消耗大且可能不精确的操作，应该尽量避免。举例来说，我们按照大学本科《线性代数》课程中的方法写程序来求解一个逆矩阵，假设使用伴随矩阵法：
A
?
1
=
1
det
(
A
)
?
adj
(
A
)
其中：
det
(
A
)
是矩阵
A
的行列式。
adj
(
A
)
是
A
的伴随矩阵。
为了求解伴随矩阵
adj
(
A
)
：
求代数余子式 (Cofactor)：对于矩阵
A
中的每一个元素
a
i
j
，计算其代数余子式
C
i
j
。
代数余子式
C
i
j
=
(
?
1
)
i
+
j
?
M
i
j
M
i
j
是删去
A
的第
i
行和第
j
列后得到的子矩阵的行列式（称为余子式）。
构造余子式矩阵：将所有代数余子式
C
i
j
按照原来的位置排列，形成一个新矩阵
C
（称为余子式矩阵）。
转置：将余子式矩阵
C
进行转置，得到的矩阵就是伴随矩阵
adj
(
A
)
。
adj
(
A
)
=
C
T
代入公式：将
det
(
A
)
和
adj
(
A
)
代入公式
A
?
1
=
1
det
(
A
)
?
adj
(
A
)
即可。
这里我们大概能估算，使用伴随矩阵法求逆矩阵的理论复杂度是
O
(
n
!
)
，这是一个阶乘级的增长，算法效率非常低。《线性代数》中介绍的另外一种算法高斯消元法也只能达到
O
(
n
3
)
，呈指数级增加。其实效率只是一方面的问题，使用计算机求解的另外一个问题是舍入误差累积：在计算机中，浮点数运算存在固有的舍入误差；求逆过程涉及大量的除法和减法运算，这些误差会在计算过程中不断累积和传播。总而言之，使用通解求解逆矩阵，可能存在不精确且性能消耗大的问题。
2.2 QR分解
那么不使用逆矩阵怎么办呢？我们需要注意的是，最小二乘问题的本质是求解，而不是求逆矩阵，因此关键是要求解正规方程：
A
T
A
θ
=
A
T
b
对矩阵
A
作QR分解：
A
=
Q
1
R
其中:
Q
1
∈
R
m
×
n
列正交，满足
Q
T
1
Q
1
=
I
n
；
R
∈
R
n
×
n
是上三角矩阵，如果
A
列满秩，则
R
的对角元均非零，可逆。
那么把
A
=
Q
1
R
代入正规方程，得到：
(
Q
1
R
)
T
(
Q
1
R
)
x
=
(
Q
1
R
)
T
b
左边整理，因为
Q
T
1
Q
1
=
I
n
：
R
T
Q
T
1
Q
1
R
x
=
R
T
R
x
右边为
R
T
Q
T
1
b
因此正规方程等价于
R
T
R
x
=
R
T
(
Q
T
1
b
)
若
R
可逆（即
A
满秩，
rank
(
A
)
=
n
），则
R
T
也可逆。左右两边左乘
(
R
T
)
?
1
，得到：
R
x
=
Q
T
1
b
.
令
c
=
Q
?
1
b
（这是一个长度为
n
的向量），于是我们得到一个简单的上三角线性系统：
R
x
=
c
,
c
=
Q
T
1
b
这就是QR方法把正规方程化简得到的核心结果：只需解上三角方程
R
x
=
Q
T
1
b
。
以上只是对
A
列满秩的情况做了推导，如果
A
列满秩，那么QR分解可以表示为
x
=
R
?
1
Q
?
1
b
；如果
A
列不满秩（
R
奇异），需要使用列主元QR方法对
R
T
R
x
=
R
T
(
Q
T
1
b
)
进行求解，或者干脆使用下面要介绍的SVD分解（奇异值分解）法。
2.3 SVD分解
另外一种求解的方法是SVD分解。对任意矩阵
A
，存在奇异值分解：
A
=
U
Σ
V
T
其中:
U
∈
R
m
×
m
为正交（列为左奇异向量），
V
∈
R
n
×
n
为正交（列为右奇异向量），
Σ
∈
R
m
×
n
为“对角块”矩阵，通常写成
Σ
=
[
Σ
r
0
0
0
]
其中
Σ
r
=
diag
(
σ
1
,
…
,
σ
r
)
，
(
σ
1
≥
σ
2
≥
?
≥
σ
r
>
0
)
，
r
=
rank
(
A
)
。
将SVD代入正规方程，先计算
A
?
A
：
A
T
A
=
(
U
Σ
V
T
)
T
(
U
Σ
V
T
)
=
V
Σ
T
U
T
U
Σ
V
T
=
V
(
Σ
T
Σ
)
V
T
.
注意
U
T
U
=
I
。而
Σ
T
Σ
是
n
×
n
的对角块矩阵，其非零对角元就是
σ
2
i
(
i
=
1..
r
)
，其余为零。
同样的，计算
A
T
b
：
A
T
b
=
V
Σ
T
U
T
b
.
于是正规方程变为：
V
(
Σ
T
Σ
)
V
T
x
=
V
Σ
T
U
T
b
.
两边左乘
V
T
，因为
V
正交，
V
T
V
=
I
，得到：
(
Σ
T
Σ
)
(
V
T
x
)
=
Σ
T
(
U
T
b
)
把
y
=
V
T
x
与
c
=
U
T
b
代入，得到更简单的对角方程：
(
Σ
T
Σ
)
y
=
Σ
T
c
接下来按奇异值分块展开对角方程，先写出
Σ
相关的形状：
Σ
=
[
Σ
r
0
0
0
]
,
Σ
?
Σ
=
[
Σ
2
r
0
0
0
]
对
y
和
c
也做相应分块：
y
=
[
y
1
y
2
]
,
c
=
[
c
1
c
2
]
其中
y
1
,
c
1
∈
R
r
对应非零奇异值，
y
2
,
c
2
对应奇异值为0的部分（维度
n
?
r
）,代入得到分块方程:
[
Σ
2
r
0
0
0
]
[
y
1
y
2
]
=
[
Σ
r
0
0
0
]
[
c
1
c
2
]
即等价于两组方程：
Σ
2
r
y
1
=
Σ
r
c
1
,
0
=
0
?
c
2
(
无约束／自由分量
)
由于
Σ
r
为对角且可逆，第一式等价于
Σ
r
y
1
=
c
1
?
y
1
=
Σ
?
1
r
c
1
.
而
y
2
（对应零奇异值的分量）在正规方程中不受约束——这反映了在列秩不足时普通最小二乘解不是唯一的（可以在零空间方向任意加解）。为得到最小范数解（惯常的选择），取
y
2
=
0
。
最后回到
x
的求解，对于
y
有：
y
=
[
Σ
?
1
r
c
1
0
]
将
c
1
与
c
=
U
?
b
关系代回：
y
=
[
Σ
?
1
r
0
0
0
]
U
T
b
由于
y
=
V
T
x
，于是：
x
=
V
y
=
V
[
Σ
?
1
r
0
0
0
]
U
T
b
定义
Σ
+
为将非零奇异值取倒数后转置得到的伪逆矩阵（对角块为
Σ
?
1
r
，其余为0)，则
x
+
=
V
Σ
+
U
T
b
这就是 最小二乘的 Moore–Penrose 伪逆解：
若
A
列满秩，则为唯一最小二乘解，由于那么
Σ
+
=
Σ
?
1
，SVD求解公式退化为常见的
x
=
V
Σ
?
1
U
T
b
若秩亏，它给出 在所有最小二乘解中范数最小的那个（minimum-norm solution）。
2.4 比较
从以上论述可以看到，SVD分解稳定且能处理秩亏的情况，但比QR分解慢，复杂度高，通常
O
(
m
n
2
)
；QR分解在列满秩、条件数不是太差时更快；若需要判定秩或求最小范数解，SVD是首选。
3. 补充
在最后补充一些基础知识，也是笔者很感兴趣的一点，那就是为什么一个矩阵A可以进行QR分解或者SVD分解呢？
3.1 QR分解
QR分解其实非常好理解，它的本质其实就是大学本科《线性代数》课程中提到的施密特（Gram–Schmidt）正交化。我们先复习一下施密特正交化相关的知识。
设有一组线性无关向量
a
1
,
a
2
,
…
,
a
n
∈
R
m
我们想把它们变成一组正交（再归一化后变成标准正交）的向量
q
1
,
q
2
,
…
,
q
n
。具体步骤如下：
取第一个向量，归一化：
q
1
=
a
1
|
a
1
|
对第 2 个向量，先减去在
q
1
上的投影：
~
q
2
=
a
2
?
(
q
T
1
a
2
)
q
1
然后归一化：
q
2
=
~
q
2
|
~
q
2
|
对第 3 个向量，减去它在前两个正交向量上的投影：
~
q
3
=
a
3
?
(
q
T
1
a
3
)
q
1
?
(
q
T
2
a
3
)
q
2
然后归一化：
q
3
=
~
q
3
|
~
q
3
|
一般地，对第
j
个向量：
~
q
j
=
a
j
?
j
?
1
∑
i
=
1
(
q
T
i
a
j
)
q
i
,
q
j
=
~
q
j
|
~
q
j
|
这样得到的
q
i
就是标准正交基，且每个
q
j
只用到了前
j
?
1
个。
现在把矩阵
A
看成由列向量组成：
A
=
[
a
1
,
a
2
,
…
,
a
n
]
∈
R
m
×
n
.
把施密特正交化写成矩阵形式，我们得到一组正交向量：
Q
1
=
[
q
1
,
q
2
,
…
,
q
n
]
∈
R
m
×
n
,
Q
T
1
Q
1
=
I
n
.
同时，原向量可以写成：
a
j
=
j
∑
i
=
1
r
i
j
q
i
其中：
r
i
j
=
q
T
i
a
j
把这些关系拼成矩阵形式：
A
=
Q
1
R
其中：
R
=
(
r
i
j
)
是
n
×
n
上三角矩阵，因为第
j
列只用到前
j
个
q
i
。
Q
1
的列正交，所以
Q
T
1
Q
1
=
I
。
3.2 SVD分解
SVD分解其实也非常有意思，同样也可以顺着《线性代数》中基础知识来进行推导。首先复习一下特征值和特征向量。对于一个方阵
A
∈
R
n
×
n
，如果存在一个非零向量
v
∈
R
n
和一个实数
λ
，使得：
A
v
=
λ
v
那么：
λ
称为 特征值（eigenvalue）
v
称为对应于
λ
的 特征向量（eigenvector）
接下来复习一下什么叫做对角化。如果一个
n
×
n
矩阵
A
可以写成：
A
=
P
D
P
?
1
其中：
D
是一个对角矩阵（只有对角线上有元素）
P
是一个可逆矩阵
我们就说
A
是 可对角化的。
而且通常：
D
的对角元素是
A
的特征值：
D
=
diag
(
λ
1
,
…
,
λ
n
)
P
的列是对应的特征向量
扫一扫
1364
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
