强撸MIT18-06灰飞烟灭(二)
第十九讲:行列式公式和代数余子式
上一讲中,我们从三个简单的性质扩展出了一些很好的推论,本讲将继续使用这三条基本性质:
- detI=1detI=1;
- 交换行行列式变号;
- 对行列式的每一行都可以单独使用线性运算,其值不变;
我们使用这三条性质推导二阶方阵行列式:
按照这个方法,我们继续计算三阶方阵的行列式,可以想到,我们保持第二、三行不变,将第一行拆分为个行列式之和,再将每一部分的第二行拆分为三部分,这样就得到九个行列式,再接着拆分这九个行列式的第三行,最终得到二十七个行列式。可以想象到,这些矩阵中有很多值为零的行列式,我们只需要找到不为零的行列式,求和即可。
同理,我们想继续推导出阶数更高的式子,按照上面的式子可知n阶行列式应该可以分解成n!个非零行列式(占据第一行的元素有n种选择,占据第二行的元素有n−1种选择,以此类推得n!):
这个公式还不完全,接下来需要考虑如何确定符号:
- 观察带有下划线的元素,它们的排列是(4,3,2,1),变为(1,2,3,4)需要两步操作,所以应取+;
- 观察带有上划线的元素,它们的排列是(3,2,1,4),变为(1,2,3,4)需要一步操作,所以应取−。
- 观察其他元素,我们无法找出除了上面两种以外的排列方式,于是该行列式值为零,这是一个奇异矩阵。
此处引入代数余子式(cofactor)的概念,它的作用是把n阶行列式化简为n−1阶行列式。
于是我们把(1)式改写为:
于是,我们可以定义aij的代数余子式:将原行列式的第i行与第j列抹去后得到的n−1阶行列式记为Cij,i+j为偶时时取+,i+j为奇时取−。
现在再来完善式子(2):将行列式A沿第一行展开:
到现在为止,我们了解了三种求行列式的方法:
- 消元,detA就是主元的乘积;
- 使用(2)式展开,求n!项之积;
- 使用代数余子式。
计算例题:
A4=|1100111001110011|沿第一行展开=|110111011|−|110011011|=−1−0=−1
第二十讲:克拉默法则、逆矩阵、体积
本讲主要介绍逆矩阵的应用。
求逆矩阵
我们从逆矩阵开始,对于二阶矩阵有[abcd]−1=1ad−bc[d−b−ca]。观察易得,系数项就是行列式的倒数,而矩阵则是由一系列代数余子式组成的。先给出公式:
观察这个公式是如何运作的,化简公式得ACT=(detA)I,写成矩阵形式有[a11a12⋯a1n⋮⋮⋱⋮an1an2⋯ann][C11⋯Cn1C12⋯Cn2⋮⋱⋮C1n⋯Cnn]=Res
对于这两个矩阵的乘积,观察其结果的元素Res11=a11C11+a12C12+⋯+a1nC1n,这正是上一讲提到的将行列式按第一行展开的结果。同理,对Res22,⋯,Resnn都有Resii=detA,即对角线元素均为detA。
再来看非对角线元素:回顾二阶的情况,如果用第一行乘以第二行的代数余子式a11C21+a12C22,得到a(−b)+ab=0。换一种角度看问题,a(−b)+ab=0也是一个矩阵的行列式值,即As=[abab]。将detAs按第二行展开,也会得到detAs=a(−b)+ab,因为行列式有两行相等所以行列式值为零。
推广到n阶,我们来看元素Res1n=a11Cn1+a12Cn2+⋯+a1nCnn,该元素是第一行与最后一行的代数余子式相乘之积。这个式子也可以写成一个特殊矩阵的行列式,即矩阵As=[a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮an−a1an−12⋯an−1na11a12⋯a1n]。计算此矩阵的行列式,将detAs按最后一行展开,也得到detAs=a11Cn1+a12Cn2+⋯+a1nCnn。同理,行列式As有两行相等,其值为零。
结合对角线元素与非对角线元素的结果,我们得到Res=[detA0⋯00detA⋯0⋮⋮⋱⋮00⋯detA],也就是(1)等式右边的(detA)I,得证。
求解Ax=b
因为我们现在有了逆矩阵的计算公式,所以对Ax=b有x=A−1b=1detACTb,这就是计算x的公式,即克莱默法则(Cramer's rule)。
现在来观察x=1detACTb,我们将得到的解拆分开来,对x的第一个分量有x1=y1detA,这里y1是一个数字,其值为y1=b1C11+b2C21+⋯+bnCn1,每当我们看到数字与代数余子式乘之积求和时,都应该联想到求行列式,也就是说y1可以看做是一个矩阵的行列式,我们设这个矩阵为B1。所以有xi=detB1detA,同理有x2=detB2detA,x2=detB2detA。
而B1是一个型为[ba2a3⋯an]的矩阵,即将矩阵A的第一列变为b向量而得到的新矩阵。其实很容易看出,detB1可以沿第一列展开得到y1=b1C11+b2C21+⋯+bnCn1。
一般的,有Bj=[a1a2⋯aj−1baj+1⋯an],即将矩阵A的第j列变为b向量而得到的新矩阵。所以,对于解的分量有xj=detBjdetA。
这个公式虽然很漂亮,但是并不方便计算。
关于体积(Volume)
先提出命题:行列式的绝对值等于一个箱子的体积。
来看三维空间中的情形,对于3阶方阵A,取第一行(a1,a2,a3),令其为三维空间中点A1的坐标,同理有点A2,A3。连接这三个点与原点可以得到三条边,使用这三条边展开得到一个平行六面体,‖detA‖就是该平行六面体的体积。
对于三阶单位矩阵,其体积为detI=1,此时这个箱子是一个单位立方体。这其实也证明了前面学过的行列式性质1。于是我们想,如果能接着证明性质2、3即可证明体积与行列式的关系。
对于行列式性质2,我们交换两行并不会改变箱子的大小,同时行列式的绝对值也没有改变,得证。
现在我们取矩阵A=Q,而Q是一个标准正交矩阵,此时这个箱子是一个立方体,可以看出其实这个箱子就是刚才的单位立方体经过旋转得到的。对于标准正交矩阵,有QTQ=I,等式两边取行列式得det(QTQ)=1=|QT||Q|,而根据行列式性质10有|QT|=|Q|,所以原式=|Q|2=1,|Q|=±1。
接下来在考虑不再是“单位”的立方体,即长方体。 假设Q矩阵的第一行翻倍得到新矩阵Q2,此时箱子变为在第一行方向上增加一倍的长方体箱子,也就是两个“标准正交箱子”在第一行方向上的堆叠。易知这个长方体箱子是原来体积的两倍,而根据行列式性质3.a有detQ2=detQ,于是体积也符合行列式的数乘性质。
我们来看二阶方阵的情形,|a+a′b+b′cd|=|abcd|+|a′b′cd|。在二阶情况中,行列式就是一个求平行四边形面积的公式,原来我们求由四个点(0,0),(a,b),(c,d),(a+c,b+d)围成的四边形的面积,需要先求四边形的底边长,再做高求解,现在只需要计算detA=ad−bc即可(更加常用的是求由(0,0),(a,b),(c,d)围成的三角形的面积,即12ad−bc)。也就是说,如果知道了歪箱子的顶点坐标,求面积(二阶情形)或体积(三阶情形)时,我们不再需要开方、求角度,只需要计算行列式的值就行了。
再多说两句我们通过好几讲得到的这个公式,在一般情形下,由点(x1,y1),(x2,y2),(x3,y3)围成的三角形面积等于12|x1y11x2y21x3y31|,计算时分别用第二行、第三行减去第一行化简到第三列只有一个1(这个操作实际作用是将三角形移动到原点),得到12|x1y11x2−x1y2−y10x3−x1y3−y10|,再按照第三列展开,得到三角形面积等于(x2−x1)(y3−y1)−(x3−x1)(y2−y1)2。
第二十一讲:特征值和特征向量
特征值、特征向量的由来
给定矩阵A,矩阵A乘以向量x,就像是使用矩阵A作用在向量x上,最后得到新的向量Ax。在这里,矩阵A就像是一个函数,接受一个向量x作为输入,给出向量Ax作为输出。
在这一过程中,我们对一些特殊的向量很感兴趣,他们在输入(x)输出(Ax)的过程中始终保持同一个方向,这是比较特殊的,因为在大多情况下,Ax与x指向不同的方向。在这种特殊的情况下,Ax平行于x,我们把满足这个条件的x成为特征向量(Eigen vector)。这个平行条件用方程表示就是:
-
对这个式子,我们试着计算特征值为0的特征向量,此时有Ax=0,也就是特征值为0的特征向量应该位于A的零空间中。
也就是说,如果矩阵是奇异的,那么它将有一个特征值为λ=0。
-
我们再来看投影矩阵P=A(ATA)−1AT的特征值和特征向量。用向量b乘以投影矩阵P得到投影向量Pb,在这个过程中,只有当b已经处于投影平面(即A的列空间)中时,Pb与b才是同向的,此时b投影前后不变(Pb=1⋅b)。
即在投影平面中的所有向量都是投影矩阵的特征向量,而他们的特征值均为1。
再来观察投影平面的法向量,也就是投影一讲中的e向量。我们知道对于投影,因为e⊥C(A),所以Pe=0e,即特征向量e的特征值为0。
于是,投影矩阵的特征值为λ=1,0。
-
再多讲一个例子,二阶置换矩阵A=[0110],经过这个矩阵处理的向量,其元素会互相交换。
那么特征值为1的特征向量(即经过矩阵交换元素前后仍然不变)应该型为[11]。
特征值为−1的特征向量(即经过矩阵交换元素前后方向相反)应该型为[1−1]。
再提前透露一个特征值的性质:对于一个n×n的矩阵,将会有n个特征值,而这些特征值的和与该矩阵对角线元素的和相同,因此我们把矩阵对角线元素称为矩阵的迹(trace)。n∑i=1λi=n∑i=1aii
在上面二阶转置矩阵的例子中,如果我们求得了一个特征值1,那么利用迹的性质,我们就可以直接推出另一个特征值是−1。
求解Ax=λx
对于方程Ax=λx,有两个未知数,我们需要利用一些技巧从这一个方程中一次解出两个未知数,先移项得(A−λI)x=0。
观察(A−λI)x=0,右边的矩阵相当于将A矩阵平移了λ个单位,而如果方程有解,则这个平移后的矩阵(A−λI)一定是奇异矩阵。根据前面学到的行列式的性质,则有det(A−λI)=0
这样一来,方程中就没有x了,这个方程也叫作特征方程(characteristic equation)。有了特征值,代回(A−λI)x=0,继续求(A−λI)的零空间即可。
-
现在计算一个简单的例子,A=[3113],再来说一点题外话,这是一个对称矩阵,我们将得到实特征值,前面还有置换矩阵、投影矩阵,矩阵越特殊,则我们得到的特征值与特征向量也越特殊。看置换矩阵中的特征值,两个实数1,−1,而且它们的特征向量是正交的。
回到例题,计算det(A−λI)=|3−λ113−λ|,也就是对角矩阵平移再取行列式。原式继续化简得(3−λ)2−1=λ2−6λ+8=0,λ1=4,λ2=2。可以看到一次项系数−6与矩阵的迹有关,常数项与矩阵的行列式有关。
继续计算特征向量,A−4I=[−111−1],显然矩阵是奇异的(如果是非奇异说明特征值计算有误),解出矩阵的零空间x1=[11];同理计算另一个特征向量,A−2I=[1111],解出矩阵的零空间x2=[1−1]。
回顾前面转置矩阵的例子,对矩阵A′=[0110]有λ1=1,x1=[11],λ2=−1,x2=[−11]。看转置矩阵A′与本例中的对称矩阵A有什么联系。
易得A=A′+3I,两个矩阵特征值相同,而其特征值刚好相差3。也就是如果给一个矩阵加上3I,则它的特征值会加3,而特征向量不变。这也很容易证明,如果Ax=λx,则(A+3I)x=λx+3x=(λ+3)x,所以x还是原来的x,而λ变为λ+3。
接下来,看一个关于特征向量认识的误区:已知Ax=λx,Bx=αx,则有(A+B)x=(λ+α)x,当B=3I时,在上例中我们看到,确实成立,但是如果B为任意矩阵,则推论不成立,因为这两个式子中的特征向量x并不一定相同,所以两个式子的通常情况是Ax=λx,By=αy,它们也就无从相加了。
-
再来看旋转矩阵的例子,旋转90∘的矩阵Q=[cos90−sin90sin90cos90]=[0−110](将每个向量旋转90∘,用Q表示因为旋转矩阵是正交矩阵中很重要的例子)。
上面提到特征值的一个性质:特征值之和等于矩阵的迹;现在有另一个性质:特征值之积等于矩阵的行列式。n∏i=1λi=detA
对于Q矩阵,有{λ1+λ2=0λ1⋅λ2=1,再来思考特征值与特征向量的由来,哪些向量旋转90∘后与自己平行,于是遇到了麻烦,并没有这种向量,也没有这样的特征值来满足前面的方程组。
我们来按部就班的计算,det(Q−λI)=|λ−11λ|=λ2+1=0,于是特征值为λ1=i,λ2=−i,我们看到这两个值满足迹与行列式的方程组,即使矩阵全是实数,其特征值也可能不是实数。本例中即出现了一对共轭负数,我们可以说,如果矩阵越接近对称,那么特征值就是实数。如果矩阵越不对称,就像本例,QT=−Q,这是一个反对称的矩阵,于是我得到了纯虚的特征值,这是极端情况,通常我们见到的矩阵是介于对称与反对称之间的。
于是我们看到,对于好的矩阵(置换矩阵)有实特征值及正交的特征向量,对于不好的矩阵(90∘旋转矩阵)有纯虚的特征值。
-
再来看一个更糟的情况,A=[3103],这是一个三角矩阵,我们可以直接得出其特征值,即对角线元素。来看如何得到这一结论的:det(A−λI)=|3−λ103−λ|=(3−λ)2=0,于是λ1=3,λ2=3。而我们说这是一个糟糕的状况,在于它的特征向量。
带入特征值计算特征向量,带入λ1=3得(A−λI)x=[0100][x1x2]=[00],算出一个特征值x1=[10],当我们带入第二个特征值λ1=3时,我们无法得到另一个与x1线性无关的特征向量了。
而本例中的矩阵A是一个退化矩阵(degenerate matrix),重复的特征值在特殊情况下可能导致特征向量的短缺。
这一讲我们看到了足够多的“不好”的矩阵,下一讲会介绍一般情况下的特征值与特征向量。
第二十二讲:对角化和A的幂
对角化矩阵
上一讲我们提到关键方程Ax=λx,通过det(A−λI)=0得到特征向量λ,再带回关键方程算出特征向量x。
在得到特征值与特征向量后,该如何使用它们?我们可以利用特征向量来对角化给定矩阵。
有矩阵A,它的特征向量为x1,x2,⋯,xn,使用特征向量作为列向量组成一个矩阵S=[x1x2⋯xn],即特征向量矩阵, 再使用公式S−1AS=Λ将A对角化。注意到公式中有S−1,也就是说特征向量矩阵S必须是可逆的,于是我们需要n个线性无关的特征向量。
现在,假设A有n个线性无关的特征向量,将它们按列组成特征向量矩阵S,则AS=A[x1x2⋯xn],当我们分开做矩阵与每一列相乘的运算时,易看出Ax1就是矩阵与自己的特征向量相乘,其结果应该等于λ1x1。那么AS=[(λ1x1)(λ2x2)⋯(λnxn)]。可以进一步化简原式,使用右乘向量按列操作矩阵的方法,将特征值从矩阵中提出来,得到[x1x2⋯xn][λ10⋯00λ2⋯0⋮⋮⋱⋮00⋯λn]=SΛ。
于是我们看到,从AS出发,得到了SΛ,特征向量矩阵又一次出现了,后面接着的是一个对角矩阵,即特征值矩阵。这样,再继续左乘S−1就得到了公式(1)。当然,所以运算的前提条件是特征向量矩阵S可逆,即矩阵A有n个线性无关的特征向量。这个式子还要另一种写法,A=SΛS−1。
我们来看如何应用这个公式,比如说要计算A2。
- 先从Ax=λx开始,如果两边同乘以A,有A2x=λAx=λ2x,于是得出结论,对于矩阵A2,其特征值也会取平方,而特征向量不变。
- 再从A=SΛS−1开始推导,则有A2=SΛS−1SΛS−1=SΛ2S−1。同样得到特征值取平方,特征向量不变。
两种方法描述的是同一个现象,即对于矩阵幂运算A2,其特征向量不变,而特征值做同样的幂运算。对角矩阵Λ2=[λ210⋯00λ22⋯0⋮⋮⋱⋮00⋯λ2n]。
特征值和特征向量给我们了一个深入理解矩阵幂运算的方法,Ak=SΛkS−1。
再来看一个矩阵幂运算的应用:如果k→∞,则Ak→0(趋于稳定)的条件是什么?从SΛkS−1易得,|λi|<1。再次强调,所有运算的前提是矩阵A存在n个线性无关的特征向量。如果没有n个线性无关的特征向量,则矩阵就不能对角化。
关于矩阵可对角化的条件:
-
如果一个矩阵有n个互不相同的特征值(即没有重复的特征值),则该矩阵具有n个线性无关的特征向量,因此该矩阵可对角化。
-
如果一个矩阵的特征值存在重复值,则该矩阵可能具有n个线性无关的特征向量。比如取10阶单位矩阵,I10具有10个相同的特征值1,但是单位矩阵的特征向量并不短缺,每个向量都可以作为单位矩阵的特征向量,我们很容易得到10个线性无关的特征向量。当然这里例子中的I10的本来就是对角矩阵,它的特征值直接写在矩阵中,即对角线元素。
同样的,如果是三角矩阵,特征值也写在对角线上,但是这种情况我们可能会遇到麻烦。矩阵A=[2102],计算行列式值det(A−λI)=|2−λ102−λ|=(2−λ)2=0,所以特征值为λ1=λ2=2,带回Ax=λx得到计算[0100]的零空间,我们发现x1=x2=[10],代数重度(algebraic multiplicity,计算特征值重复次数时,就用代数重度,就是它作为多项式根的次数,这里的多项式就是(2−λ)2)为2,这个矩阵无法对角化。这就是上一讲的退化矩阵。
我们不打算深入研究有重复特征值的情形。
求uk+1=Auk
从u1=Au0开始,u2=A2u0,所有uk=Aku0。下一讲涉及微分方程(differential equation),会有求导的内容,本讲先引入简单的差分方程(difference equation)。本例是一个一阶差分方程组(first order system)。
要解此方程,需要将u0展开为矩阵A特征向量的线性组合,即u0=c1x1+c2x2+⋯+cnxn=[x1x2⋯xn][c1c2⋮cn]=Sc。于是Au0=c1Ax1+c2Ax2+⋯+cnAxn=c1λ1x1+c2λ2x2+⋯+cnλnxn。继续化简原式,Au0=[x1x2⋯xn][λ10⋯00λ2⋯0⋮⋮⋱⋮00⋯λn][c1c2⋮cn]=SΛc。用矩阵的方式同样可以得到该式:Au0=SΛS−1u0=SΛS−1Sc=SΛc。
那么如果我们要求A100u0,则只需要将λ变为λ100,而系数c与特征向量x均不变。
当我们真的要计算A100u0时,就可以使用SΛ100c=c1λ1001x1+c2λ1002x2+⋯+cnλ100nxn。
接下来看一个斐波那契数列(Fibonacci sequence)的例子:
0,1,1,2,3,5,8,13,⋯,F100=?,我们要求第一百项的公式,并观察这个数列是如何增长的。可以想象这个数列并不是稳定数列,因此无论如何该矩阵的特征值并不都小于一,这样才能保持增长。而他的增长速度,则有特征值来决定。
已知Fk+2=Fk1+Fk,但这不是uk+1=Auk的形式,而且我们只要一个方程,而不是方程组,同时这是一个二阶差分方程(就像含有二阶导数的微分方程,希望能够化简为一阶倒数,也就是一阶差分)。
使用一个小技巧,令uk=[Fk+1Fk],再追加一个方程组成方程组:{Fk+2=Fk+1+FkFk+1=Fk+1,再把方程组用矩阵表达得到[Fk+2Fk+1]=[1110][Fk+1Fk],于是我们得到了uk+1=Auk,A=[1110]。我们把二阶标量方程(second-order scalar problem)转化为一阶向量方程组(first-order system)。
我们的矩阵A=[1110]是一个对称矩阵,所以它的特征值将会是实数,且他的特征向量将会互相正交。因为是二阶,我们可以直接利用迹与行列式解方程组{λ1+λ2=1λ1⋅λ2=−1。在求解之前,我们先写出一般解法并观察|A−λI|=|1−λ11−λ|=λ2−λ−1=0,与前面斐波那契数列的递归式Fk+2=Fk+1+Fk→Fk+2−Fk+1−Fk=0比较,我们发现这两个式子在项数与幂次上非常相近。
- 用求根公式解特征值得{λ1=12(1+√5)≈1.618λ2=12(1−√5)≈−0.618,得到两个不同的特征值,一定会有两个线性无关的特征向量,则该矩阵可以被对角化。
我们先来观察这个数列是如何增长的,数列增长由什么来控制?——特征值。哪一个特征值起决定性作用?——较大的一个。
F100=c1(1+√52)100+c2(1−√52)100≈c1(1+√52)100,由于−0.618在幂增长中趋近于0,所以近似的忽略该项,剩下较大的项,我们可以说数量增长的速度大约是1.618。可以看出,这种问题与求解Ax=b不同,这是一个动态的问题,A的幂在不停的增长,而问题的关键就是这些特征值。
- 继续求解特征向量,A−λI=[1−λ111−λ],因为有根式且矩阵只有二阶,我们直接观察[1−λ111−λ][??]=0,由于λ2−λ−1=0,则其特征向量为[λ1],即x1=[λ11],x2=[λ21]。
最后,计算初始项u0=[F1F0]=[10],现在将初始项用特征向量表示出来[10]=c1x1+c2x2,计算系数得c1=√55,c2=−√55。
来回顾整个问题,对于动态增长的一阶方程组,初始向量是u0,关键在于确定A的特征值及特征向量。特征值将决定增长的趋势,发散至无穷还是收敛于某个值。接下来需要找到一个展开式,把u0展开成特征向量的线性组合。
-
再下来就是套用公式,即A的k次方表达式Ak=SΛkS−1,则有u99=Au98=⋯=A99u0=SΛ99S−1Sc=SΛ99c,代入特征值、特征向量得u99=[F100F99]=[1+√521−√5211][(1+√52)9900(1−√52)99][√55−√55]=[c1λ1001+c2λ1002c1λ991+c2λ992],最终结果为F100=c1λ1001+c2λ1002。
-
原式的通解为uk=c1λkx1+c2λkx2。
下一讲将介绍求解微分方程。
第二十三讲:微分方程和eAt
微分方程dudt=Au
本讲主要讲解解一阶方程(first-order system)一阶倒数(first derivative)常系数(constant coefficient)线性方程,上一讲介绍了如何计算矩阵的幂,本讲将进一步涉及矩阵的指数形式。我们通过解一个例子来详细介绍计算方法。
有方程组{du1dt=−u1+2u2du2dt=u1−2u2,则系数矩阵是A=[−121−2],设初始条件为在0时刻u(0)=[u1u2]=[10]。
-
这个初始条件的意义可以看做在开始时一切都在u1中,但随着时间的推移,将有du2dt>0,因为u1项初始为正,u1中的事物会流向u2。随着时间的发展我们可以追踪流动的变化。
-
根据上一讲所学的知识,我们知道第一步需要找到特征值与特征向量。A=[−121−2],很明显这是一个奇异矩阵,所以第一个特征值是λ1=0,另一个特征向量可以从迹得到tr(A)=−3。当然我们也可以用一般方法计算|A−λI|=|−1−λ21−2−λ|=λ2+3λ=0。
(教授提前剧透,特征值λ2=−3将会逐渐消失,因为答案中将会有一项为e−3t,该项会随着时间的推移趋近于0。答案的另一部分将有一项为e0t,该项是一个常数,其值为1,并不随时间而改变。通常含有0特征值的矩阵会随着时间的推移达到稳态。)
-
求特征向量,λ1=0时,即求A的零空间,很明显x1=[21];λ2=−3时,求A+3I的零空间,[2211]的零空间为x2=[1−1]。
-
则方程组的通解为:u(t)=c1eλ1tx1+c2eλ2tx2,通解的前后两部分都是该方程组的纯解,即方程组的通解就是两个与特征值、特征向量相关的纯解的线性组合。我们来验证一下,比如取u=eλ1tx1带入dudt=Au,对时间求导得到λ1eλ1tx1=Aeλ1tx1,化简得λ1x1=Ax1。
对比上一讲,解uk+1=Auk时得到uk=c1λkx1+c2λkx2,而解dudt=Au我们得到u(t)=c1eλ1tx1+c2eλ2tx2。
-
继续求c1,c2,u(t)=c1⋅1⋅[21]+c2⋅e−3t⋅[1−1],已知t=0时,[10]=c1[21]+c2[1−1](Sc=u(0)),所以c1=13,c2=13。
-
于是我们写出最终结果,u(t)=13[21]+13e−3t[1−1]。
稳定性:这个流动过程从u(0)=[10]开始,初始值1的一部分流入初始值0中,经过无限的时间最终达到稳态u(∞)=[2313]。所以,要使得u(t)→0,则需要负的特征值。但如果特征值为复数呢?如λ=−3+6i,我们来计算|e(−3+6i)t|,其中的|e6it|部分为|cos6t+isin6t|=1,因为这部分的模为cos2α+sin2α=1,这个虚部就在单位圆上转悠。所以只有实数部分才是重要的。所以我们可以把前面的结论改为需要实部为负数的特征值。实部会决定最终结果趋近于0或∞,虚部不过是一些小杂音。
收敛态:需要其中一个特征值实部为0,而其他特征值的实部皆小于0。
发散态:如果某个特征值实部大于0。上面的例子中,如果将A变为−A,特征值也会变号,结果发散。
再进一步,我们想知道如何从直接判断任意二阶矩阵的特征值是否均小于零。对于二阶矩阵A=[abcd],矩阵的迹为a+d=λ1+λ2,如果矩阵稳定,则迹应为负数。但是这个条件还不够,有反例迹小于0依然发散:[−2001],迹为−1但是仍然发散。还需要加上一个条件,因为detA=λ1⋅λ2,所以还需要行列式为正数。
总结:原方程组有两个相互耦合的未知函数,u1,u2相互耦合,而特征值和特征向量的作则就是解耦,也就是对角化(diagonalize)。回到原方程组dudt=Au,将u表示为特征向量的线性组合u=Sv,代入原方程有Sdvdt=ASv,两边同乘以S−1得dvdt=S−1ASv=Λv。以特征向量为基,将u表示为Sv,得到关于v的对角化方程组,新方程组不存在耦合,此时{dv1dt=λ1v1dv2dt=λ2v2⋮⋮dvndt=λnvn,这是一个各未知函数间没有联系的方程组,它们的解的一般形式为v(t)=eΛtv(0),则原方程组的解的一般形式为u(t)=eAtu(0)=SeΛtS−1u(0)。这里引入了指数部分为矩阵的形式。
指数矩阵eAt
在上面的结论中,我们见到了eAt。这种指数部分带有矩阵的情况称为指数矩阵(exponential matrix)。
理解指数矩阵的关键在于,将指数形式展开称为幂基数形式,就像ex=1+x22+x36+⋯一样,将eAt展开成幂级数的形式为:
再说些题外话,有两个极具美感的泰勒级数:ex=∑xnn!与11−x=∑xn,如果把第二个泰勒级数写成指数矩阵形式,有(I−At)−1=I+At+(At)2+(At)3+⋯,这个式子在t非常小的时候,后面的高次项近似等于零,所以可以用来近似I−At的逆矩阵,通常近似为I+At,当然也可以再加几项。第一个级数对我们而言比第二个级数好,因为第一个级数总会收敛于某个值,所以ex总会有意义,而第二个级数需要A特征值的绝对值小于1(因为涉及矩阵的幂运算)。我们看到这些泰勒级数的公式对矩阵同样适用。
回到正题,我们需要证明SeΛtS−1=eAt,继续使用泰勒级数:
需要注意的是,eAt的泰勒级数展开是恒成立的,但我们推出的版本却需要矩阵可对角化这个前提条件。
最后,我们来看看什么是eΛt,我们将eAt变为对角矩阵就是因为对角矩阵简单、没有耦合,eΛt=[eλ1t0⋯00eλ2t⋯0⋮⋮⋱⋮00⋯eλnt]。
有了u(t)=SeΛtS−1u(0),再来看矩阵的稳定性可知,所有特征值的实部均为负数时矩阵收敛,此时对角线上的指数收敛为0。如果我们画出复平面,则要使微分方程存在稳定解,则特征值存在于复平面的左侧(即实部为负);要使矩阵的幂收敛于0,则特征值存在于单位圆内部(即模小于1),这是幂稳定区域。(上一讲的差分方程需要计算矩阵的幂。)
同差分方程一样,我们来看二阶情况如何计算,有y″+by′+k=0。我们也模仿差分方程的情形,构造方程组{y″=−by′−kyy′=y′,写成矩阵形式有[y″y′]=[−b−k10][y′y],令u′=[y″y′], u=[y′y]。
继续推广,对于5阶微分方程y′′′′′+by⁗+cy‴+dy″+ey′+f=0,则可以写作[y′′′′′y⁗y‴y″y′]=[−b−c−d−e−f10000010000010000010][y⁗y‴y″y′y],这样我们就把一个五阶微分方程化为5×5一阶方程组了,然后就是求特征值、特征向量了步骤了。
第二十四讲:马尔科夫矩阵、傅里叶级数
马尔科夫矩阵
马尔科夫矩阵(Markov matrix)是指具有以下两个特性的矩阵:
- 矩阵中的所有元素大于等于0;(因为马尔科夫矩阵与概率有关,而概率是非负的。)
- 每一列的元素之和为1
对于马尔科夫矩阵,我们关心幂运算过程中的稳态(steady state)。与上一讲不同,指数矩阵关系特征值是否为0,而幂运算要达到稳态需要特征值为1。
根据上面两条性质,我们可以得出两个推论:
- 马尔科夫矩阵必有特征值为1;
- 其他的特征值的绝对值皆小于1。
使用第二十二讲中得到的公式进行幂运算uk=Aku0=SΛkS−1u0=SΛkS−1Sc=SΛkc=c1λk1x1+c2λk2x2+⋯+cnλknxn,从这个公式很容易看出幂运算的稳态。比如我们取λ1=1,其他的特征值绝对值均小于1,于是在经过k次迭代,随着时间的推移,其他项都趋近于0,于是在k→∞时,有稳态uk=c1x1,这也就是初始条件u0的第1个分量。
我们来证明第一个推论,取A=[0.10.010.30.20.990.30.700.4],则A−I=[−0.90.010.30.2−0.010.30.70−0.6]。观察A−I易知其列向量中元素之和均为0,因为马尔科夫矩阵的性质就是各列向量元素之和为1,现在我们从每一列中减去了1,所以这是很自然的结果。而如果列向量中元素和为0,则矩阵的任意行都可以用“零减去其他行之和”表示出来,即该矩阵的行向量线性相关。
用以前学过的子空间的知识描述,当n阶方阵各列向量元素之和皆为1时,则有[11⋮1]在矩阵A−I左零空间中,即(A−I)T行向量线性相关。而A特征值1所对应的特征向量将在A−I的零空间中,因为Ax=x→(A−I)x=0。
另外,特征值具有这样一个性质:矩阵与其转置的特征值相同。因为我们在行列式一讲了解了性质10,矩阵与其转置的行列式相同,那么如果det(A−λI)=0,则有det(A−λI)T=0,根据矩阵转置的性质有det(AT−λIT)=0,即det(AT−λI)=0。这正是AT特征值的计算式。
然后计算特征值λ1=1所对应的特征向量,(A−I)x1=0,得出x1=[0.6330.7],特征向量中的元素皆为正。
接下来介绍马尔科夫矩阵的应用,我们用麻省和加州这两个州的人口迁移为例:
[ucalumass]k+1[0.90.20.10.8][ucalumass]k,元素非负,列和为一。这个式子表示每年有10的人口从加州迁往麻省,同时有20的人口从麻省迁往加州。注意使用马尔科夫矩阵的前提条件是随着时间的推移,矩阵始终不变。
设初始情况[ucalumass]0=[01000],我们先来看第一次迁徙后人口的变化情况:[ucalumass]1=[0.90.20.10.8][01000]=[200800],随着时间的推移,会有越来越多的麻省人迁往加州,而同时又会有部分加州人迁往麻省。
计算特征值:我们知道马尔科夫矩阵的一个特征值为λ1=1,则另一个特征值可以直接从迹算出λ2=0.7。
计算特征向量:带入λ1=1求A−I的零空间有[−0.10.20.1−0.2],则x1=[21],此时我们已经可以得出无穷步后稳态下的结果了。u∞=c1[21]且人口总数始终为1000,则c1=10003,稳态时[ucalumass]∞=[2000310003]。注意到特征值为1的特征向量元素皆为正。
为了求每一步的结果,我们必须解出所有特征向量。带入λ2=0.7求A−0.7I的零空间有[0.20.20.10.1],则x2=[−11]。
通过u0解出c1,c2,uk=c11k[21]+c20.7k[−11],带入k=0得u0=[01000]=c1[21]+c2[−11],解出c1=10003,c2=20003。
另外,有时人们更喜欢用行向量,此时将要使用行向量乘以矩阵,其行向量各分量之和为1。
傅里叶级数
在介绍傅里叶级数(Fourier series)之前,先来回顾一下投影。
设q1,q2,⋯qn为一组标准正交基,则向量v在该标准正交基上的展开为v=x1q1+x2q2+⋯+xnqn,此时我们想要得到各系数xi的值。比如求x1的值,我们自然想要消掉除x1q1外的其他项,这时只需要等式两边同乘以qT1,因为的qi向量相互正交且长度为1,则qTiqj=0,q2i=1所以原式变为qT1v=x1。
写为矩阵形式有[q1 q2 ⋯ qn][x1x2⋮xn]=v,即Qx=v。所以有x=Q−1v,而在第十七讲我们了解到标准正交基有QT=Q−1,所以我们不需要计算逆矩阵可直接得出x=QTv。此时对于x的每一个分量有xi=qTiv。
接下来介绍傅里叶级数。先写出傅里叶级数的展开式:
傅里叶发现,如同将向量v展开(投影)到向量空间的一组标准正交基中,在函数空间中,我们也可以做类似的展开。将函数f(x)投影在一系列相互正交的函数中。函数空间中的f(x)就是向量空间中的v;函数空间中的1,cosx,sinx,cos2x,sin2x,⋯就是向量空间中的q1,q2,⋯,qn;不同的是,函数空间是无限维的而我们以前接触到的向量空间通常是有限维的。
再来介绍何为“函数正交”。对于向量正交我们通常使用两向量内积(点乘)为零判断。我们知道对于向量v,w的内积为vTw=v1w1+v2w2+⋯+vnwn=0,也就是向量的每个分量之积再求和。而对于函数f(x)⋅g(x)内积,同样的,我们需要计算两个函数的每个值之积而后求和,由于函数取值是连续的,所以函数内积为:
在本例中,由于傅里叶级数使用正余弦函数,它们的周期都可以算作2π,所以本例的函数点积可以写作fTg=∫2π0f(x)g(x)dx。我来检验一个内积∫2π0sinxcosxdx=12sin2x|2π0=0,其余的三角函数族正交性结果可以参考傅里叶级数的“希尔伯特空间的解读”一节。
最后我们来看cosx项的系数是多少(a0是f(x)的平均值)。同向量空间中的情形一样,我们在等式两边同时做cosx的内积,原式变为∫2π0f(x)cosxdx=a1∫2π0cos2xdx,因为正交性等式右边仅有cosx项不为零。进一步化简得a1π=∫2π0f(x)cosxdx→a1=1π∫2π0f(x)cosxdx。
于是,我们把函数f(x)展开到了函数空间的一组标准正交基上。
第二十五讲:复习二
- 我们学习了正交性,有矩阵Q=[q1 q2 ⋯ qn],若其列向量相互正交,则该矩阵满足QTQ=I。
- 进一步研究投影,我们了解了Gram-Schmidt正交化法,核心思想是求法向量,即从原向量中减去投影向量E=b−P,P=Ax=ATbATA⋅A。
- 接着学习了行列式,根据行列式的前三条性质,我们拓展出了性质4-10。
- 我们继续推导出了一个利用代数余子式求行列式的公式。
- 又利用代数余子式推导出了一个求逆矩阵的公式。
- 接下来我们学习了特征值与特征向量的意义:Ax=λx,进而了解了通过det(A−λI)=0求特征值、特征向量的方法。
- 有了特征值与特征向量,我们掌握了通过公式AS=ΛS对角化矩阵,同时掌握了求矩阵的幂Ak=SΛkS−1。
微分方程不在本讲的范围内。下面通过往年例题复习上面的知识。
-
求a=[212]的投影矩阵P:(由a⊥(b−p)→AT(b−Aˆx)=0得到ˆx=(ATA)−1ATb,求得p=Aˆx=A(ATA)−1ATb=Pb最终得到P)P=A(ATA)−1AT_a=aaTaTa=19[424212424]。
求P矩阵的特征值:观察矩阵易知矩阵奇异,且为秩一矩阵,则其零空间为2维,所以由Px=0x得出矩阵的两个特征向量为λ1=λ2=0;而从矩阵的迹得知trace(P)=1=λ1+λ2+λ3=0+0+1,则第三个特征向量为λ3=1。
求λ3=1的特征向量:由Px=x我们知道经其意义为,x过矩阵P变换后不变,又有P是向量a的投影矩阵,所以任何向量经过P变换都会落在a的列空间中,则只有已经在a的列空间中的向量经过P的变换后保持不变,即其特征向量为x=a=[212],也就是Pa=a。
有差分方程uk+1=Puk, u0=[990],求解uk:我们先不急于解出特征值、特征向量,因为矩阵很特殊(投影矩阵)。首先观察u1=Pu0,式子相当于将u0投影在了a的列空间中,计算得u1=aaTu0aTa=3a=[636](这里的3相当于做投影时的系数ˆx),其意义为u1在a上且距离u0最近。再来看看u2=Pu1,这个式子将u1再次投影到a的列空间中,但是此时的u1已经在该列空间中了,再次投影仍不变,所以有uk=Pku0=Pu0=[636]。
上面的解法利用了投影矩阵的特殊性质,如果在一般情况下,我们需要使用AS=SΛ→A=SΛS−1→uk+1=Auk=Ak+1u0,u0=Sc→uk+1=SΛk+1S−1Sc=SΛk+1c,最终得到公式Aku0=c1λk1x1+c2λk2x2+⋯+cnλknxn。题中P的特殊性在于它的两个“零特征值”及一个“一特征值”使得式子变为Aku0=c3x3,所以得到了上面结构特殊的解。
-
将点(1,4), (2,5), (3,8)拟合到一条过零点的直线上:设直线为y=Dt,写成矩阵形式为[123]D=[458],即AD=b,很明显D不存在。利用公式ATAˆD=ATb得到14D=38, ˆD=3814,即最佳直线为y=3814t。这个近似的意义是将b投影在了A的列空间中。
-
求a1=[123] a2=[111]的正交向量:找到平面A=[a1,a2]的正交基,使用Gram-Schmidt法,以a1为基准,正交化a2,也就是将a2中平行于a1的分量去除,即a2−xa1=a2−aT1a2aT1a1a1=[111]−614[123]。
-
有4×4矩阵A,其特征值为λ1,λ2,λ3,λ4,则矩阵可逆的条件是什么:矩阵可逆,则零空间中只有零向量,即Ax=0x没有非零解,则零不是矩阵的特征值。
detA−1是什么:detA−1=1detA,而detA=λ1λ2λ3λ4,所以有detA−1=1λ1λ2λ3λ4。
trace(A+I)的迹是什么:我们知道trace(A)=a11+a22+a33+a44=λ1+λ2+λ3+λ4,所以有trace(A+I)=a11+1+a22+1+a33+1+a44+1=λ1+λ2+λ3+λ4+4。
-
有矩阵A4=[1100111001110011],求Dn=?Dn−1+?Dn−2:求递归式的系数,使用代数余子式将矩阵安第一行展开得detA4=1⋅|110111011|−1⋅|110011011|=1⋅|110111011|−1⋅|1111|=detA3−detA2。则可以看出有规律Dn=Dn−1−Dn−2,D1=1,D2=0。
使用我们在差分方程中的知识构建方程组{Dn=Dn−1−Dn−2Dn−1=Dn−1,用矩阵表达有[DnDn−1]=[1−110][Dn−1Dn−2]。计算系数矩阵Ac的特征值,|1−λ11−λ|=λ2−λ+1=0,解得λ1=1+√3i2,λ2=1−√3i2,特征值为一对共轭复数。
要判断递归式是否收敛,需要计算特征值的模,即实部平方与虚部平方之和14+34=1。它们是位于单位圆eiθ上的点,即cosθ+isinθ,从本例中可以计算出θ=60∘,也就是可以将特征值写作λ1=eiπ/3,λ2=e−iπ/3。注意,从复平面单位圆上可以看出,这些特征值的六次方将等于一:e2πi=e2πi=1。继续深入观察这一特性对矩阵的影响,λ61=λ6=1,则对系数矩阵有A6c=I。则系数矩阵Ac服从周期变化,既不发散也不收敛。
-
有这样一类矩阵A4=[0100102002030030],求投影到A3列空间的投影矩阵:有A3=[010102020],按照通常的方法求P=A(ATA)AT即可,但是这样很麻烦。我们可以考察这个矩阵是否可逆,因为如果可逆的话,R4空间中的任何向量都会位于A4的列空间,其投影不变,则投影矩阵为单位矩阵I。所以按行展开求行列式detA4=−1⋅−1⋅−3⋅−3=9,所以矩阵可逆,则P=I。
求A3的特征值及特征向量:|A3−λI|=|−λ101−λ202−λ|=−λ3+5λ=0,解得λ1=0,λ2=√5,λ3=−√5。
我们可以猜测这一类矩阵的规律:奇数阶奇异,偶数阶可逆。
第二十六讲:对称矩阵及正定性
对称矩阵
前面我们学习了矩阵的特征值与特征向量,也了解了一些特殊的矩阵及其特征值、特征向量,特殊矩阵的特殊性应该会反映在其特征值、特征向量中。如马尔科夫矩阵,有一特征值为1,本讲介绍(实)对称矩阵。
先提前介绍两个对称矩阵的特性:
- 特征值为实数;(对比第二十一讲介绍的旋转矩阵,其特征值为纯虚数。)
- 特征向量相互正交。(当特征值重复时,特征向量也可以从子空间中选出相互正交正交的向量。)
典型的状况是,特征值不重复,特征向量相互正交。
- 那么在通常(可对角化)情况下,一个矩阵可以化为:A=SΛS−1;
- 在矩阵对称的情况下,通过性质2可知,由特征向量组成的矩阵S中的列向量是相互正交的,此时如果我们把特征向量的长度统一化为1,就可以得到一组标准正交的特征向量。则对于对称矩阵有A=QΛQ−1,而对于标准正交矩阵,有Q=QT,所以对称矩阵可以写为A=QΛQT
观察(1)式,我们发现这个分解本身就代表着对称,(QΛQT)T=(QT)TΛTQT=QΛQT。(1)式在数学上叫做谱定理(spectral theorem),谱就是指矩阵特征值的集合。(该名称来自光谱,指一些纯事物的集合,就像将特征值分解成为特征值与特征向量。)在力学上称之为主轴定理(principle axis theorem),从几何上看,它意味着如果给定某种材料,在合适的轴上来看,它就变成对角化的,方向就不会重复。
-
现在我们来证明性质1。对于矩阵Ax=λx_,对于其共轭部分总有ˉAˉx=ˉλˉx,根据前提条件我们只讨论实矩阵,则有Aˉx=ˉλˉx,将等式两边取转置有¯ˉxTA=ˉxTˉλ。将“下划线”式两边左乘ˉxT有ˉxTAx=ˉxTλx,“上划线”式两边右乘x有ˉxTAx=ˉxTˉλx,观察发现这两个式子左边是一样的,所以ˉxTλx=ˉxTˉλx,则有λ=ˉλ(这里有个条件,ˉxTx≠0),证毕。
观察这个前提条件,ˉxTx=[ˉx1ˉx2⋯ˉxn][x1x2⋮xn]=ˉx1x1+ˉx2x2+⋯+ˉxnxn,设x1=a+ib,ˉx1=a−ib则ˉx1x1=a2+b2,所以有ˉxTx>0。而ˉxTx就是x长度的平方。
拓展这个性质,当A为复矩阵,根据上面的推导,则矩阵必须满足A=ˉAT时,才有性质1、性质2成立(教授称具有这种特征值为实数、特征向量相互正交的矩阵为“好矩阵”)。
继续研究A=QΛQT=[q1 q2 ⋯ qn][λ1⋯λ2⋯⋮⋮⋱⋮⋯λn][qT1qT1⋮qT1]=λ1q1qT1+λ2q2qT2+⋯+λnqnqTn,注意这个展开式中的qqT,q是单位列向量所以qTq=1,结合我们在第十五讲所学的投影矩阵的知识有qqTqTq=qqT是一个投影矩阵,很容易验证其性质,比如平方它会得到qqTqqT=qqT于是多次投影不变等。
每一个对称矩阵都可以分解为一系列相互正交的投影矩阵。
在知道对称矩阵的特征值皆为实数后,我们再来讨论这些实数的符号,因为特征值的正负号会影响微分方程的收敛情况(第二十三讲,需要实部为负的特征值保证收敛)。用消元法取得矩阵的主元,观察主元的符号,主元符号的正负数量与特征向量的正负数量相同。
正定性
如果对称矩阵是“好矩阵”,则正定矩阵(positive definite)是其一个更好的子类。正定矩阵指特征值均为正数的矩阵(根据上面的性质有矩阵的主元均为正)。
举个例子,[5223],由行列式消元知其主元为5,115,按一般的方法求特征值有|5−λ223−lambda|=λ2−8λ+11=0,λ=4±√5。
正定矩阵的另一个性质是,所有子行列式为正。对上面的例子有|5|=5,|5223|=11。
我们看到正定矩阵将早期学习的的消元主元、中期学习的的行列式、后期学习的特征值结合在了一起。
第二十七讲:复数矩阵和快速傅里叶变换
本讲主要介绍复数向量、复数矩阵的相关知识(包括如何做复数向量的点积运算、什么是复数对称矩阵等),以及傅里叶矩阵(最重要的复数矩阵)和快速傅里叶变换。
复数矩阵运算
先介绍复数向量,我们不妨换一个字母符号来表示:z=[z1z2⋮zn],向量的每一个分量都是复数。此时z不再属于Rn实向量空间,它现在处于Cn复向量空间。
计算复向量的模
对比实向量,我们计算模只需要计算|v|=√vTv即可,而如果对复向量使用zTz则有zTz=[z1z2⋯zn][z1z2⋮zn]=z21+z22+⋯+z2n,这里zi是复数,平方后虚部为负,求模时本应相加的运算变成了减法。(如向量[1i],右乘其转置后结果为0,但此向量的长度显然不是零。)
根据上一讲我们知道,应使用|z|=√ˉzTz,即[ˉz1ˉz2⋯ˉzn][z1z2⋮zn],即使用向量共轭的转置乘以原向量即可。(如向量[1i],右乘其共轭转置后结果为[1−i][1i]=2。)
我们把共轭转置乘以原向量记为zHz,H读作埃尔米特(人名为Hermite,形容词为Hermitian)
计算向量的内积
有了复向量模的计算公式,同理可得,对于复向量,内积不再是实向量的yTx形式,复向量内积应为yHx。
对称性
对于实矩阵,AT=A即可表达矩阵的对称性。而对于复矩阵,我们同样需要求一次共轭ˉAT=A。举个例子[23+i3−i5]是一个复数情况下的对称矩阵。这叫做埃尔米特矩阵,有性质AH=A。
正交性
在第十七讲中,我们这样定义标准正交向量:qTiqj={0i≠j1i=j。现在,对于复向量我们需要求共轭:ˉqTiqj=qHiqj={0i≠j1i=j。
第十七讲中的标准正交矩阵:Q=[q1 q2 ⋯ qn]有QTQ=I。现在对于复矩阵则有QHQ=I。
就像人们给共轭转置起了个“埃尔米特”这个名字一样,正交性(orthogonal)在复数情况下也有了新名字,酉(unitary),酉矩阵(unitary matrix)与正交矩阵类似,满足QHQ=I的性质。而前面提到的傅里叶矩阵就是一个酉矩阵。
傅里叶矩阵
n阶傅里叶矩阵Fn=[111⋯11ww2⋯wn−11w2w4⋯w2(n−1)⋮⋮⋮⋱⋮1wn−1w2(n−1)⋯w(n−1)2],对于每一个元素有(Fn)ij=wiji,j=0,1,2,⋯,n−1。矩阵中的w是一个非常特殊的值,满足wn=1,其公式为w=ei2π/n。易知w在复平面的单位圆上,w=cos2πn+isin2πn。
在傅里叶矩阵中,当我们计算w的幂时,w在单位圆上的角度翻倍。比如在6阶情形下,w=e2π/6,即位于单位圆上60∘角处,其平方位于单位圆上120∘角处,而w6位于1处。从开方的角度看,它们是1的6个六次方根,而一次的w称为原根。
-
我们现在来看4阶傅里叶矩阵,先计算w有w=i, w2=−1, w3=−i, w4=1,F4=[11111ii2i31i2i4i61i3i6i9]=[11111i−1−i1−11−11−i−1i]。
矩阵的四个列向量正交,我们验证一下第二列和第四列,¯c2Tc4=1−0+1−0=0,正交。不过我们应该注意到,F4的列向量并不是标准的,我们可以给矩阵乘上系数12(除以列向量的长度)得到标准正交矩阵F4=12[11111i−1−i1−11−11−i−1i]。此时有FH4F4=I,于是该矩阵的逆矩阵也就是其共轭转置FH4。
快速傅里叶变换(Fast Fourier transform/FFT)
对于傅里叶矩阵,F6, F3、F8, F4、F64, F32之间有着特殊的关系。
举例,有傅里叶矩阵F64,一般情况下,用一个列向量右乘F64需要约642次计算,显然这个计算量是比较大的。我们想要减少计算量,于是想要分解F64,联系到F32,有[F64]=[IDI−D][F3200F32][1⋯0⋯0⋯1⋯1⋯0⋯0⋯1⋯⋱⋱⋱⋱⋯1⋯0⋯0⋯1]。
我们分开来看等式右侧的这三个矩阵:
-
第一个矩阵由单位矩阵I和对角矩阵D=[1ww2⋱w31]组成,我们称这个矩阵为修正矩阵,显然其计算量来自D矩阵,对角矩阵的计算量约为32即这个修正矩阵的计算量约为32,单位矩阵的计算量忽略不计。
-
第二个矩阵是两个F32与零矩阵组成的,计算量约为2×322。
-
第三个矩阵通常记为P矩阵,这是一个置换矩阵,其作用是讲前一个矩阵中的奇数列提到偶数列之前,将前一个矩阵从[x0 x1 ⋯]变为[x0 x2 ⋯ x1 x3 ⋯],这个置换矩阵的计算量也可以忽略不计。(这里教授似乎在黑板上写错了矩阵,可以参考FFT、How the FFT is computed做进一步讨论。)
所以我们把642复杂度的计算化简为2×322+32复杂度的计算,我们可以进一步化简F32得到与F16有关的式子[I32D32I32−D32][I16D16I16−D16I16D16I16−D16][F16F16F16F16][P16P16][ P32 ]。而322的计算量进一步分解为2×162+16的计算量,如此递归下去我们最终得到含有一阶傅里叶矩阵的式子。
来看化简后计算量,2(2(2(2(2(2(1)2+1)+2)+4)+8)+16)+32,约为6×32=log264×642,算法复杂度为n2log2n。
于是原来需要n2的运算现在只需要n2log2n就可以实现了。不妨看看n=10的情况,不使用FFT时需要n2=1024×1024次运算,使用FFT时只需要n2log2n=5×1024次运算,运算量大约是原来的1200。
下一讲将继续介绍特征值、特征向量及正定矩阵。
第二十八讲:正定矩阵和最小值
本讲我们会了解如何完整的测试一个矩阵是否正定,测试xTAx是否具有最小值,最后了解正定的几何意义——椭圆(ellipse)和正定性有关,双曲线(hyperbola)与正定无关。另外,本讲涉及的矩阵均为实对称矩阵。
正定性的判断
我们仍然从二阶说起,有矩阵A=[abbd],判断其正定性有以下方法:
- 矩阵的所有特征值大于零则矩阵正定:λ1>0, λ2>0;
- 矩阵的所有顺序主子阵(leading principal submatrix)的行列式(即顺序主子式,leading principal minor)大于零则矩阵正定:a>0, ac−b2>0;
- 矩阵消元后主元均大于零:a>0, ac−b2a>0;
- xTAx>0;
大多数情况下使用4来定义正定性,而用前三条来验证正定性。
来计算一个例子:A=[266?],在?处填入多少才能使矩阵正定?
-
来试试18,此时矩阵为A=[26618],detA=0,此时的矩阵成为半正定矩阵(positive semi-definite)。矩阵奇异,其中一个特征值必为0,从迹得知另一个特征值为20。矩阵的主元只有一个,为2。
计算xTAx,得[x1x2][26618][x1x2]=2x21+12x1x2+18x22这样我们得到了一个关于x1,x2的函数f(x1,x2)=2x21+12x1x2+18x22,这个函数不再是线性的,在本例中这是一个纯二次型(quadratic)函数,它没有线性部分、一次部分或更高次部分(Ax是线性的,但引入xT后就成为了二次型)。
当?取18时,判定1、2、3都是“刚好不及格”。
-
我们可以先看“一定不及格”的样子,令?=7,矩阵为A=[2667],二阶顺序主子式变为−22,显然矩阵不是正定的,此时的函数为f(x1,x2)=2x21+12x1x2+7x22,如果取x1=1,x2=−1则有f(1,−1)=2−12+7<0。
如果我们把z=2x2+12xy+7y2放在直角坐标系中,图像过原点z(0,0)=0,当y=0或x=0或x=y时函数为开口向上的抛物线,所以函数图像在某些方向上是正值;而在某些方向上是负值,比如x=−y,所以函数图像是一个马鞍面(saddle),(0,0,0)点称为鞍点(saddle point),它在某些方向上是极大值点,而在另一些方向上是极小值点。(实际上函数图像的最佳观测方向是沿着特征向量的方向。)
-
再来看一下“一定及格”的情形,令?=20,矩阵为A=[26620],行列式为detA=4,迹为trace(A)=22,特征向量均大于零,矩阵可以通过测试。此时的函数为f(x1,x2)=2x21+12x1x2+20x22,函数在除(0,0)外处处为正。我们来看看z=2x2+12xy+20y2的图像,式子的平方项均非负,所以需要两个平方项之和大于中间项即可,该函数的图像为抛物面(paraboloid)。在(0,0)点函数的一阶偏导数均为零,二阶偏导数均为正(马鞍面的一阶偏导数也为零,但二阶偏导数并不均为正,所以),函数在改点取极小值。
在微积分中,一元函数取极小值需要一阶导数为零且二阶导数为正dudx=0,d2udx2>0。在线性代数中我们遇到了了多元函数f(x1,x2,⋯,xn),要取极小值需要二阶偏导数矩阵为正定矩阵。
在本例中(即二阶情形),如果能用平方和的形式来表示函数,则很容易看出函数是否恒为正,f(x,y)=2x2+12xy+20y2=2(x+3y)2+2y2。另外,如果是上面的?=7的情形,则有f(x,y)=2(x+3y)2−11y2,如果是?=18的情形,则有f(x,y)=2(x+3y)2。
如果令z=1,相当于使用z=1平面截取该函数图像,将得到一个椭圆曲线。另外,如果在?=7的马鞍面上截取曲线将得到一对双曲线。
再来看这个矩阵的消元,[26620]=[10−31][2602],这就是A=LU,可以发现矩阵L中的项与配平方中未知数的系数有关,而主元则与两个平方项外的系数有关,这也就是为什么正数主元得到正定矩阵。
上面又提到二阶导数矩阵,这个矩阵型为[fxxfxyfyxfyy],显然,矩阵中的主对角线元素(纯二阶导数)必须为正,并且主对角线元素必须足够大来抵消混合导数的影响。同时还可以看出,因为二阶导数的求导次序并不影响结果,所以矩阵必须是对称的。现在我们就可以计算n×n阶矩阵了。
接下来计算一个三阶矩阵,A=[2−10−12−10−12],它是正定的吗?函数xTAx是多少?函数在原点去最小值吗?图像是什么样的?
- 先来计算矩阵的顺序主子式,分别为2,3,4;再来计算主元,分别为2,32,43;计算特征值,λ1=2−√2,λ2=2,λ3=2+√2。
- 计算xTAx=2x21+2x22+2x23−2x1x2−2x2x3。
- 图像是四维的抛物面,当我们在f(x1,x2,x3)=1处截取该面,将得到一个椭圆体。一般椭圆体有三条轴,特征值的大小决定了三条轴的长度,而特征向量的方向与三条轴的方向相同。
现在我们将矩阵A分解为A=QΛQT,可以发现上面说到的各种元素都可以表示在这个分解的矩阵中,我们称之为主轴定理(principal axis theorem),即特征向量说明主轴的方向、特征值说明主轴的长度。
A=QΛQT是特征值相关章节中最重要的公式。
第二十九讲:相似矩阵和若尔当形
在本讲的开始,先接着上一讲来继续说一说正定矩阵。
-
正定矩阵的逆矩阵有什么性质?我们将正定矩阵分解为A=SΛS−1,引入其逆矩阵A−1=SΛ−1S−1,我们知道正定矩阵的特征值均为正值,所以其逆矩阵的特征值也必为正值(即原矩阵特征值的倒数)所以,正定矩阵的逆矩阵也是正定的。
-
如果A, B均为正定矩阵,那么A+B呢?我们可以从判定xT(A+B)x入手,根据条件有xTAx>0, xTBx>0,将两式相加即得到xT(A+B)x>0。所以正定矩阵之和也是正定矩阵。
-
再来看有m×n矩阵A,则ATA具有什么性质?我们在投影部分经常使用ATA,这个运算会得到一个对称矩阵,这个形式的运算用数字打比方就像是一个平方,用向量打比方就像是向量的长度平方,而对于矩阵,有ATA正定:在式子两边分别乘向量及其转置得到xTATAx,分组得到(Ax)T(Ax),相当于得到了向量Ax的长度平方,则|Ax|2≥0。要保证模不为零,则需要Ax的零空间中仅有零向量,即A的各列线性无关(rank(A)=n)即可保证|Ax|2>0,ATA正定。
-
另外,在矩阵数值计算中,正定矩阵消元不需要进行“行交换”操作,也不必担心主元过小或为零,正定矩阵具有良好的计算性质。
接下来进入本讲的正题。
相似矩阵
先列出定义:矩阵A, B对于某矩阵M满足B=M−1AM时,成A, B互为相似矩阵。
对于在对角化一讲(第二十二讲)中学过的式子S−1AS=Λ,则有A相似于Λ。
-
举个例子,A=[2112],容易通过其特征值得到相应的对角矩阵Λ=[3001],取M=[1401],则B=M−1AM=[1−401][2112][1401]=[−2−1516]。
我们来计算这几个矩阵的的特征值(利用迹与行列式的性质),λΛ=3, 1、λA=3, 1、λB=3, 1。
所以,相似矩阵有相同的特征值。
- 继续上面的例子,特征值为3, 1的这一族矩阵都是相似矩阵,如[3701]、[1703],其中最特殊的就是Λ。
现在我们来证明这个性质,有Ax=λx, B=M−1AM,第一个式子化为AMM−1x=λx,接着两边同时左乘M−1得M−1AMM−1x=λM−1x,进行适当的分组得(M−1AM)M−1x=λM−1x即BM−1x=λM−1x。
BM−1=λM−1x可以解读成矩阵B与向量M−1x之积等于λ与向量M−1x之积,也就是B的仍为λ,而特征向量变为M−1x。
以上就是我们得到的一族特征值为3, 1的矩阵,它们具有相同的特征值。接下来看特征值重复时的情形。
- 特征值重复可能会导致特征向量短缺,来看一个例子,设λ1=λ2=4,写出具有这种特征值的矩阵中的两个[4004],[4104]。其实,具有这种特征值的矩阵可以分为两族,第一族仅有一个矩阵[4004],它只与自己相似(因为M−1[4004]M=4M−1IM=4I=[4004],所以无论M如何取值该对角矩阵都只与自己相似);另一族就是剩下的诸如[4104]的矩阵,它们都是相似的。在这个“大家族”中,[4104]是“最好”的一个矩阵,称为若尔当形。
若尔当形在过去是线性代数的核心知识,但现在不是了(现在是下一讲的奇异值分解),因为它并不容易计算。
- 继续上面的例子,我们在在出几个这一族的矩阵[4104], [51−13], [40174],我们总是可以构造出一个满足trace(A)=8, detA=16的矩阵,这个矩阵总是在这一个“家族”中。
若尔当形
再来看一个更加“糟糕”的矩阵:
-
矩阵[0100001000000000],其特征值为四个零。很明显矩阵的秩为2,所以其零空间的维数为4−2=2,即该矩阵有两个特征向量。可以发现该矩阵在主对角线的上方有两个1,在对角线上每增加一个1,特征向量个个数就减少一个。
-
令一个例子,[0100000000010000],从特征向量的数目看来这两个矩阵是相似的,其实不然。
若尔当认为第一个矩阵是由一个3×3的块与一个1×1的块组成的 [0100000000010000],而第二个矩阵是由两个2×2矩阵组成的[0100000000010000],这些分块被称为若尔当块。
若尔当块的定义型为Ji=[λi1⋯λi1⋯λi⋯⋮⋮⋮⋱λi],它的对角线上只为同一个数,仅有一个特征向量。
所有有,每一个矩阵A都相似于一个若尔当矩阵,型为J=[J1J2⋱Jd]。注意,对角线上方还有1。若尔当块的个数即为矩阵特征值的个数。
在矩阵为“好矩阵”的情况下,n阶矩阵将有n个不同的特征值,那么它可以对角化,所以它的若尔当矩阵就是Λ,共n个特征向量,有n个若尔当块。
第三十讲:奇异值分解
本讲我们介绍将一个矩阵写为A=UΣVT,分解的因子分别为正交矩阵、对角矩阵、正交矩阵,与前面几讲的分解不同的是,这两个正交矩阵通常是不同的,而且这个式子可以对任意矩阵使用,不仅限于方阵、可对角化的方阵等。
-
在正定一讲中(第二十八讲)我们知道一个正定矩阵可以分解为A=QΛQT的形式,由于A对称性其特征向量是正交的,且其Λ矩阵中的元素皆为正,这就是正定矩阵的奇异值分解。在这种特殊的分解中,我们只需要一个正交矩阵Q就可以使等式成立。
-
在对角化一讲中(第二十二讲),我们知道可对角化的矩阵能够分解为A=SΛST的形式,其中S的列向量由A的特征向量组成,但S并不是正交矩阵,所以这不是我们希望得到的奇异值分解。
我们现在要做的是,在A的列空间中找到一组特殊的正交基v1,v2,⋯,vr,这组基在A的作用下可以转换为A的行空间中的一组正交基u1,u2,⋯,ur。
用矩阵语言描述为A[v1 v2 ⋯ vr]=[σ1u1 σ2u2 ⋯ σrur]=[u1 u2 ⋯ ur][σ1σ2⋱σn],即Av1=σ1u1, Av2=σ2u2,⋯,Avr=σrur,这些σ是缩放因子,表示在转换过程中有拉伸或压缩。而A的左零空间和零空间将体现在σ的零值中。
另外,如果算上左零、零空间,我们同样可以对左零、零空间取标准正交基,然后写为A[v1 v2 ⋯ vr vr+1 ⋯ vm]=[u1 u2 ⋯ ur ur+1 ⋯ un][σ1⋱σr[0]],此时U是m×m正交矩阵,Σ是m×n对角矩阵,VT是n×n正交矩阵。
最终可以写为AV=UΣ,可以看出这十分类似对角化的公式,矩阵A被转化为对角矩阵Σ,我们也注意到U, V是两组不同的正交基。(在正定的情况下,U, V都变成了Q。)。进一步可以写作A=UΣV−1,因为V是标准正交矩阵所以可以写为A=UΣVT
计算一个例子,A=[44−33],我们需要找到:
- 行空间R2的标准正交基v1,v2;
- 列空间R2的标准正交基u1,u2;
- σ1>0,σ2>0。
在A=UΣVT中有两个标准正交矩阵需要求解,我们希望一次只解一个,如何先将U消去来求V?
这个技巧会经常出现在长方形矩阵中:求ATA,这是一个对称正定矩阵(至少是半正定矩阵),于是有ATA=VΣTUTUΣVT,由于U是标准正交矩阵,所以UTU=I,而ΣTΣ是对角线元素为σ2的对角矩阵。
现在有ATA=V[σ1σ2⋱σn]VT,这个式子中V即是ATA的特征向量矩阵而Σ2是其特征值矩阵。
同理,我们只想求U时,用AAT消掉V即可。
我们来计算ATA=[4−343][44−33]=[257725],对于简单的矩阵可以直接观察得到特征向量ATA[11]=32[11], ATA[1−1]=18[1−1],化为单位向量有σ1=32, v1=[1√21√2], σ2=18, v2=[1√2−1√2]。
到目前为止,我们得到[44−33]=[u?u?u?u?][√3200√18][1√21√21√2−1√2],接下来继续求解U。
AAT=UΣVTVΣTUT=UΣ2UT,求出AAT的特征向量即可得到U,[44−33][4−343]=[320018],观察得AAT[10]=32[10], AAT[01]=18[01]。但是我们不能直接使用这一组特征向量,因为式子AV=UΣ明确告诉我们,一旦V确定下来,U也必须取能够满足该式的向量,所以此处Av2=[0−√18]=u2σ2=[0−1]√18,则u1=[10], u2=[0−1]。(这个问题在本讲的官方笔记中有详细说明。)
-
补充:AB的特征值与BA的特征值相同,证明来自Are the eigenvalues of AB equal to the eigenvalues of BA? (Citation needed!):
取λ≠0,v是AB在特征值取λ时的的特征向量,则有Bv≠0,并有λBv=B(λv)=B(ABv)=(BA)Bv,所以Bv是BA在特征值取同一个λ时的特征向量。
再取AB的特征值λ=0,则0=detAB=detAdetB=detBA,所以λ=0也是BA的特征值,得证。
最终,我们得到[44−33]=[100−1][√3200√18][1√21√21√2−1√2]。
再做一个例子,A=[4386],这是个秩一矩阵,有零空间。A的行空间为[43]的倍数,A的列空间为[48]的倍数。
- 标准化向量得v1=[0.80.6], u1=1√5[12]。
- ATA=[4836][4386]=[80606045],由于A是秩一矩阵,则ATA也不满秩,所以必有特征值0,则另特征值一个由迹可知为125。
- 继续求零空间的特征向量,有v2=[0.6−0,8], u1=1√5[2−1]
最终得到[4386]=[12_2−1_][√125000_][0.80.60.6_−0.8_],其中下划线部分都是与零空间相关的部分。
- v1, ⋯, vr是行空间的标准正交基;
- u1, ⋯, ur是列空间的标准正交基;
- vr+1, ⋯, vn是零空间的标准正交基;
- ur+1, ⋯, um是左零空间的标准正交基。
通过将矩阵写为Avi=σiui形式,将矩阵对角化,向量u, v之间没有耦合,A乘以每个v都能得到一个相应的u。
第三十一讲:线性变换及对应矩阵
如何判断一个操作是不是线性变换?线性变换需满足以下两个要求:
即变换T需要同时满足加法和数乘不变的性质。将两个性质合成一个式子为:T(cv+dw)=cT(v)+dT(w)
例1,二维空间中的投影操作,T:R2→R2,它可以将某向量投影在一条特定直线上。检查一下投影操作,如果我们将向量长度翻倍,则其投影也翻倍;两向量相加后做投影与两向量做投影再相加结果一致。所以投影操作是线性变换。
“坏”例1,二维空间的平移操作,即平面平移:
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
fig = plt.figure()
sp1 = plt.subplot(221)
vectors_1 = np.array([[0,0,3,2],])
X_1, Y_1, U_1, V_1 = zip(*vectors_1)
plt.axhline(y=0, c='black')
plt.axvline(x=0, c='black')
sp1.quiver(X_1, Y_1, U_1, V_1, angles='xy', scale_units='xy', scale=1)
sp1.set_xlim(0, 10)
sp1.set_ylim(0, 5)
sp1.set_xlabel("before shifted")
sp2 = plt.subplot(222)
vector_2 = np.array([[0,0,3,2],
[3,2,2,0],
[0,0,5,2],
[0,0,10,4]])
X_2,Y_2,U_2,V_2 = zip(*vector_2)
plt.axhline(y=0, c='black')
plt.axvline(x=0, c='black')
sp2.quiver(X_2, Y_2, U_2, V_2, angles='xy', scale_units='xy', scale=1)
sp2.set_xlim(0, 10)
sp2.set_ylim(0, 5)
sp2.set_xlabel("shifted by horizontal 2 then double")
sp3 = plt.subplot(223)
vectors_1 = np.array([[0,0,6,4],])
X_1, Y_1, U_1, V_1 = zip(*vectors_1)
plt.axhline(y=0, c='black')
plt.axvline(x=0, c='black')
sp3.quiver(X_1, Y_1, U_1, V_1, angles='xy', scale_units='xy', scale=1)
sp3.set_xlim(0, 10)
sp3.set_ylim(0, 5)
sp3.set_xlabel("double the vector")
sp4 = plt.subplot(224)
vector_2 = np.array([[0,0,6,4],
[6,4,2,0],
[0,0,8,4]])
X_2,Y_2,U_2,V_2 = zip(*vector_2)
plt.axhline(y=0, c='black')
plt.axvline(x=0, c='black')
sp4.quiver(X_2, Y_2, U_2, V_2, angles='xy', scale_units='xy', scale=1)
sp4.set_xlim(0, 10)
sp4.set_ylim(0, 5)
sp4.set_xlabel("doubled vector shifted by horizontal 2")
plt.subplots_adjust(hspace=0.33)
plt.draw()
plt.close(fig)
比如,上图中向量长度翻倍,再做平移,明显与向量平移后再翻倍的结果不一致。
有时我们也可以用一个简单的特例判断线性变换,检查T(0)?=0。零向量平移后结果并不为零。
所以平面平移操作并不是线性变换。
“坏”例2,求模运算,T(v)=‖v‖, T:R3→R1,这显然不是线性变换,比如如果我们将向量翻倍则其模翻倍,但如果我将向量翻倍取负,则其模依然翻倍。所以T(−v)≠−T(v)
例2,旋转45∘操作,T:R2→R2,也就是将平面内一个向量映射为平面内另一个向量。检查可知,如果向量翻倍,则旋转后同样翻倍;两个向量先旋转后相加,与这两个向量先相加后旋转得到的结果一样。
所以从上面的例子我们知道,投影与旋转都是线性变换。
例3,矩阵乘以向量,T(v)=Av,这也是一个(一系列)线性变换,不同的矩阵代表不同的线性变换。根据矩阵的运算法则有A(v+w)=A(v)+A(w), A(cv)=cAv。比如取A=[100−1],作用于平面上的向量v,会导致v的x分量不变,而y分量取反,也就是图像沿x轴翻转。
线性变换的核心,就是该变换使用的相应的矩阵。
比如我们需要做一个线性变换,将一个三维向量降至二维,T:R3→R2,则在T(v)=Av中,v∈R3, T(v)∈R2,所以A应当是一个2×3矩阵。
如果我们希望知道线性变换T对整个输入空间Rn的影响,我们可以找到空间的一组基v1, v2, ⋯, vn,检查T对每一个基的影响T(v1), T(v2), ⋯, T(vn),由于输入空间中的任意向量都满足:
所以我们可以根据T(v)推出线性变换T对空间内任意向量的影响,得到:
现在我们需要考虑,如何把一个与坐标无关的线性变换变成一个与坐标有关的矩阵呢?
在1式中,c1,c2,⋯,cn就是向量v在基v1,v2,⋯,vn上的坐标,比如分解向量v=[324]=3[100]+2[010]+4[001],式子将向量v分解在一组标准正交基[100],[010],[001]上。当然,我们也可以选用矩阵的特征向量作为基向量,基的选择是多种多样的。
我们打算构造一个矩阵A用以表示线性变换T:Rn→Rm。我们需要两组基,一组用以表示输入向量,一组用以表示输出向量。令v1,v2,⋯,vn为输入向量的基,这些向量来自Rn;w1,w2,⋯,wm作为输出向量的基,这些向量来自Rm。
我们用二维空间的投影矩阵作为例子:
fig = plt.figure()
vectors_1 = np.array([[0, 0, 3, 2],
[0, 0, -2, 3]])
X_1, Y_1, U_1, V_1 = zip(*vectors_1)
plt.axis('equal')
plt.axhline(y=0, c='black')
plt.axvline(x=0, c='black')
plt.quiver(X_1, Y_1, U_1, V_1, angles='xy', scale_units='xy', scale=1)
plt.plot([-6, 12], [-4, 8])
plt.annotate('$v_1=w_1$', xy=(1.5, 1), xytext=(10, -20), textcoords='offset points', size=14, arrowprops=dict(arrowstyle="->"))
plt.annotate('$v_2=w_2$', xy=(-1, 1.5), xytext=(-60, -20), textcoords='offset points', size=14, arrowprops=dict(arrowstyle="->"))
plt.annotate('project line', xy=(4.5, 3), xytext=(-90, 10), textcoords='offset points', size=14, arrowprops=dict(arrowstyle="->"))
ax = plt.gca()
ax.set_xlim(-5, 5)
ax.set_ylim(-4, 4)
ax.set_xlabel("Project Example")
plt.draw()
plt.close(fig)
从图中可以看到,设输入向量的基为v1,v2,v1就在投影上,而v2垂直于投影方向,输出向量的基为w1,w2,而v1=w1,v2=w2。那么如果输入向量为v=c1v1+c2v2,则输出向量为T(v)=c1v1,也就是线性变换去掉了法线方向的分量,输入坐标为(c1,c2),输出坐标变为(c1,0)。
找出这个矩阵并不困难,Av=w,则有[1000][c1c2]=[c10]。
本例中我们选取的基极为特殊,一个沿投影方向,另一个沿投影法线方向,其实这两个向量都是投影矩阵的特征向量,所以我们得到的线性变换矩阵是一个对角矩阵,这是一组很好的基。
所以,如果我们选取投影矩阵的特征向量作为基,则得到的线性变换矩阵将是一个包含投影矩阵特征值的对角矩阵。
继续这个例子,我们不再选取特征向量作为基,而使用标准基v1=[10],v2=[01],我们继续使用相同的基作为输出空间的基,即v1=w1,v2=w2。此时投影矩阵为P=aaTaTa=[12121212],这个矩阵明显没有上一个矩阵“好”,不过这个矩阵也是一个不错的对称矩阵。
总结通用的计算线性变换矩阵A的方法:
- 确定输入空间的基v1,v2,⋯,vn,确定输出空间的基w1,w2,⋯,wm;
- 计算T(v1)=a11w1+a21w2+⋯+am1wm,求出的系数ai1就是矩阵A的第一列;
- 继续计算T(v2)=a12w1+a22w2+⋯+am2wm,求出的系数ai2就是矩阵A的第二列;
- 以此类推计算剩余向量直到vn;
- 最终得到矩阵A=[a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮am1am2⋯amn]。
最后我们介绍一种不一样的线性变换,T=ddx:
-
设输入为c1+c2x+c3x3,基为1,x,x2;
-
则输出为导数:c2+2c3x,基为1,x;
所以我们需要求一个从三维输入空间到二维输出空间的线性变换,目的是求导。求导运算其实是线性变换,因此我们只要知道少量函数的求导法则(如sinx,cosx,ex),就能求出它们的线性组合的导数。
有A[c1c2c3]=[c22c3],从输入输出的空间维数可知,A是一个2×3矩阵,A=[010002]。
最后,矩阵的逆相当于对应线性变换的逆运算,矩阵的乘积相当于线性变换的乘积,实际上矩阵乘法也源于线性变换。
第三十二讲:基变换和图像压缩
图像压缩
本讲我们介绍一种图片有损压缩的一种方法:JPEG。
假设我们有一张图片,长宽皆为512个像素,我们用xi来表示第i个像素,如果是灰度照片,通常xi可以在[0,255]上取值,也就是8 bits。对于这承载这张图片信息的向量x来说,有x∈Rn,n=5122。而如果是彩色照片,通常需要三个量来表示一个像素,则向量长度也会变为现在的三倍。
如此大的数据不经过压缩很难广泛传播。教学录像采用的压缩方法就是JPEG(Joint Photographic Expert Group,联合图像专家组),该方法采用的就是基变换的方式压缩图像。比如说一块干净的黑白,其附近的像素值应该非常接近,此时如果一个像素一个像素的描述黑白灰度值就太浪费空间了,所以标准基在这种情况下并不能很好的利用图片的特性。
我们知道,标准基是 [10⋮0][01⋮0]⋯[00⋮1],我们想寻找一个更好的基。
我们试试使用别的基描述图片,比如:
- 基中含有的一个向量 [11⋯1]T,即分量全为1的向量,一个向量就可以完整的给出所有“像素一致图像”的信息;
- 另一个向量 [1−1⋯1−1]T,正负交替出现,比如描述国际象棋棋盘;
- 第三个个向量 [11⋯−1−1]T,一半正一半负,比如描述一半亮一半暗的图片;
傅里叶基
现在我们来介绍傅里叶基,以8×8傅里叶基为例(这表示我们每次只处理8×8像素的一小块图像):
Fn=[111⋯11ww2⋯wn−11w2w4⋯w2(n−1)⋮⋮⋮⋱⋮1wn−1w2(n−1)⋯w(n−1)2], w=ei2π/n, n=8,我们不需要深入8阶傅里叶基的细节,先看看使用傅里叶基的思路是怎样的。
每次处理8×8的一小块时,会遇到64个像素,也就是64个基向量,64个系数,在64维空间中利用傅里叶向量做基变换:
-
输入信号x为64维向量基变换→输出信号c为x在傅里叶基下的64个系数。
注意前面做的都是无损的步骤,我们只是选了R64的一组基,接着把信号用这组基表达出来。
接下来的步骤就涉及到压缩和损失了:
-
一种方法是扔掉较小的系数,这叫做阈值量化(thresholding),我们设定一个阈值,任何不在阈值范围内的基向量、系数都将丢弃,虽然有信息损失,但是只要阈值设置合理,肉眼几乎无法区别压缩前后的图片。经由此步处理,向量c变为ˆc,而ˆc将有很多0。
通常[11⋯1]T向量很难被丢弃,它通常具有较大的系数。但是[1−1⋯1−1]T向量在平滑信号中的可能性就很小了。前一个的向量称作低频信号,频率为0,后一个向量称作高频信号,也是我们能够得到的最高频率的信号,如果是噪音或抖动输出的就是它。
比如讲课的视频图像信号,这种平滑的情形下输出的大多是低频信号,很少出现噪音。
-
接着我们用这些系数ˆc来重构信号,用这些系数乘以对应的基向量ˆx=∑ˆcivi,但是这个求和不再是64项求和了,因为压缩后的系数中有很多零存在,比如说我们压缩后ˆc中仅有三个非零项,那么压缩比将近达到21:1。
我们再来提一下视频压缩:视频是一系列连续图像,且相近的帧非常接近,而我们的压缩算法就需要利用这个相近性质。在实际生活中,从时间与空间的角度讲,事物不会瞬间改变。
小波基
接下来介绍另一组基,它是傅里叶基的竞争对手,名为小波(wavelets),同样以8×8为例:
[11111111][1111−1−1−1−1][11−1−10000][000011−1−1][1−1000000][001−10000][00001−100][0000001−1]。
可以看出傅里叶基中频率最高的向量为小波后四个基向量之和。
在标准基下的一组(按八个一组计算,P∈R8)像素P=[p1p2⋮p8]=c1w1+c2w2+⋯+cnwn=[w1 w2 ⋯ wn][c1c2⋮cn],即P=WC,我们需要计算像素向量在另一组基下系数,所以有C=W−1P。
此时我们发现,如果选取“好的基”会使得逆矩阵的求解过程变简单,所谓“好的基”:
-
计算快;
我们需要大量使用P=WC来计算整幅图在另一个基下的表达,在傅里叶变换中我们学习了快速傅里叶变换(FFT),同样的在小波变换中也有快速小波变换(FWT);
另外的,我们需要计算其逆矩阵,所以这个逆矩阵计算也必须快,观察小波基不难发现基向量相互正交,假设我们已经对小波基做了标准化处理,则小波基是一组标准正交基,所以有W−1=WT。
-
仅需少量向量即可最大限度的重现图像。
因为在图像压缩时,我们会舍弃较小的系数,比如c5,c6,c7,c8,所以后四个的基向量都会被舍弃,重现图像时仅使用前四个基向量的线性组合,而好的基选取会在使用较少基的前提下保证图像质量不会有较大损失。
题外话:JPEG2000标准会将小波基纳入压缩算法。我们上面介绍的是最简单的一组小波基,而FBI的指纹识别或JPEG2000的压缩算法纳入的是更加平滑的小波基,不会使用像上面介绍的那种直接从1变为−1的基。
要想继续了解小波基,可以参考一篇非常精彩的文章能不能通俗的讲解下傅立叶分析和小波分析之间的关系?——“咚懂咚懂咚“的答案
基变换
前面介绍小波基的时候我们就已经做了一次基变换。
将目标基的向量按列组成矩阵W,则基变换就是[x]x=Wc→[c]。
看一个例子,有线性变换T:R8→R8,在第一组基v1,v2,⋯,v8上计算得到矩阵A,在第二组基w1,w2,⋯,wn上计算得到矩阵B。先说结论,矩阵A,B是相似的,也就是有B=M−1AM,而M就是基变换矩阵。
进行基变换时会发生两件事:
-
每个向量都会有一组新的坐标,而x=Wc就是新旧坐标的关系;
-
每个线性变换都会有一个新的矩阵,而B=M−1AM就是新旧矩阵的关系。
再来看什么是A矩阵?
对于第一组基v1,v2,⋯,v8,要完全了解线性变换T,只需要知道T作用在基的每一个向量上会产生什么结果即可。因为在这个基下的每一个向量都可以写成x=c1v1+c2v2+⋯+c8v8的形式,所以T(x)=c1T(v1)+c2T(v2)+⋯+c8T(v8)。
而且T(v1)=a11v1+a21v2+⋯+a81v8, T(v2)=a12v1+a22v2+⋯+a82v8, ⋯,则矩阵[A]=[a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮am1am2⋯amn]
这些都是上一讲结尾所涉及的知识。
最后我们以一个更加特殊的基收场,设v1,v2,⋯,vn是一组特征向量,也就是T(vi)=λ1vi,那么问题就是矩阵A是什么?
继续使用线性变换中学到的,输入的第一个向量v1经由T加工后得到λ1v1,第二个向量v2T→λ2v2,继续做下去,最终有vn=vnT→λnvn。除了λivi外的其他基向量都变为0,那么矩阵A=[λ1λ2⋱λn]。
这是一个非常完美的基,我们在图像处理中最想要的就是这种基,但是找出像素矩阵的特征向量代价太大,所以我们找了一些代价小同时效果也不错的基,比如小波基、傅里叶基等等。
第三十三讲:单元检测3复习
在上一次复习中,我们已经涉及了求特征值与特征向量(通过解方程det(A−λI)=0得出λ,再将λ带入A−λI求其零空间得到x)。
接下的章节来我们学习了:
- 解微分方程dudt=Au,并介绍了指数矩阵eAt;
- 介绍了对称矩阵的性质A=AT,了解了其特征值均为实数且总是存在足量的特征向量(即使特征值重复特征向量也不会短缺,总是可以对角化);同时对称矩阵的特征向量正交,所以对称矩阵对角化的结果可以表示为A=QΛQT;
- 接着我们学习了正定矩阵;
- 然后学习了相似矩阵,B=M−1AM,矩阵A,B特征值相同,其实相似矩阵是用不同的基表示相同的东西;
- 最后我们学习了奇异值分解A=UΣVT。
现在,我们继续通过例题复习这些知识点。
-
解方程dudt=Au=[0−1010−1010]u。
首先通过A的特征值/向量求通解u(t)=c1eλ1tx1+c2eλ2tx2+c3eλ3tx3,很明显矩阵是奇异的,所以有λ1=0;
继续观察矩阵会发现AT=−A,这是一个反对称矩阵(anti-symmetric)或斜对陈矩阵(skew-symmetric),这与我们在第二十一讲介绍过的旋转矩阵类似,它的特征值应该为纯虚数(特征值在虚轴上),所以我们猜测其特征值应为0⋅i, b⋅i, −b⋅i。通过解det(A−λI)=0验证一下:[−λ−101−λ−101λ]=λ3+2λ=0,λ2=√2i,λ3=−√2i。
此时u(t)=c1+c2e√2itx2+c3e−√2itx3,e√2it始终在复平面单位圆上,所以u(t)及不发散也不收敛,它只是具有周期性。当t=0时有u(0)=c1+c2+c3,如果使e√2iT=1即√2iT=2πi则也能得到u(T)=c1+c2+c3,周期T=π√2。
另外,反对称矩阵同对称矩阵一样,具有正交的特征向量。当矩阵满足什么条件时,其特征向量相互正交?答案是必须满足AAT=ATA。所以对称矩阵A=AT满足此条件,同时反对称矩阵A=−AT也满足此条件,而正交矩阵Q−1=QT同样满足此条件,这三种矩阵的特征向量都是相互正交的。
上面的解法并没有求特征向量,进而通过u(t)=eAtu(0)得到通解,现在我们就来使用指数矩阵来接方程。如果矩阵可以对角化(在本例中显然可以),则A=SΛS−1,eAt=SeΛtS−1=S[eλ1teλ1t⋱eλ1t]S−1,这个公式在能够快速计算S,λ时很方便求解。
-
已知矩阵的特征值λ1=0,λ2=c,λ3=2,特征向量x1=[111],x2=[1−10],x3=[11−2]:
c如何取值才能保证矩阵可以对角化?其实可对角化只需要有足够的特征向量即可,而现在特征向量已经足够,所以c可以取任意值。
c如何取值才能保证矩阵对称?我们知道,对称矩阵的特征值均为实数,且注意到给出的特征向量是正交的,有了实特征值及正交特征向量,我们就可以得到对称矩阵。
c如何取值才能使得矩阵正定?已经有一个零特征值了,所以矩阵不可能是正定的,但可以是半正定的,如果c去非负实数。
c如何取值才能使得矩阵是一个马尔科夫矩阵?在第二十四讲我们知道马尔科夫矩阵的性质:必有特征值等于1,其余特征值均小于1,所以A不可能是马尔科夫矩阵。
c取何值才能使得P=A2是一个投影矩阵?我们知道投影矩阵的一个重要性质是P2=P,所以有对其特征值有λ2=λ,则c=0,2。
题设中的正交特征向量意义重大,如果没有正交这个条件,则矩阵A不会是对称、正定、投影矩阵。因为特征向量的正交性我们才能直接去看特征值的性质。
-
复习奇异值分解,A=UΣVT:
先求正交矩阵V:ATA=VΣTUTUΣVT=V(ΣTΣ)VT,所以V是矩阵ATA的特征向量矩阵,而矩阵ΣTΣ是矩阵ATA的特征值矩阵,即ATA的特征值为σ2。
接下来应该求正交矩阵U:AAT=UΣTVTVΣUT=U(ΣTΣ)UT,但是请注意,我们在这个式子中无法确定特征向量的符号,我们需要使用Avi=σiui,通过已经求出的vi来确定ui的符号(因为AV=UΣ),进而求出U。
已知A=[u1 u2][3002][v1 v2]T
从已知的Σ矩阵可以看出,A矩阵是非奇异矩阵,因为它没有零奇异值。另外,如果把Σ矩阵中的2改成−5,则题目就不再是奇异值分解了,因为奇异值不可能为负;如果将2变为0,则A是奇异矩阵,它的秩为1,零空间为1维,v2在其零空间中。
-
A是正交对称矩阵,那么它的特征值具有什么特点?
首先,对于对称矩阵,有特征值均为实数;然后是正交矩阵,直觉告诉我们|λ|=1。来证明一下,对于Qx=λx,我们两边同时取模有‖Qx‖=|λ|‖x‖,而正交矩阵不会改变向量长度,所以有‖x‖=|λ|‖x‖,因此λ=±1。
A是正定的吗?并不一定,因为特征向量可以取−1。
A的特征值没有重复吗?不是,如果矩阵大于2阶则必定有重复特征值,因为只能取±1。
A可以被对角化吗?是的,任何对称矩阵、任何正交矩阵都可以被对角化。
A是非奇异矩阵吗?是的,正交矩阵都是非奇异矩阵。很明显它的特征值都不为零。
证明P=12(A+I)是投影矩阵。
我们使用投影矩阵的性质验证,首先由于A是对称矩阵,则P一定是对称矩阵;接下来需要验证P2=P,也就是14(A2+2A+I)=12(A+I)。来看看A2是什么,A是正交矩阵则AT=A−1,而A又是对称矩阵则A=AT=A−1,所以A2=I。带入原式有14(2A+2I)=12(A+I),得证。
我们可以使用特征值验证,A的特征值可以取±1,则A+I的特征值可以取0,2,12(A+I)的特征值为0,1,特征值满足投影矩阵且它又是对称矩阵,得证。
第三十四讲:左右逆和伪逆
前面我们涉及到的逆(inverse)都是指左、右乘均成立的逆矩阵,即A−1A=I=AA−1。在这种情况下,m×n矩阵A满足m=n=rank(A),也就是满秩方阵。
左逆(left inserve)
记得我们在最小二乘一讲(第十六讲)介绍过列满秩的情况,也就是列向量线性无关,但行向量通常不是线性无关的。常见的列满秩矩阵A满足m>n=rank(A)。
列满秩时,列向量线性无关,所以其零空间中只有零解,方程Ax=b可能有一个唯一解(b在A的列空间中,此特解就是全部解,因为通常的特解可以通过零空间中的向量扩展出一组解集,而此时零空间只有列向量),也可能无解(b不在A的列空间中)。
另外,此时行空间为Rn,也正印证了与行空间互为正交补的零空间中只有列向量。
现在来观察ATA,也就是在m>n=rank(A)的情况下,n×m矩阵乘以m×n矩阵,结果为一个满秩的n×n矩阵,所以ATA是一个可逆矩阵。也就是说(ATA)−1AT⏟A=I成立,而大括号部分的(ATA)−1AT称为长方形矩阵A的左逆
顺便复习一下最小二乘一讲,通过关键方程ATAˆx=ATb,A−1left被当做一个系数矩阵乘在b向量上,求得b向量投影在A的列空间之后的解ˆx=(ATA)−1ATb。如果我们强行给左逆左乘矩阵A,得到的矩阵就是投影矩阵P=A(ATA)−1AT,来自p=Aˆx=A(ATA)−1AT,它将右乘的向量b投影在矩阵A的列空间中。
再来观察AAT矩阵,这是一个m×m矩阵,秩为rank(AAT)=n<m,也就是说AAT是不可逆的,那么接下来我们看看右逆。
右逆(right inverse)
可以与左逆对称的看,右逆也就是研究m×n矩阵A行满秩的情况,此时n>m=rank(A)。对称的,其左零空间中仅有零向量,即没有行向量的线性组合能够得到零向量。
行满秩时,矩阵的列空间将充满向量空间C(A)=Rm,所以方程Ax=b总是有解集,由于消元后有n−m个自由变量,所以方程的零空间为n−m维。
与左逆对称,再来观察AAT,在n>m=rank(A)的情况下,m×n矩阵乘以n×m矩阵,结果为一个满秩的m×m矩阵,所以此时AAT是一个满秩矩阵,也就是AAT可逆。所以AAT(AAT)⏟=I,大括号部分的AT(AAT)称为长方形矩阵的右逆
同样的,如果我们强行给右逆右乘矩阵A,将得到另一个投影矩阵P=AT(AAT)A,与上一个投影矩阵不同的是,这个矩阵的A全部变为AT了。所以这是一个能够将右乘的向量b投影在A的行空间中。
前面我们提及了逆(方阵满秩),并讨论了左逆(矩阵列满秩)、右逆(矩阵行满秩),现在看一下第四种情况,m×n矩阵A不满秩的情况。
伪逆(pseudo inverse)
有m×n矩阵A,满足rank(A)<min(m, n),则
- 列空间C(A)∈Rm, dimC(A)=r,左零空间N(AT)∈Rm, dimN(AT)=m−r,列空间与左零空间互为正交补;
- 行空间C(AT)∈Rn, dimC(AT)=r,零空间N(A)∈Rn, dimN(A)=n−r,行空间与零空间互为正交补。
现在任取一个向量x,乘上A后结果Ax一定落在矩阵A的列空间C(A)中。而根据维数,x∈Rn, Ax∈Rm,那么我们现在猜测,输入向量x全部来自矩阵的行空间,而输出向量Ax全部来自矩阵的列空间,并且是一一对应的关系,也就是Rn的r维子空间到Rm的r维子空间的映射。
而矩阵A现在有这些零空间存在,其作用是将某些向量变为零向量,这样Rn空间的所有向量都包含在行空间与零空间中,所有向量都能由行空间的分量和零空间的分量构成,变换将零空间的分量消除。但如果我们只看行空间中的向量,那就全部变换到列空间中了。
那么,我们现在只看行空间与列空间,在行空间中任取两个向量x, y∈C(AT),则有Ax≠Ay。所以从行空间到列空间,变换A是个不错的映射,如果限制在这两个空间上,A可以说“是个可逆矩阵”,那么它的逆就称作伪逆,而这个伪逆的作用就是将列空间的向量一一映射到行空间中。通常,伪逆记作A+,因此Ax=(Ax), y=A+(Ay)。
现在我们来证明对于x,y∈C(AT), x≠y,有Ax,Ay∈C(A), Ax≠Ay:
- 反证法,设Ax=Ay,则有A(x−y)=0,即向量x−y∈N(A);
- 另一方面,向量x,y∈C(AT),所以两者之差x−y向量也在C(AT)中,即x−y∈C(AT);
- 此时满足这两个结论要求的仅有一个向量,即零向量同时属于这两个正交的向量空间,从而得到x=y,与题设中的条件矛盾,得证。
伪逆在统计学中非常有用,以前我们做最小二乘需要矩阵列满秩这一条件,只有矩阵列满秩才能保证ATA是可逆矩阵,而统计中经常出现重复测试,会导致列向量线性相关,在这种情况下ATA就成了奇异矩阵,这时候就需要伪逆。
接下来我们介绍如何计算伪逆A+:
其中一种方法是使用奇异值分解,A=UΣVT,其中的对角矩阵型为Σ=[σ1⋱σ2[0]],对角线非零的部分来自ATA, AAT比较好的部分,剩下的来自左/零空间。
我们先来看一下Σ矩阵的伪逆是多少,这是一个m×n矩阵,rank(Σ)=r,Σ+=[1σ1⋱1σr[0]],伪逆与原矩阵有个小区别:这是一个n×m矩阵。则有ΣΣ+=[1⋱1[0]]m×m,Σ+Σ=[1⋱1[0]]n×n。
观察ΣΣ+和Σ+Σ不难发现,ΣΣ+是将向量投影到列空间上的投影矩阵,而Σ+Σ是将向量投影到行空间上的投影矩阵。我们不论是左乘还是右乘伪逆,得到的不是单位矩阵,而是投影矩阵,该投影将向量带入比较好的空间(行空间和列空间,而不是左/零空间)。
接下来我们来求A的伪逆:
第三十五讲:期末复习
依然是从以往的试题入手复习知识点。
-
已知m×n矩阵A,有Ax=[100]无解;Ax=[010]仅有唯一解,求关于m,n,rank(A)的信息。
首先,最容易判断的是m=3;而根据第一个条件可知,矩阵不满秩,有r<m;根据第二个条件可知,零空间仅有零向量,也就是矩阵消元后没有自由变量,列向量线性无关,所以有r=n。
综上,有m=3>n=r。
根据所求写出一个矩阵A的特例:A=[001001]。
detATA?=detAAT:不相等,因为ATA可逆而AAT不可逆,所以行列式不相等。(但是对于方阵,detAB=detBA恒成立。)
ATA可逆吗?是,因为r=n,矩阵列向量线性无关,即列满秩。
AAT正定吗?否,因为AAT是3×n矩阵与n×3矩阵之积,是一个三阶方阵,而AAT秩为2,所以不是正定矩阵。(不过AAT一定是半正定矩阵。)
求证ATy=c至少有一个解:因为A的列向量线性无关,所以AT的行向量线性无关,消元后每行都有主元,且总有自由变量,所以零空间中有非零向量,零空间维数是m−r(可以直接从dimN(AT)=m−r得到结论)。
-
设A=[v1 v2 v3],对于Ax=v1−v2+v3,求x。
按列计算矩阵相乘,有x=[1−11]。
若Ax=v_1-v_2+v_3=0,则解是唯一的吗?为什么。如果解释唯一的,则零空间中只有零向量,而在此例中x=[1−11]就在零空间中,所以解不唯一。
若v1,v2,v3是标准正交向量,那么怎样的线性组合c1v1+c2v2能够最接近v3?此问是考察投影概念,由于是正交向量,所以只有0向量最接近v3。
-
矩阵A=[.2.4.3.4.2.3.4.4.4],求稳态。
这是个马尔科夫矩阵,前两之和为第三列的两倍,奇异矩阵总有一个特征值为0,而马尔科夫矩阵总有一个特征值为1,剩下一个特征值从矩阵的迹得知为−.2。
再看马尔科夫过程,设从u(0)开始,uk=Aku0,u0=[0100]。先代入特征值λ1=0, λ2=1, λ3=−.2查看稳态uk=c1λk1x1+c2λk2x2+c3λk3x3,当k→∞,第一项与第三项都会消失,剩下u∞=c2x2。
到这里我们只需求出λ2对应的特征向量即可,带入特征值求解(A−I)x=0,有[−.8.4.3.4−.8.3.4.4−.6][???]=[000],可以消元得,也可以直接观察得到x2=[334]。
剩下就是求c2了,可以通过u0一一解出每个系数,但是这就需要解出每一个特征值。另一种方法,我们可以通过马尔科夫矩阵的特性知道,对于马尔科夫过程的每一个uk都有其分量之和与初始值分量之和相等,所以对于x2=[334],有c2=1。所以最终结果是u∞=[334]。
-
对于二阶方阵,回答以下问题:
求投影在直线a=[4−3]上的投影矩阵:应为P=aaTaTa。
已知特征值λ1=2, x1=[12]λ2=3, x2=[21]求原矩阵A:从对角化公式得A=SΛS−1=[1221][0003][1221]−1,解之即可。
A是一个实矩阵,且对任意矩阵B,A都不能分解成A=BTB,给出A的一个例子:我们知道BTB是对称的,所以给出一个非对称矩阵即可。
矩阵A有正交的特征向量,但不是对称的,给出一个A的例子:我们在三十三讲提到过,反对称矩阵,因为满足AAT=ATA而同样具有正交的特征向量,所以有A=[01−10]或旋转矩阵[cosθ−sinθsinθcosθ],这些矩阵都具有复数域上的正交特征向量组。 -
最小二乘问题,因为时间的关系直接写出计算式和答案,[101112][CD]=[341](Ax=b),解得[ˆCˆD]=[113−1]。
求投影后的向量p:向量p就是向量b在矩阵A列空间中的投影,所以p=[p1p2p3]=[101112][ˆCˆD]。
求拟合直线的图像:x=0,1,2时y=p1,p2,p2所在的直线的图像,y=ˆC+ˆDx即y=113−x。
%matplotlib inline
import matplotlib.pyplot as plt
from sklearn import linear_model
import numpy as np
import pandas as pd
import seaborn as sns
x = np.array([0, 1, 2]).reshape((-1,1))
y = np.array([3, 4, 1]).reshape((-1,1))
predict_line = np.array([-1, 4]).reshape((-1,1))
regr = linear_model.LinearRegression()
regr.fit(x, y)
ey = regr.predict(x)
fig = plt.figure()
plt.axis('equal')
plt.axhline(y=0, c='black')
plt.axvline(x=0, c='black')
plt.scatter(x, y, c='r')
plt.scatter(x, regr.predict(x), s=20, c='b')
plt.plot(predict_line, regr.predict(predict_line), c='g', lw='1')
[ plt.plot([x[i], x[i]], [y[i], ey[i]], 'r', lw='1') for i in range(len(x))]
plt.draw()
plt.close(fig)
-
接上面的题目
求一个向量b使得最小二乘求得的[ˆCˆD]=[00]:我们知道最小二乘求出的向量[ˆCˆD]使得A列向量的线性组合最接近b向量(即b在A列空间中的投影),如果这个线性组合为0向量(即投影为0),则b向量与A的列空间正交,所以可以取b=[1−21]同时正交于A的两个列向量。
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 震惊!C++程序真的从main开始吗?99%的程序员都答错了
· 单元测试从入门到精通
· 【硬核科普】Trae如何「偷看」你的代码?零基础破解AI编程运行原理
· 上周热点回顾(3.3-3.9)
· winform 绘制太阳,地球,月球 运作规律