关于为什么最近开了这个坑:
因为发现学校里教的线代毕竟是比较理论的东西,教的计算方法写成程序拿来实战百分之百是会被残酷的数据碾爆的。毕竟现实中的矩阵一多半是大而优雅或者小而丑陋或者大而丑陋,哪个都不是大一线代课上讲的一个小小的高斯消元就能解决的。要么时间爆炸,要么数值误差爆炸,至少有一个适合你。
另一方面,计算机图形学中,无论是物理模拟还是几何处理都不可避免地要与大规模的矩阵打交道,因此很有必要建立自己的数值线性代数库,学习矩阵计算的方法。
这是线性代数剩下的还需要学习的内容。 这个系列博客会用来记录数值线性代数中的一些理论与应用。
没错你已经发现了这不是讲线性代数的而是讲数值线性代数的
关于数值线代在网上并没有找到什么特别猛料的网课(也可能是我手段不足?),不过MIT的Strang老爷子新开的"数据分析中的矩阵方法"还是非常不错,就是太简单了,属于是讲给宝宝的矩阵计算入门特别适合下饭或者赶路。这一次我们参考数值线代的圣经 Matrix Computation 来学习。
2. 矩阵分析基础
很多东西在别的地方都有提过,无论是数学物理方法还是数值分析都有讲过向量范数,矩阵范数等等的概念,这些概念在此当然不会再全部扯一遍当然在那些地方我也不会再扯一遍。这里只是过一些简单的结论和话题。这一章的主要目的应当还是掌握必备的误差估计的知识。
2.1. First Things First
很多线性代数里老生常谈的东西这里也就不多加叙述了。但是有些东西,学院派的线性代数可能不怎么管它,在数值方法里却还是有点用处的。
Sherman-Morrison-Woodbury公式
对于可逆矩阵A,B来说,有B−1=A−1−B−1(B−A)A−1。验证起来还是很容易的。
没人会蠢到用这个东西去求逆,但是如果我们用A+UVT代替这里的B,就可以得到 Sherman-Morrison-Woodbury公式:
(A+UVT)−1=A−1−(A+UVT)−1UVTA−1=A−1−A−1(E+UVTA−1)−1UVTA−1
=A−1−A−1(U−1+U−1UVTA−1)−1VTA−1=A−1−A−1(U−1+VTA−1)−1VTA−1
=A−1−A−1U(Ek+VTA−1U)−1VTA−1
这里唯一要运用的就是矩阵乘积的逆的性质,并不难(但是我居然还想了一会儿,丢人)。
这个公式看起来更丑了,但是它好就好在U,V可以不是方阵,事实上U,V∈Rn×k就可以,而k可以是任意的(对于任意的k的情况可以由方阵情况推出)。特别的我们取k=1,这时U,V退化为向量u,v。取α=1+vTA−1u就有
(A+uvT)−1=A−1−α1A−1uvTA−1
这个式子被称为Sherman-Morrison公式,很有用,在Strang老爷子的课里也讲过这个东西,可以用于对矩阵的低秩扰动的求逆,后者就更加有用了。
矩阵和向量范数
向量范数再稍微提一个性质也就是著名的赫尔德不等式:
∣xTy∣≤∥x∥p∥y∥q,p1+q1=1
首先说明一点,在有限维空间中,范数都是等价的,因此这里我们取各种乱七八糟的范数,数值效果可能有区别,不过至少收敛性上都是没啥问题的。但是在无穷维空间上就没那么幸运了,这些事情就需要留给犯寒泛函分析,这里当然就不提了。
Frobenius范数自不必多说,这里再简单回顾一下矩阵的算子范数。
矩阵A的算子范数定义为∥A∥=∥x∥=1max∥Ax∥,当然,这是和选取的向量范数有关的。
关于算子范数有一个挺不错的性质:
∥AB∥≤∥A∥∥B∥
这对任意的矩阵A,B都成立,当然实际上考虑到矩阵范数作为函数,对于不同形状的矩阵是不同的,这个不等式比看起来的还是要更好,它甚至能揭示出三个完全不一样的函数的关系。
还有另外一种矩阵范数的定义,更宽泛:
∥A∥α,β=∥x∥α=1max∥Ax∥β
矩阵的2-范数
Thm. 存在单位向量x使得ATAz=μ2x,μ=∥A∥2
这个定理的证明还是得依据定义来。。。假设存在这样的z就一定有∥Az∥2=∥A∥2
令g(x)=2∥x∥22∥Ax∥22
z最小化g(x)于是就可以求梯度并且令其为0,最终能得到
ATAz=(zTATAz)z
比较一下就知道∥Az∥2=μ
这其实揭示了一个事情:∥A∥22是p(λ)=det(ATA−λI)的零点,也就是说
∥A∥2=λmax(ATA)
逆矩阵的一些误差估计
Lemma. 设F∈Rn×n,∥F∥p<1,则I−F可逆,且有
∥(I−F)−1∥p≤1−∥F∥1
可逆性的证明是经典的习题,因为我们能发现I−F是一个对角占优矩阵,对角占优矩阵是可逆的。不过这里还有另一个结论:
(I−F)−1=k=0∑∞Fk
这么写当然是OK的,因为F的范数小于1,这一定是收敛的。
两侧取范数,右边变成无穷等比数列求和,就得到了范数不等式。
这个引理表明了另一个事情:
∥(I−F)−1−I∥p≤1−∥F∥p∥F∥p
考虑∥F∥p∼O(ϵ)且ϵ≪1,这意味着,如果我们对单位矩阵进行O(ϵ)的扰动,其逆的变化也是O(ϵ)的。
有了这个引理就能得到下面的定理:
Thm. 若A可逆,且r=∥A−1E∥p<1,那么A+E可逆且
(A+E)−1−A−1≤1−r∥E∥p∥A−1∥p2
Proof. 设F=−EA−1,容易验证
(A+E)−1−A−1=−A−1EA−1(I+F)−1
两侧取范数即可
The Eckhart-Young Theorem
这个定理的前置是SVD,SVD到处都有讲,本来应该说是个很重要的概念,但是因为实在是过于普遍所以这里就不多提了。
所幸网上已有很多资料,读者可以自行学习。这里至少应该给出个结论但是我连结论也不会给的。
这个定理在Strang老爷子的课里也专门讲过,但是Strang没有讲证明实际上他那个课都不怎么讲证明。这个定理是矩阵低秩近似的理论基础,还是相当重要的。
Thm. 若k<r=rk(A),令Ak=i=1∑kσiuiviT,则有
rk(B)=kmin∥A−B∥2=∥A−Ak∥2=σk+1
当然,这里σk+1说的是奇异值,并且σ是从大往小排列的。
第二个等式相当容易验证的,重点是第一个。
首先很容易发现Ak的确是个秩为k的矩阵,我们要证明没有秩为k矩阵能比Ak更懂如何逼近逼近得更好,即要证明
∥A−B∥22≥σk+12
首先寻找矩阵B的零空间null(B)的一组正交基{x1,x2,...xn−k}
运用维数公式,我们有
dim(null(B))+dim(<v1,v2,...,vk+1>)−dim(null(B)∩<v1,v2,...,vk+1>)≤n
即
n−k+k+1−dim(null(B)∩<v1,v2,...,vk+1>)≤n
即有dim(null(B)∩<v1,v2,...,vk+1>)≥1
这说明S=null(B)∩<v1,v2,...,vk+1>={0},线性代数老师应该也提醒过这个地方绝对不能把{0}写成空集。
设z是S中的单位向量,于是有
Az=i=1∑k+1uiσiviTz=i=1∑k+1σiui(zTvi)=i=1∑k+1σi(viTz)ui
∥Az−Bz∥22=∥Az∥22=∥i=1∑k+1σi(viTz)ui∥22
首先由于ui是单位正交基,我们可以直接把范数和求和交换并且化简。于是
∥(A−B)z∥22=i=1∑k+1σi2(viTz)2
其次,由于z∈<v1,v2,...,vk+1>并且∥z∥22=1,而且v1,v2,...,vk+1都是单位正交基,有
∥z∥22=i=1∑k+1(viTz)2=1
所以说∥(A−B)z∥22=i=1∑k+1σi2(viTz)2是σ12,σ22,...σk+12的凸组合,它肯定比最小的那一个即σk+12要大,由于对任意满足条件的∥(A−B)z∥22都成立,我们取其上确界应当还是成立的,由于这个上确界一定不比∥A−B∥22大(因为z取不到全空间)。就有
∥A−B∥22≥σk+12
2.2. 子空间度量
很多时候我们希望能够衡量两个子空间有多“接近”,或者说两个子空间有多少“重合”。这个时候考虑子空间的正交基就是非常好的,因为正交基拥有很好的运算性质。
2.2.1. 正交投影
其实要说吧,这个概念线代课上也讲过,不过这里还是重新提一下,至少给出一些重要结论,尤其是一些矩阵结论坦白说这一块儿我是没学好的最后还是滚去问了线代老师。
正交投影表示把一个向量投影到子空间上,所谓投影也就是说将向量分解为在子空间和其正交补空间中的两部分然后取前者。正交投影其作为算子肯定有其对应的方阵,依此定义:
如果P∈Rn×n满足ran(P)=S,且P2=P,PT=P称其为子空间S上的正交投影。
这里ran就是把矩阵作为线性映射的像空间,或者说矩阵的列空间。
ran(P)是列空间。本来投影算子的定义其实可以没有对称性的,但是这里有。这一点。。。挺迷的
一个重要的结论是,如果子空间S的一组单位正交基拼成的矩阵是V=[v1∣v2∣...∣vk],那么P上的正交投影矩阵P=VVT。
这个玩意儿。。。很容易验证它确实满足上面的定义。但是这个结构我不是很能用线性映射去理解它。。。
在此特别感谢我的线性代数老师,拜托他解读了一下其中的线性映射理解,写得非常清晰。
有了这些东西,再结合SVD,我们就能够对子空间进行各种拆解。
考虑奇异值分解A=UΣVT,设rank(A)=k,将U写为[Uk∣Uk~],其中Uk=[u1,u2,...,uk],用同样的方法将V写为[Vk∣Vk~]
于是我们就知道,ran(A)=ran(UΣVT)=ran(ΣVT)=ran(VΣT)=ran(Vk)。第二个等号成立是因为U是满秩矩阵。
这说明Vk是ran(A)的一组单位正交基,由于由上面的结论我们就可以把各种正交投影矩阵写出来。
例如,假设做奇异值分解后A=UΣVT,设rank(A)=k,且U=[Uk~∣Un−k~],V=[Vk~∣Vn−k~],我们就能根据上面的结论写出null(A),ran(A)等等的正交投影矩阵。
2.2.2. 子空间的距离
其实距离好定义,就是两个子空间(维度相同)中点的距离的最小值。我们同样可以用矩阵来给出这个度量。
Thm.设W=[W1∣W2],Z=[Z1∣Z2]都是n×n方阵,W1,Z1∈Rn×k,若S1=ran(W1),S2=ran(Z1),则有
dist(S1,S2)=∥WTZ2∥2=∥Z1TW2∥2
证明主要运用分块矩阵,关键是发现dist(S1,S2)=∥W1W1T−Z1Z1T∥2
2.2.3. CS分解
一个简化版的结论是:
Thm. 考虑矩阵
Q=[Q1Q2]
且Q1∈Rm1×n,Q2∈Rm2×n
如果Q的列向量正交,则存在正交矩阵U1,U2,V使得
[U100U2][Q1Q2]V=[CS]
其中C,S都是对角矩阵,且C=diag(cos(θ1),...cos(θn)),S=C=diag(sin(θ1),...sin(θn))
证明方法是首先对Q1做SVD得到U1,V,然后构造U2,分离出Q2的一个可以同时被V对角化的子空间。证明涉及到很多量化的性质。
一般性的结论是:
Thm.若Q=[Q11Q21Q12Q22]是正交矩阵,就可以把它写成一个很丑的但是只含C,S,I,0块的矩阵。。。
如果去看看这个形式其实能发现两个子空间的正交基确实有一定联系但是不知道有啥用。。。
2.3. 线性系统的敏感性
实话说。。。上面那一节真挺无聊也不知道有啥用差点让我还没起步就弃坑了,但是这个话题就很有用了。
这一节我们的研究对象就只有一个Ax=b,A∈Rn×n,b∈Rn
2.3.1. 条件数
这其实是数值分析里讲过的话题。。。但是我很幸运地看完了那一章然后一点都不记得(
通过奇异值分解,我们可以写出来A=i=1∑nσiuiviT
于是解得x=i=1∑nσiuiTbvi
可以发现,如果不幸的是某些σi极其小,那么在这些方向上对A或者b一点小扰动就会让解产生很大的偏移。
为了更加精细地衡量这个偏移,我们考虑一个微扰的线性系统
(A+ϵF)x(ϵ)=b+ϵf,x(0)=x
假设A是可逆的,那么x(ϵ)在原点至少一个小邻域内是可积的,而且ϵ很小,我们可以在原点附近用泰勒展开来近似x(ϵ),但是要求x(ϵ)的导数直接算并不是很直观,所以我们用隐函数求导会比较方便。
Ax˙(ϵ)+Fx(ϵ)+ϵFx˙(ϵ)=f
代入ϵ=0得x˙(ϵ)=A−1(f−Fx)
现在再用泰勒展开
x(ϵ)−x=A−1(f−Fx)ϵ+O(ϵ2)
为了估计误差,我们将其取范数并在两侧除以∥x∥,得到
∥x∥∥x(ϵ)−x∥≤ϵ∥A−1∥{∥x∥∥f∥+∥F∥}+O(ϵ2)
对于方阵A,我们定义A的条件数为κ(A)=∥A∥∥A−1∥
如果A不可逆,我们可以顺利成章地定义κ(A)=+∞
由于∥b∥∥A∥−1≤∥x∥,我们可以把上面的不等式里的∥x∥放缩掉,来把右边整理成和A,b的相对误差有关的不等式
∥x∥∥x(ϵ)−x∥≤κ(A)(ρA+ρb)+O(ϵ2)
其中ρA=ϵ∥A∥∥F∥,ρA=ϵ∥b∥∥f∥
我们就发现,x的相对误差差不多是A,b的相对误差的κ(A)倍。于是κ(A)就衡量了Ax=b的敏感性。根据前面的范数定义可以写出2-范数下定义的条件数和矩阵奇异值的关系为
κ(A)=σminσmax
这上面的条件数都是在2-范数下定义的,实际上条件数的定义不依赖于范数,定义也有很多,这里不全列出。但是要记住我们是通过微分引入条件数的。
条件数是衡量矩阵的数值性质的非常好的量,很自然地会想到行列式是不是也是这么样一个量,但是事实上矩阵的行列式和线性系统的解的数值性质关系并不大。
最后,有了条件数,我们就可以用它来较精确地估计线性系统解的偏移量。
Lemma. 设方程Ax=b,(A+ΔA)y=b+Δb,且∥ΔA∥≤ϵ∥A∥,∥Δb∥≤ϵ∥b∥,则有
∥x∥∥y∥≤1−ϵκ(A)1+ϵκ(A)
对微扰后的方程左乘A−1然后运用条件数的定义和范数性质放缩可得。
Thm. 如果引理中条件满足,那么
∥x∥∥y−x∥≤1−ϵκ(A)2ϵκ(A)
证明不难。这个界还行,但还可以更加精细化:
∥x∥∞∥y−x∥∞≤1−ϵκ(A)2ϵ∥∣A−1∣∣A∣∥∞
其中∣A∣表示把|A|的元素全部取绝对值。∥∣A−1∣∣A∣∥∞被称为Skeel条件数。
2.4. 有限精度的计算
这个很重要,因为计算机里用的浮点数都是精度有限的。
IEEE标准这种人尽皆知的东西就无需多言了。
2.4.1. 舍入误差
在IEEE标准里,用一定位数的二进制位编码指数,这个位数决定了表示的最低精度,更高的精度就必须舍入,通常采用的是舍入到最近浮点数,也就是所谓的“四舍六入五成双”。
用round(x)表示x舍入后的结果,容易发现,如果我们使用e位编码指数,那么就有
∣x∣∣round(x)−x∣≤2−e−1
如果考虑到这个层次的细节,那也太过于底层了。计算机按照IEEE标准实现的浮点运算还是值得相信的,为了方便这个层面上的表示,用fl(x)表示x作为一个浮点数的存储或者计算。
于是我们就知道fl(x)=x(1+δ),∣δ∣<2−e−1
2−e−1在书中被记作u,对单精度浮点来说差不多是10−7,对双精度浮点来说差不多是10−16,可以发现如果把这两种浮点数用的编码位带进去,的确是差不多的数量级。
2.4.2. 浮点数注意事项
计算顺序很重要
这是计算机内容,神书csapp强调过这个事情,无需多言。
逼近精度小不一定是真的小
例如,如果我们用一阶差分算sinx的数值导数,取步长h=u能得到近似最小的误差。
更细节地说一下这个事情,我们如果试图用一阶差分也就是
d=hsin(x+h)−sinx来近似计算导数。假设sin的计算没有误差,那么由数值分析知识(其实就是一个泰勒展开),我们有$|d - \cos x|\sim O(h) $
事实上,正弦函数的数值计算本身也有误差,这是我们无法改变的,我们除以h的操作把这个误差放大了h1倍,这意味着h绝非取得越小越好,因此实际上我们取h=u才能近似地最小化此误差。
2.4.3. 应用:存储矩阵
很容易就能知道,如果把一个矩阵A存在计算机里,会有
∣fl(A)−A∣≤u∣A∣
∥fl(A)−A∥1≤u∥A∥1
这边很多结论就不细说了,基本上都是一个放缩就能放出来的事情。大部分结论这里会跳过,后面需要用到的时候直接给出结果,不给出证明。
但是有个结论还是在这里提一嘴。
Lemma. 如果(1+α)=k=1∏n(1+αk)且∣αk∣≤u且nu≤0.01,那么∣α∣≤1.01nu
这个结果还是比较恐怖的,它表明如果两个向量的夹角比较接近垂直的话,向量点积的相对误差可能会很大。
2.4.4. 前向和后向误差分析
舍入误差拿来计算相对误差时,既可以是和实际计算结果作比,也可以是和输入数据作比。前者就是前向误差分析,后者就是后向误差分析。
原书在这里还分析了一下Strassen矩阵乘法,略去。
2.4.5. 分析理想的解方程组
我们重新回到数值线性代数的一个重要课题(也是下一个课题):解方程组Ax=b.
假设我们解出来的解是x^,那么它满足
(A+E)x^=b+e
这里的E不是单位矩阵,而是误差矩阵。有∥E∥∞≤u∥A∥∞,∥e∥∞≤u∥b∥∞ 并且fl(b)=b+e,fl(A)=A+E
假设uκ∞(A)≤21,那么利用前面在矩阵条件数那里得到的结果,我们就知道
∥x∥∞∥x−x^∥∞≤4uκ∞(A)
Fin.