继续MIT笔记的内容,前面讨论了AX=b的解情况,我着重画了一些解空间的图形。后来我们将这些子空间分为四个部分,即 A A A的列空间、行空间和零空间,以及 A T A^T AT的零空间 ,这四个空间的关系见课本的figure4.2,其实由上次的最后一个例子,我们就能发现这样的一种垂直关系。应该指出,垂直关系是一种很不错的关系,本次主要介绍投影矩阵P和标准正交矩阵Q。
“Follow the rules”
老爷子在讲正交向量和正交子空间时,问道零向量是否与任意向量均正交,他给出了这样一个建议“The one thing about math is you’re supposed to follow the rules.” 这个建议的味道在于,凡事难以决定乃至不可思议时,你应该回归定义,倘若定义中允许它发生,而它也发生了,那它就是合理的,你也该去正视它。所以,零向量是否与任意向量均正交,答案是“sure”。那么来看正交(或言之为空间关系中的垂直)的定义:
若内积空间中两向量的内积为0,则它们正交。
举figure4.2来说,已知零空间的定义为
零空间是在线性映射(即矩阵)的背景下出现的,指像为零的原像空间,即 { x ∣ A x = 0 x| Ax=0 x∣Ax=0} 。
可见零空间的定义中恰好就内嵌了正交的定义。于是看 A X = 0 AX=0 AX=0, A A A通过 X X X表征的零空间与 A A A表征的行空间做内积为0,也就是 A A A的零空间与 A A A的行空间垂直,即figure4.2的左边部分; 再看 A T Y = 0 A^TY=0 ATY=0, A T A^T AT通过 Y Y Y表征的零空间与 A T A^T AT 的行空间(由于转置,也即 A A A 的列空间)垂直,展示为figure4.2的右边部分。
投影矩阵P
什么是投影矩阵?我们结合figure4.6说明这个概念。
从向量入手
先看左边部分,这是一个将向量
b
⃗
\vec{b}
b 投影到
a
⃗
\vec{a}
a 上的问题,结合正交的定义,我们有
e
⃗
⋅
a
⃗
=
0
⇒
(
b
⃗
−
p
⃗
)
⋅
a
⃗
=
0
⇒
b
⃗
⋅
a
⃗
−
p
⃗
⋅
a
⃗
=
0
⇒
p
⃗
⋅
a
⃗
=
b
⃗
⋅
a
⃗
\vec{e} \cdot \vec{a} = 0 \Rightarrow (\vec{b}-\vec{p}) \cdot \vec{a} = 0 \Rightarrow \vec{b} \cdot \vec{a} - \vec{p} \cdot \vec{a} = 0 \Rightarrow \vec{p} \cdot \vec{a} = \vec{b} \cdot \vec{a}
e⋅a=0⇒(b−p)⋅a=0⇒b⋅a−p⋅a=0⇒p⋅a=b⋅a
两边同时右乘
a
⃗
\vec{a}
a ,注意向量乘法不满足结合律,矩阵乘法才满足,但矩阵乘法不满足交换律,向量乘法满足。于是就有
p
⃗
∣
a
∣
2
=
b
⃗
⋅
a
⃗
⋅
a
⃗
⇒
p
⃗
=
b
⃗
⋅
a
⃗
∣
a
∣
2
a
⃗
=
a
⃗
⋅
b
⃗
∣
a
∣
2
a
⃗
\vec{p} |a|^2 = \vec{b} \cdot \vec{a} \cdot \vec{a} \Rightarrow \vec{p} = \frac{\vec{b} \cdot \vec{a}}{|a|^2}\vec{a}= \frac{\vec{a} \cdot \vec{b} }{|a|^2}\vec{a}
p∣a∣2=b⋅a⋅a⇒p=∣a∣2b⋅aa=∣a∣2a⋅ba
我们用矩阵形式表示上述向量式,
a
T
a
a^Ta
aTa表示对
a
⃗
\vec{a}
a 所表示列向量的模值平方,且注意到
a
T
b
a^Tb
aTb也为一数,就有:
p
=
a
T
b
a
T
a
a
=
a
a
T
b
a
T
a
=
a
a
T
a
T
a
b
p = \frac{a^Tb}{a^Ta}a= a\frac{a^Tb}{a^Ta}=\frac{a a^T}{a^Ta}b
p=aTaaTba=aaTaaTb=aTaaaTb
投影矩阵于是起到这样一个作用,找出了向量
b
⃗
\vec{b}
b 在向量
a
⃗
\vec{a}
a 上的分量,也就是
p
=
P
b
=
a
x
^
p = Pb=a\hat{x}
p=Pb=ax^ 。 对于左边部分,投影矩阵即
P
=
a
a
T
a
T
a
P = \frac{a a^T}{a^Ta}
P=aTaaaT
再看矩阵
但看这个过程,我们很容易发现,
a
T
a
a^Ta
aTa 和
a
T
b
a^Tb
aTb均为数这个条件是苛刻的,针对的是一维的列矩阵;对于一般性的矩阵,这两个的结果都是矩阵而非一个数。于是我们继续来看右边部分,将 b 投影到一个用二维列向量表示的平面上。尽管方程变复杂了,但得到方程仍旧是一样的,即用平面中的一组基底与
e
⃗
\vec{e}
e作内积后结果为0。那么对于
A
A
A而言,其两个基底所在的平面即其列空间,若将e视作一个列向量,即解方程:
A
T
e
=
0
A^Te=0
ATe=0
我们从投影矩阵P的角度去产生e:
e
=
b
−
p
=
b
−
P
b
=
b
−
A
x
^
e = b-p=b-Pb=b-A\hat{x}
e=b−p=b−Pb=b−Ax^
反代后,就有:
A
T
(
b
−
A
x
^
)
=
A
T
b
−
A
T
A
x
^
=
0
⇒
A
T
A
x
^
=
A
T
b
⇒
x
^
=
(
A
T
A
)
−
1
A
T
b
⇒
p
=
A
x
^
=
A
(
A
T
A
)
−
1
A
T
b
A^T(b-A\hat{x})=A^Tb-A^TA\hat{x}=0 \Rightarrow A^T A\hat{x} = A^Tb \Rightarrow \hat{x} = (A^TA)^{-1}A^Tb \Rightarrow p = A\hat{x} = A(A^TA)^{-1}A^Tb
AT(b−Ax^)=ATb−ATAx^=0⇒ATAx^=ATb⇒x^=(ATA)−1ATb⇒p=Ax^=A(ATA)−1ATb
最后推得投影矩阵P
P
=
A
(
A
T
A
)
−
1
A
T
P = A(A^TA)^{-1}A^T
P=A(ATA)−1AT
但这种推导,显然基于一个假设,即 “If A has independent columns, then
A
T
A
A^TA
ATA is invertible.” 或言之,对于A而言,存在左逆矩阵。
最小二乘法
一个有趣的关于投影矩阵的用法是导出最小二乘法,其常用于线性回归,我们既可以用一般性的梯度下降法交给计算机迭代导数求解,也可以直接从矩阵出发求解。我们不如来看一个一般性的例子:
A
x
=
[
a
11
a
12
a
21
a
22
a
31
a
32
.
.
.
.
.
.
]
[
C
D
]
=
[
b
1
b
2
b
3
.
.
.
]
=
b
Ax = \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \\ a_{31} & a_{32} \\ . & .\\. & .\\ . &.\end{bmatrix} \begin{bmatrix} C \\ D \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ b_3 \\ . \\ . \\ . \end{bmatrix} =b
Ax=⎣⎢⎢⎢⎢⎢⎢⎡a11a21a31...a12a22a32...⎦⎥⎥⎥⎥⎥⎥⎤[CD]=⎣⎢⎢⎢⎢⎢⎢⎡b1b2b3...⎦⎥⎥⎥⎥⎥⎥⎤=b
A
x
=
b
Ax= b
Ax=b对于数据点个数(方程数)远大于拟合的参数个数(2个未知数)时,显然要么有解要么就无解。按照最小二乘的想法,面对无解的情况时,我们应求下列函数的最小值,即
m
i
n
(
a
11
C
+
a
12
D
−
b
1
)
2
+
(
a
21
C
+
a
22
D
−
b
2
)
2
+
(
a
31
C
+
a
32
D
−
b
3
)
2
+
.
.
.
min \ \ (a_{11}C+a_{12}D-b_1)^2+(a_{21}C+a_{22}D-b_2)^2+(a_{31}C+a_{32}D-b_3)^2+...
min (a11C+a12D−b1)2+(a21C+a22D−b2)2+(a31C+a32D−b3)2+...
我们对C和D分别求偏导(已略去系数2),就有
a
11
(
a
11
C
+
a
12
D
−
b
1
)
+
a
21
(
a
21
C
+
a
22
D
−
b
2
)
+
a
31
(
a
31
C
+
a
32
D
−
b
3
)
+
.
.
.
=
0
a
12
(
a
11
C
+
a
12
D
−
b
1
)
+
a
22
(
a
21
C
+
a
22
D
−
b
2
)
+
a
32
(
a
31
C
+
a
32
D
−
b
3
)
+
.
.
.
=
0
a_{11} (a_{11}C+a_{12}D-b_1)+a_{21}(a_{21}C+a_{22}D-b_2)+a_{31}(a_{31}C+a_{32}D-b_3)+... = 0 \\ a_{12} (a_{11}C+a_{12}D-b_1)+a_{22}(a_{21}C+a_{22}D-b_2)+a_{32}(a_{31}C+a_{32}D-b_3)+... = 0 \\
a11(a11C+a12D−b1)+a21(a21C+a22D−b2)+a31(a31C+a32D−b3)+...=0a12(a11C+a12D−b1)+a22(a21C+a22D−b2)+a32(a31C+a32D−b3)+...=0
我们将两个式子相加,不难发现,C的系数即是
∑
a
i
1
2
+
a
i
1
a
i
2
\sum a_{i1}^2+a_{i1}a_{i2}
∑ai12+ai1ai2,D的系数即是
∑
a
i
2
2
+
a
i
1
a
i
2
\sum a_{i2}^2+a_{i1}a_{i2}
∑ai22+ai1ai2, 而
b
i
b_i
bi前的系数相当于
a
i
1
+
a
i
2
a_{i1}+a_{i2}
ai1+ai2。于是我们可以写出这个方程
A
T
A
[
C
D
]
=
A
T
b
A^TA \begin{bmatrix} C \\ D \end{bmatrix} = A^Tb
ATA[CD]=ATb
于是A的左逆存在时,就有:
[
C
D
]
=
(
A
T
A
)
−
1
A
T
b
\begin{bmatrix} C \\ D \end{bmatrix} =(A^TA)^{-1}A^Tb
[CD]=(ATA)−1ATb
那么问题来了,如何从P导出这个式子呢?我们知道,当P存在时,我们就能将不在A列空间中的b投影到A所在的列空间中,这就能形成有解的保证。也就是说:
P
b
=
A
x
^
⇒
A
(
A
T
A
)
−
1
A
T
b
=
A
x
^
Pb = A \hat{x} \Rightarrow A(A^TA)^{-1}A^Tb = A\hat{x}
Pb=Ax^⇒A(ATA)−1ATb=Ax^
于是我们就有
x
^
=
(
A
T
A
)
−
1
A
T
b
\hat{x} = (A^TA)^{-1}A^Tb
x^=(ATA)−1ATb
从这种角度去看,最小二乘法的实质就是让向量投影到平面上所产生的误差和
∑
e
i
\sum e_i
∑ei(向量和) 最小 。我们可以用matlab验证其正确性,即 先用ployfit函数求解,再用矩阵验证。
% 用matlab自带的ployfit函数拟合曲线
x=[9,13,15,17,18.6,20,23,29,31.7,35];
y=[-8,-6.45,-5.1,-4,-3,-1.95,-1.5,-0.4,0.2,-0.75];
coefficient=polyfit(x,y,1); %用一次函数拟合曲线,想用几次函数拟合,就把n设成那个数
y1=polyval(coefficient,x);
plot(x,y,'o',x,y1,'-');
legend('散点',['拟合曲线y=',num2str(coefficient(1)),'*x',num2str(coefficient(2))]);
% 验证
c = ones(1,length(x)); % 常数列
A = [c;x]'; % 矩阵A
est = inv(A'*A)*A'*y'; % 估计值
text(15,-6,['常数项系数=',num2str(est(1)),' 一次项系数=',num2str(est(2))]);
标准正交矩阵
性质
所谓标准,即组成矩阵的向量组中每个向量模值为1,所谓正交,即向量组中的任意两个向量的内积为0,即一开始就提醒的““Follow the rules””。规定这样一个特殊的矩阵必然有一些优良的性质,由于任意两个向量的内积为0,我们很容易明白其都是线性无关的。根据标准正交矩阵的定义,我们很容易推得下式:
Q
Q
T
=
I
QQ^T = I
QQT=I
因为Q的各列只是与自己是相关的,且模为1,与其余各列都是无关的,内积均为0,所以结果将是一个单位阵。进一步对于方阵而言,我们知道由于各列线性无关,逆必然存在,而逆的定义似乎也就是这里
Q
T
Q^T
QT所在的位置,于是,就有
Q
−
1
=
Q
T
Q^{-1}=Q^{T}
Q−1=QT
标准正交方阵的逆就是其转置,这份礼物让人足够开心。
过程
我们来展示将矩阵标准正交化的方法,即Gram-Schmidt正交化,该步骤分为两步,第一步是正交化,第二步是标准化。理解正交化的原理同样可借助于我们投影矩阵P的那个图,即我们现在不求p二是求正交的e。
借助figure4.6,我们知道e应该是这么求的:
B
=
b
−
P
b
=
b
−
a
(
a
T
a
)
−
1
a
T
b
=
b
−
a
T
b
a
T
a
a
B= b-Pb=b-a(a^Ta)^{-1}a^Tb = b-\frac{ a^T b}{a^Ta}a
B=b−Pb=b−a(aTa)−1aTb=b−aTaaTba
类似的,如果是对于一个平面而言,则应该是减去前面已正交化的两个基底,即:
C
=
c
−
P
1
c
−
P
2
c
=
c
−
a
T
c
a
T
a
a
−
B
T
c
B
T
B
B
C= c-P_1c-P_2 c=c-\frac{ a^T c}{a^Ta}a -\frac{ B^T c}{B^TB}B
C=c−P1c−P2c=c−aTaaTca−BTBBTcB
可以验证 a T B = 0 , a T C = 0 , B T a = 0 , B T C = 0 a^TB=0,a^TC=0,B^Ta=0, B^TC=0 aTB=0,aTC=0,BTa=0,BTC=0,这说明我们这种正交化方法的正确性。
根据图示我们着重来看列3正交化的过程,由前文已知,c减去的两项实际是c在a和B上的投影,如果我们先将这个投影相加,那么结果就是蓝色平面中的红色对角线ON,相当于要说明:
列
3
−
红
色
对
交
线
=
正
交
化
的
列
3
?
⇒
M
N
⊥
平
面
O
P
N
Q
?
列3 - 红色对交线 = 正交化的列3 ? \Rightarrow MN \perp 平面OPNQ ?
列3−红色对交线=正交化的列3?⇒MN⊥平面OPNQ?
实际上,这个结论是显然的:
O
P
⊥
P
M
,
O
P
⊥
P
N
(
O
Q
)
⇒
O
P
⊥
M
N
O
Q
⊥
Q
M
,
O
Q
⊥
Q
N
(
O
P
)
⇒
O
Q
⊥
M
N
M
N
⊥
O
P
,
M
N
⊥
O
P
⇒
M
N
⊥
平
面
O
P
N
Q
OP \perp PM,OP \perp PN(OQ) \Rightarrow OP \perp MN \\ OQ \perp QM,OQ \perp QN(OP) \Rightarrow OQ \perp MN \\ MN \perp OP , MN \perp OP \Rightarrow MN \perp 平面OPNQ \\
OP⊥PM,OP⊥PN(OQ)⇒OP⊥MNOQ⊥QM,OQ⊥QN(OP)⇒OQ⊥MNMN⊥OP,MN⊥OP⇒MN⊥平面OPNQ
由此可合理外推,对于超平面而言,第n个向量应减去已正交化的n-1个基底。当然,这种基于空间垂直出发的视角对于矩阵来说显得繁琐,同矩阵的
A
=
L
U
A=LU
A=LU变换一样,矩阵的标准正交化也对应着一个矩阵R,即
A
=
Q
R
A=QR
A=QR。我们容易看出这个R是一个上三角矩阵(验证矩阵Q的标准正交性质),以三维空间为例
[
a
b
c
]
=
[
q
1
q
2
q
3
]
[
q
1
T
a
q
1
T
b
q
1
T
c
0
q
2
T
b
q
2
T
c
0
0
q
3
T
c
]
\begin{bmatrix} a & b &c \end{bmatrix}= \begin{bmatrix} q_1 & q_2 & q_3 \end{bmatrix}\begin{bmatrix} q_1^Ta & q_1^Tb & q_1^Tc \\ 0 & q_2^Tb & q_2^Tc \\ 0 & 0 &q_3^Tc \end{bmatrix}
[abc]=[q1q2q3]⎣⎡q1Ta00q1Tbq2Tb0q1Tcq2Tcq3Tc⎦⎤
由A (=a) 、B、C推得
q
1
q_1
q1、
q
2
q_2
q2、
q
3
q_3
q3 的过程还差一步,即标准化,只需将已正交化的向量除以其模即可:
q
1
=
A
∣
A
∣
q
2
=
B
∣
B
∣
q
3
=
C
∣
C
∣
q _1 = \frac{A}{|A|} \ \ \ q _2 = \frac{B}{|B|} \ \ \ q _3 = \frac{C}{|C|}
q1=∣A∣A q2=∣B∣B q3=∣C∣C
例子
用一个例子作为结束:
A
=
[
1
3
1
2
2
2
3
1
1
]
=
[
a
1
a
2
a
3
]
A = \begin{bmatrix} 1 &3 & 1 \\ 2 & 2 &2 \\ 3 & 1 & 1 \end{bmatrix} =\begin{bmatrix} a_1 & a_2 & a_3 \end{bmatrix}
A=⎣⎡123321121⎦⎤=[a1a2a3]
先来看个空间直观一下:
% 画向量图
quiver3(0,0,0,1,2,3,'m'); hold on; quiver3(0,0,0,3,2,1,'black'); % 基底
hold on; quiver3(0,0,0,1,2,1,'r'); % 列3
V1 = [1;2;3]; V2 = [3;2;1];
% 求法向量
Vn = cross(V1,V2);
Vn = Vn/norm(Vn);
% 画单位法向量
hold on;quiver3(0,0,0,Vn(1),Vn(2),Vn(3),'g');
% 画平面
syms x1;syms x2;syms x3;
plane = -(x1*Vn(1)+x2*Vn(2))/Vn(3);
hold on;
fmesh(plane);
% 标注
legend('列1','列2','列3','单位法向量','列1和列2所在平面');
首先进行正交化,即
A
=
a
1
=
[
1
2
3
]
B
=
a
2
−
A
T
a
2
A
T
A
A
=
[
3
2
1
]
−
10
14
[
1
2
3
]
=
[
16
/
7
4
/
7
−
8
/
7
]
C
=
a
3
−
A
T
a
3
A
T
A
A
−
B
T
a
3
B
T
B
B
=
[
1
2
1
]
−
8
14
[
1
2
3
]
−
112
336
[
16
/
7
4
/
7
−
8
/
7
]
=
[
−
1
/
3
2
/
3
−
1
/
3
]
A = a_1= \begin{bmatrix} 1 \\2 \\ 3 \end{bmatrix} \\ B = a_2 - \frac{ A^T a_2}{A^TA}A =\begin{bmatrix} 3 \\2 \\ 1 \end{bmatrix}-\frac{10}{14}\begin{bmatrix} 1 \\2 \\ 3 \end{bmatrix}=\begin{bmatrix} 16/7\\ 4/7 \\ -8/7 \end{bmatrix} \\ C = a_3-\frac{ A^T a_3}{A^TA}A-\frac{ B^T a_3}{B^TB}B=\begin{bmatrix} 1 \\2 \\ 1 \end{bmatrix}-\frac{8}{14}\begin{bmatrix} 1 \\2 \\ 3 \end{bmatrix} -\frac{112}{336}\begin{bmatrix} 16/7\\ 4/7 \\ -8/7 \end{bmatrix} = \begin{bmatrix}-1/3 \\ 2/3 \\-1/3\end{bmatrix}
A=a1=⎣⎡123⎦⎤B=a2−ATAATa2A=⎣⎡321⎦⎤−1410⎣⎡123⎦⎤=⎣⎡16/74/7−8/7⎦⎤C=a3−ATAATa3A−BTBBTa3B=⎣⎡121⎦⎤−148⎣⎡123⎦⎤−336112⎣⎡16/74/7−8/7⎦⎤=⎣⎡−1/32/3−1/3⎦⎤
然后进行标准化,即
q
1
=
A
∣
A
∣
=
[
1
2
3
]
/
14
=
[
0.2673
0.5345
0.8018
]
q
2
=
B
∣
B
∣
=
[
16
/
7
4
/
7
−
8
/
7
]
/
336
/
49
=
[
0.8729
0.2182
−
0.4362
]
q
3
=
C
∣
C
∣
=
[
−
1
/
3
2
/
3
−
1
/
3
]
/
6
/
9
=
[
−
0.4082
0.8165
−
0.4082
]
q _1 = \frac{A}{|A|} = \begin{bmatrix} 1 \\2 \\ 3 \end{bmatrix}/ \sqrt{14} =\begin{bmatrix} 0.2673 \\ 0.5345 \\ 0.8018 \end{bmatrix} \\ q _2 = \frac{B}{|B|} = \begin{bmatrix} 16/7\\ 4/7 \\ -8/7 \end{bmatrix} / \sqrt{336/49} = \begin{bmatrix} 0.8729 \\ 0.2182 \\ -0.4362 \end{bmatrix}\\ q _3 = \frac{C}{|C|} = \begin{bmatrix}-1/3 \\ 2/3 \\-1/3\end{bmatrix}/\sqrt{6/9}=\begin{bmatrix} -0.4082 \\ 0.8165 \\ -0.4082\end{bmatrix}
q1=∣A∣A=⎣⎡123⎦⎤/14=⎣⎡0.26730.53450.8018⎦⎤q2=∣B∣B=⎣⎡16/74/7−8/7⎦⎤/336/49=⎣⎡0.87290.2182−0.4362⎦⎤q3=∣C∣C=⎣⎡−1/32/3−1/3⎦⎤/6/9=⎣⎡−0.40820.8165−0.4082⎦⎤
上述施密特正交化方法的编程实现见下,matlab中也可自带的orth函数可直接求一组标准正交基,可以验证结果的正确性。
a = [1,3,1;2,2,2;3,1,1];
%% 施密特正交化方法
[m,n] = size(a);
if(m<n)
error('行小于列,无法计算,请转置后重新输入');
end
b=zeros(m,n);
%正交化
b(:,1)=a(:,1);
for i=2:n
for j=1:i-1
b(:,i)=b(:,i)-dot(a(:,i),b(:,j))/dot(b(:,j),b(:,j))*b(:,j);
end
b(:,i)=b(:,i)+a(:,i);
end
%单位化
for k=1:n
b(:,k)=b(:,k)/norm(b(:,k));
end
%% 直接用orth方法
result = orth(a);
b
result
b =
0.2673 0.8729 -0.4082
0.5345 0.2182 0.8165
0.8018 -0.4364 -0.4082
result =
-0.5494 -0.7071 -0.4451
-0.6295 -0.0000 0.7770
-0.5494 0.7071 -0.4451