其他分享
首页 > 其他分享> > 强撸MIT18-06灰飞烟灭(二)

强撸MIT18-06灰飞烟灭(二)

作者:互联网

第十九讲:行列式公式和代数余子式

上一讲中,我们从三个简单的性质扩展出了一些很好的推论,本讲将继续使用这三条基本性质:

  1. \(\det I=1\);
  2. 交换行行列式变号;
  3. 对行列式的每一行都可以单独使用线性运算,其值不变;

我们使用这三条性质推导二阶方阵行列式:

\[\begin{vmatrix}a&b\\c&d\end{vmatrix}=\begin{vmatrix}a&0\\c&d\end{vmatrix}+\begin{vmatrix}0&b\\c&d\end{vmatrix}=\begin{vmatrix}a&0\\c&0\end{vmatrix}+\begin{vmatrix}a&0\\0&d\end{vmatrix}+\begin{vmatrix}0&b\\c&0\end{vmatrix}+\begin{vmatrix}0&b\\0&d\end{vmatrix}=ad-bc \]

按照这个方法,我们继续计算三阶方阵的行列式,可以想到,我们保持第二、三行不变,将第一行拆分为个行列式之和,再将每一部分的第二行拆分为三部分,这样就得到九个行列式,再接着拆分这九个行列式的第三行,最终得到二十七个行列式。可以想象到,这些矩阵中有很多值为零的行列式,我们只需要找到不为零的行列式,求和即可。

\[\begin{vmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\end{vmatrix}=\begin{vmatrix}a_{11}&0&0\\0&a_{22}&0\\0&0&a_{33}\end{vmatrix}+\begin{vmatrix}a_{11}&0&0\\0&0&a_{23}\\0&a_{32}&0\end{vmatrix}+\begin{vmatrix}0&a_{12}&0\\a_{21}&0&0\\0&0&a_{33}\end{vmatrix}+\begin{vmatrix}0&a_{12}&0\\0&0&a_{23}\\a_{31}&0&0\end{vmatrix}+\begin{vmatrix}0&0&a_{13}\\a_{21}&0&0\\0&a_{32}&0\end{vmatrix}+\begin{vmatrix}0&0&a_{13}\\0&a_{22}&0\\a_{31}&0&0\end{vmatrix} \]

\[原式=a_{11}a_{22}a_{33}-a_{11}a_{23}a_{32}-a_{12}a_{21}a_{33}+a_{12}a_{23}a_{31}+a_{13}a_{21}a_{32}-a_{13}a_{22}a_{31}\tag{1} \]

同理,我们想继续推导出阶数更高的式子,按照上面的式子可知\(n\)阶行列式应该可以分解成\(n!\)个非零行列式(占据第一行的元素有\(n\)种选择,占据第二行的元素有\(n-1\)种选择,以此类推得\(n!\)):

\[\det A=\sum_{n!} \pm a_{1\alpha}a_{2\beta}a_{3\gamma}\cdots a_{n\omega}, (\alpha, \beta, \gamma, \omega)=P_n^n\tag{2} \]

这个公式还不完全,接下来需要考虑如何确定符号:

\[\begin{vmatrix}0&0&\overline 1&\underline 1\\0&\overline 1&\underline 1&0\\\overline 1&\underline 1&0&0\\\underline 1&0&0&\overline 1\end{vmatrix} \]

此处引入代数余子式(cofactor)的概念,它的作用是把\(n\)阶行列式化简为\(n-1\)阶行列式。

于是我们把\((1)\)式改写为:

\[a_{11}(a_{22}a_{33}-a_{23}a_{32})+a_{12}(a_{21}a_{33}-a_{23}a_{31})+a_{13}(a_{21}a_{32}-a_{22}a_{31}) \]

\[\begin{vmatrix}a_{11}&0&0\\0&a_{22}&a_{23}\\0&a_{32}&a_{33}\end{vmatrix}+\begin{vmatrix}0&a_{12}&0\\a_{21}&0&a_{23}\\a_{31}&0&a_{33}\end{vmatrix}+\begin{vmatrix}0&0&a_{13}\\a_{21}&a_{22}&0\\a_{31}&a_{32}&0\end{vmatrix} \]

于是,我们可以定义\(a_{ij}\)的代数余子式:将原行列式的第\(i\)行与第\(j\)列抹去后得到的\(n-1\)阶行列式记为\(C_{ij}\),\(i+j\)为偶时时取\(+\),\(i+j\)为奇时取\(-\)。

现在再来完善式子\((2)\):将行列式\(A\)沿第一行展开:

\[\det A=a_{11}C_{11}+a_{12}C_{12}+\cdots+a_{1n}C_{1n} \]

到现在为止,我们了解了三种求行列式的方法:

  1. 消元,\(\det A\)就是主元的乘积;
  2. 使用\((2)\)式展开,求\(n!\)项之积;
  3. 使用代数余子式。

计算例题:
\(A_4=\begin{vmatrix}1&1&0&0\\1&1&1&0\\0&1&1&1\\0&0&1&1\end{vmatrix}\stackrel{沿第一行展开}{=}\begin{vmatrix}1&1&0\\1&1&1\\0&1&1\end{vmatrix}-\begin{vmatrix}1&1&0\\0&1&1\\0&1&1\end{vmatrix}=-1-0=-1\)

第二十讲:克拉默法则、逆矩阵、体积

本讲主要介绍逆矩阵的应用。

求逆矩阵

我们从逆矩阵开始,对于二阶矩阵有\(\begin{bmatrix}a&b\\c&d\end{bmatrix}^{-1}=\frac{1}{ad-bc}\begin{bmatrix}d&-b\\-c&a\end{bmatrix}\)。观察易得,系数项就是行列式的倒数,而矩阵则是由一系列代数余子式组成的。先给出公式:

\[A^{-1}=\frac{1}{\det A}C^T \tag{1} \]

观察这个公式是如何运作的,化简公式得\(AC^T=(\det A)I\),写成矩阵形式有\(\begin{bmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\\vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{bmatrix}\begin{bmatrix}C_{11}&\cdots&C_{n1}\\C_{12}&\cdots&C_{n2}\\\vdots&\ddots&\vdots\\C_{1n}&\cdots&C_{nn}\end{bmatrix}=Res\)

对于这两个矩阵的乘积,观察其结果的元素\(Res_{11}=a_{11}C_{11}+a_{12}C_{12}+\cdots+a_{1n}C_{1n}\),这正是上一讲提到的将行列式按第一行展开的结果。同理,对\(Res_{22}, \cdots, Res_{nn}\)都有\(Res_{ii}=\det A\),即对角线元素均为\(\det A\)。

再来看非对角线元素:回顾二阶的情况,如果用第一行乘以第二行的代数余子式\(a_{11}C_{21}+a_{12}C_{22}\),得到\(a(-b)+ab=0\)。换一种角度看问题,\(a(-b)+ab=0\)也是一个矩阵的行列式值,即\(A_{s}=\begin{bmatrix}a&b\\a&b\end{bmatrix}\)。将\(\det A_{s}\)按第二行展开,也会得到\(\det A_{s}=a(-b)+ab\),因为行列式有两行相等所以行列式值为零。

推广到\(n\)阶,我们来看元素\(Res_{1n}=a_{11}C_{n1}+a_{12}C_{n2}+\cdots+a_{1n}C_{nn}\),该元素是第一行与最后一行的代数余子式相乘之积。这个式子也可以写成一个特殊矩阵的行列式,即矩阵\(A_{s}=\begin{bmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{n-a1}&a_{n-12}&\cdots&a_{n-1n}\\a_{11}&a_{12}&\cdots&a_{1n}\end{bmatrix}\)。计算此矩阵的行列式,将\(\det A_{s}\)按最后一行展开,也得到\(\det A_{s}=a_{11}C_{n1}+a_{12}C_{n2}+\cdots+a_{1n}C_{nn}\)。同理,行列式\(A_{s}\)有两行相等,其值为零。

结合对角线元素与非对角线元素的结果,我们得到\(Res=\begin{bmatrix}\det A&0&\cdots&0\\0&\det A&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\det A\end{bmatrix}\),也就是\((1)\)等式右边的\((\det A)I\),得证。

求解\(Ax=b\)

因为我们现在有了逆矩阵的计算公式,所以对\(Ax=b\)有\(x=A^{-1}b=\frac{1}{\det A}C^Tb\),这就是计算\(x\)的公式,即克莱默法则(Cramer's rule)。

现在来观察\(x=\frac{1}{\det A}C^Tb\),我们将得到的解拆分开来,对\(x\)的第一个分量有\(x_1=\frac{y_1}{\det A}\),这里\(y_1\)是一个数字,其值为\(y_1=b_1C_{11}+b_2C_{21}+\cdots+b_nC_{n1}\),每当我们看到数字与代数余子式乘之积求和时,都应该联想到求行列式,也就是说\(y_1\)可以看做是一个矩阵的行列式,我们设这个矩阵为\(B_1\)。所以有\(x_i=\frac{\det B_1}{\det A}\),同理有\(x_2=\frac{\det B_2}{\det A}\),\(x_2=\frac{\det B_2}{\det A}\)。

而\(B_1\)是一个型为\(\Bigg[b a_2 a_3 \cdots a_n\Bigg]\)的矩阵,即将矩阵\(A\)的第一列变为\(b\)向量而得到的新矩阵。其实很容易看出,\(\det B_1\)可以沿第一列展开得到\(y_1=b_1C_{11}+b_2C_{21}+\cdots+b_nC_{n1}\)。

一般的,有\(B_j=\Bigg[a_1 a_2 \cdots a_{j-1} b a_{j+1} \cdots a_n\Bigg]\),即将矩阵\(A\)的第\(j\)列变为\(b\)向量而得到的新矩阵。所以,对于解的分量有\(x_j=\frac{\det B_j}{\det A}\)。

这个公式虽然很漂亮,但是并不方便计算。

关于体积(Volume)

先提出命题:行列式的绝对值等于一个箱子的体积。

来看三维空间中的情形,对于\(3\)阶方阵\(A\),取第一行\((a_1,a_2,a_3)\),令其为三维空间中点\(A_1\)的坐标,同理有点\(A_2, A_3\)。连接这三个点与原点可以得到三条边,使用这三条边展开得到一个平行六面体,\(\left\|\det A\right\|\)就是该平行六面体的体积。

对于三阶单位矩阵,其体积为\(\det I=1\),此时这个箱子是一个单位立方体。这其实也证明了前面学过的行列式性质1。于是我们想,如果能接着证明性质2、3即可证明体积与行列式的关系。

对于行列式性质2,我们交换两行并不会改变箱子的大小,同时行列式的绝对值也没有改变,得证。

现在我们取矩阵\(A=Q\),而\(Q\)是一个标准正交矩阵,此时这个箱子是一个立方体,可以看出其实这个箱子就是刚才的单位立方体经过旋转得到的。对于标准正交矩阵,有\(Q^TQ=I\),等式两边取行列式得\(\det(Q^TQ)=1=\left|Q^T\right|\left|Q\right|\),而根据行列式性质10有\(\left|Q^T\right|=\left|Q\right|\),所以\(原式=\left|Q\right|^2=1, \left|Q\right|=\pm 1\)。

接下来在考虑不再是“单位”的立方体,即长方体。 假设\(Q\)矩阵的第一行翻倍得到新矩阵\(Q_2\),此时箱子变为在第一行方向上增加一倍的长方体箱子,也就是两个“标准正交箱子”在第一行方向上的堆叠。易知这个长方体箱子是原来体积的两倍,而根据行列式性质3.a有\(\det Q_2=\det Q\),于是体积也符合行列式的数乘性质。

我们来看二阶方阵的情形,\(\begin{vmatrix}a+a'&b+b'\\c&d\end{vmatrix}=\begin{vmatrix}a&b\\c&d\end{vmatrix}+\begin{vmatrix}a'&b'\\c&d\end{vmatrix}\)。在二阶情况中,行列式就是一个求平行四边形面积的公式,原来我们求由四个点\((0,0), (a,b), (c,d), (a+c,b+d)\)围成的四边形的面积,需要先求四边形的底边长,再做高求解,现在只需要计算\(\det A=ad-bc\)即可(更加常用的是求由\((0,0), (a,b), (c,d)\)围成的三角形的面积,即\(\frac{1}{2}ad-bc\))。也就是说,如果知道了歪箱子的顶点坐标,求面积(二阶情形)或体积(三阶情形)时,我们不再需要开方、求角度,只需要计算行列式的值就行了。

再多说两句我们通过好几讲得到的这个公式,在一般情形下,由点\((x_1,y_1), (x_2,y_2), (x_3,y_3)\)围成的三角形面积等于\(\frac{1}{2}\begin{vmatrix}x_1&y_1&1\\x_2&y_2&1\\x_3&y_3&1\end{vmatrix}\),计算时分别用第二行、第三行减去第一行化简到第三列只有一个\(1\)(这个操作实际作用是将三角形移动到原点),得到\(\frac{1}{2}\begin{vmatrix}x_1&y_1&1\\x_2-x_1&y_2-y_1&0\\x_3-x_1&y_3-y_1&0\end{vmatrix}\),再按照第三列展开,得到三角形面积等于\(\frac{(x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1)}{2}\)。

第二十一讲:特征值和特征向量

特征值、特征向量的由来

给定矩阵\(A\),矩阵\(A\)乘以向量\(x\),就像是使用矩阵\(A\)作用在向量\(x\)上,最后得到新的向量\(Ax\)。在这里,矩阵\(A\)就像是一个函数,接受一个向量\(x\)作为输入,给出向量\(Ax\)作为输出。

在这一过程中,我们对一些特殊的向量很感兴趣,他们在输入(\(x\))输出(\(Ax\))的过程中始终保持同一个方向,这是比较特殊的,因为在大多情况下,\(Ax\)与\(x\)指向不同的方向。在这种特殊的情况下,\(Ax\)平行于\(x\),我们把满足这个条件的\(x\)成为特征向量(Eigen vector)。这个平行条件用方程表示就是:

\[Ax=\lambda x\tag{1} \]

再提前透露一个特征值的性质:对于一个\(n\times n\)的矩阵,将会有\(n\)个特征值,而这些特征值的和与该矩阵对角线元素的和相同,因此我们把矩阵对角线元素称为矩阵的迹(trace)。$$\sum_{i=1}^n \lambda_i=\sum_{i=1}^n a_{ii}$$

在上面二阶转置矩阵的例子中,如果我们求得了一个特征值\(1\),那么利用迹的性质,我们就可以直接推出另一个特征值是\(-1\)。

求解\(Ax=\lambda x\)

对于方程\(Ax=\lambda x\),有两个未知数,我们需要利用一些技巧从这一个方程中一次解出两个未知数,先移项得\((A-\lambda I)x=0\)。

观察\((A-\lambda I)x=0\),右边的矩阵相当于将\(A\)矩阵平移了\(\lambda\)个单位,而如果方程有解,则这个平移后的矩阵\((A-\lambda I)\)一定是奇异矩阵。根据前面学到的行列式的性质,则有$$\det{(A-\lambda{I})}=0\tag{2}$$

这样一来,方程中就没有\(x\)了,这个方程也叫作特征方程(characteristic equation)。有了特征值,代回\((A-\lambda I)x=0\),继续求\((A-\lambda I)\)的零空间即可。

接下来,看一个关于特征向量认识的误区:已知\(Ax=\lambda x, Bx=\alpha x\),则有\((A+B)x=(\lambda+\alpha)x\),当\(B=3I\)时,在上例中我们看到,确实成立,但是如果\(B\)为任意矩阵,则推论不成立,因为这两个式子中的特征向量\(x\)并不一定相同,所以两个式子的通常情况是\(Ax=\lambda x, By=\alpha y\),它们也就无从相加了。

这一讲我们看到了足够多的“不好”的矩阵,下一讲会介绍一般情况下的特征值与特征向量。

第二十二讲:对角化和\(A\)的幂

对角化矩阵

上一讲我们提到关键方程\(Ax=\lambda x\),通过\(\det(A-\lambda I)=0\)得到特征向量\(\lambda\),再带回关键方程算出特征向量\(x\)。

在得到特征值与特征向量后,该如何使用它们?我们可以利用特征向量来对角化给定矩阵。

有矩阵\(A\),它的特征向量为\(x_1, x_2, \cdots, x_n\),使用特征向量作为列向量组成一个矩阵\(S=\Bigg[x_1x_2\cdots x_n\Bigg]\),即特征向量矩阵, 再使用公式$$S^{-1}AS=\Lambda\tag{1}$$将\(A\)对角化。注意到公式中有\(S^{-1}\),也就是说特征向量矩阵\(S\)必须是可逆的,于是我们需要\(n\)个线性无关的特征向量。

现在,假设\(A\)有\(n\)个线性无关的特征向量,将它们按列组成特征向量矩阵\(S\),则\(AS=A\Bigg[x_1x_2\cdots x_n\Bigg]\),当我们分开做矩阵与每一列相乘的运算时,易看出\(Ax_1\)就是矩阵与自己的特征向量相乘,其结果应该等于\(\lambda_1x_1\)。那么\(AS=\Bigg[(\lambda_1x_1)(\lambda_2x_2)\cdots(\lambda_nx_n)\Bigg]\)。可以进一步化简原式,使用右乘向量按列操作矩阵的方法,将特征值从矩阵中提出来,得到\(\Bigg[x_1x_2\cdots x_n\Bigg]\begin{bmatrix}\lambda_1&0&\cdots&0\\0&\lambda_2&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\lambda_n\end{bmatrix}=S\Lambda\)。

于是我们看到,从\(AS\)出发,得到了\(S\Lambda\),特征向量矩阵又一次出现了,后面接着的是一个对角矩阵,即特征值矩阵。这样,再继续左乘\(S^{-1}\)就得到了公式\((1)\)。当然,所以运算的前提条件是特征向量矩阵\(S\)可逆,即矩阵\(A\)有\(n\)个线性无关的特征向量。这个式子还要另一种写法,\(A=S\Lambda S^{-1}\)。

我们来看如何应用这个公式,比如说要计算\(A^2\)。

两种方法描述的是同一个现象,即对于矩阵幂运算\(A^2\),其特征向量不变,而特征值做同样的幂运算。对角矩阵\(\Lambda^2=\begin{bmatrix}\lambda_1^2&0&\cdots&0\\0&\lambda_2^2&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\lambda_n^2\end{bmatrix}\)。

特征值和特征向量给我们了一个深入理解矩阵幂运算的方法,\(A^k=S\Lambda^kS^{-1}\)。

再来看一个矩阵幂运算的应用:如果\(k\to\infty\),则\(A^k\to 0\)(趋于稳定)的条件是什么?从\(S\Lambda^kS^{-1}\)易得,\(|\lambda_i|<1\)。再次强调,所有运算的前提是矩阵\(A\)存在\(n\)个线性无关的特征向量。如果没有\(n\)个线性无关的特征向量,则矩阵就不能对角化。

关于矩阵可对角化的条件:

我们不打算深入研究有重复特征值的情形。

求\(u_{k+1}=Au_k\)

从\(u_1=Au_0\)开始,\(u_2=A^2u_0\),所有\(u_k=A^ku_0\)。下一讲涉及微分方程(differential equation),会有求导的内容,本讲先引入简单的差分方程(difference equation)。本例是一个一阶差分方程组(first order system)。

要解此方程,需要将\(u_0\)展开为矩阵\(A\)特征向量的线性组合,即\(u_0=c_1x_1+c_2x_2+\cdots+c_nx_n=\Bigg[x_1x_2\cdots x_n\Bigg]\begin{bmatrix}c_1\\c_2\\\vdots\\c_n\end{bmatrix}=Sc\)。于是\(Au_0=c_1Ax_1+c_2Ax_2+\cdots+c_nAx_n=c_1\lambda_1x_1+c_2\lambda_2x_2+\cdots+c_n\lambda_nx_n\)。继续化简原式,\(Au_0=\Bigg[x_1x_2\cdots x_n\Bigg]\begin{bmatrix}\lambda_1&0&\cdots&0\\0&\lambda_2&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\lambda_n\end{bmatrix}\begin{bmatrix}c_1\\c_2\\\vdots\\c_n\end{bmatrix}=S\Lambda c\)。用矩阵的方式同样可以得到该式:\(Au_0=S\Lambda S^{-1}u_0=S\Lambda S^{-1}Sc=S\Lambda c\)。

那么如果我们要求\(A^{100}u_0\),则只需要将\(\lambda\)变为\(\lambda^{100}\),而系数\(c\)与特征向量\(x\)均不变。

当我们真的要计算\(A^{100}u_0\)时,就可以使用\(S\Lambda^{100}c=c_1\lambda_1^{100}x_1+c_2\lambda_2^{100}x_2+\cdots+c_n\lambda_n^{100}x_n\)。

接下来看一个斐波那契数列(Fibonacci sequence)的例子:

\(0,1,1,2,3,5,8,13,\cdots,F_{100}=?\),我们要求第一百项的公式,并观察这个数列是如何增长的。可以想象这个数列并不是稳定数列,因此无论如何该矩阵的特征值并不都小于一,这样才能保持增长。而他的增长速度,则有特征值来决定。

已知\(F_{k+2}=F_{k_1}+F_{k}\),但这不是\(u_{k+1}=Au_{k}\)的形式,而且我们只要一个方程,而不是方程组,同时这是一个二阶差分方程(就像含有二阶导数的微分方程,希望能够化简为一阶倒数,也就是一阶差分)。

使用一个小技巧,令\(u_{k}=\begin{bmatrix}F_{k+1}\\F_{k}\end{bmatrix}\),再追加一个方程组成方程组:\(\begin{cases}F_{k+2}&=F_{k+1}+F_{k}\\F_{k+1}&=F_{k+1}\end{cases}\),再把方程组用矩阵表达得到\(\begin{bmatrix}F_{k+2}\\F_{k+1}\end{bmatrix}=\begin{bmatrix}1&1\\1&0\end{bmatrix}\begin{bmatrix}F_{k+1}\\F_{k}\end{bmatrix}\),于是我们得到了\(u_{k+1}=Au_{k}, A=\begin{bmatrix}1&1\\1&0\end{bmatrix}\)。我们把二阶标量方程(second-order scalar problem)转化为一阶向量方程组(first-order system)。

我们的矩阵\(A=\begin{bmatrix}1&1\\1&0\end{bmatrix}\)是一个对称矩阵,所以它的特征值将会是实数,且他的特征向量将会互相正交。因为是二阶,我们可以直接利用迹与行列式解方程组\(\begin{cases}\lambda_1+\lambda_2&=1\\\lambda_1\cdot\lambda_2&=-1\end{cases}\)。在求解之前,我们先写出一般解法并观察\(\left|A-\lambda I\right|=\begin{vmatrix}1-\lambda&1\\1&-\lambda\end{vmatrix}=\lambda^2-\lambda-1=0\),与前面斐波那契数列的递归式\(F_{k+2}=F_{k+1}+F_{k}\rightarrow F_{k+2}-F_{k+1}-F_{k}=0\)比较,我们发现这两个式子在项数与幂次上非常相近。

我们先来观察这个数列是如何增长的,数列增长由什么来控制?——特征值。哪一个特征值起决定性作用?——较大的一个。

\(F_{100}=c_1\left(\frac{1+\sqrt{5}}{2}\right)^{100}+c_2\left(\frac{1-\sqrt{5}}{2}\right)^{100}\approx c_1\left(\frac{1+\sqrt{5}}{2}\right)^{100}\),由于\(-0.618\)在幂增长中趋近于\(0\),所以近似的忽略该项,剩下较大的项,我们可以说数量增长的速度大约是\(1.618\)。可以看出,这种问题与求解\(Ax=b\)不同,这是一个动态的问题,\(A\)的幂在不停的增长,而问题的关键就是这些特征值。

最后,计算初始项\(u_0=\begin{bmatrix}F_1\\F_0\end{bmatrix}=\begin{bmatrix}1\\0\end{bmatrix}\),现在将初始项用特征向量表示出来\(\begin{bmatrix}1\\0\end{bmatrix}=c_1x_1+c_2x_2\),计算系数得\(c_1=\frac{\sqrt{5}}{5}, c_2=-\frac{\sqrt{5}}{5}\)。

来回顾整个问题,对于动态增长的一阶方程组,初始向量是\(u_0\),关键在于确定\(A\)的特征值及特征向量。特征值将决定增长的趋势,发散至无穷还是收敛于某个值。接下来需要找到一个展开式,把\(u_0\)展开成特征向量的线性组合。

下一讲将介绍求解微分方程。

第二十三讲:微分方程和\(e^{At}\)

微分方程\(\frac{\mathrm{d}u}{\mathrm{d}t}=Au\)

本讲主要讲解解一阶方程(first-order system)一阶倒数(first derivative)常系数(constant coefficient)线性方程,上一讲介绍了如何计算矩阵的幂,本讲将进一步涉及矩阵的指数形式。我们通过解一个例子来详细介绍计算方法。

有方程组\(\begin{cases}\frac{\mathrm{d}u_1}{\mathrm{d}t}&=-u_1+2u_2\\\frac{\mathrm{d}u_2}{\mathrm{d}t}&=u_1-2u_2\end{cases}\),则系数矩阵是\(A=\begin{bmatrix}-1&2\\1&-2\end{bmatrix}\),设初始条件为在\(0\)时刻\(u(0)=\begin{bmatrix}u_1\\u_2\end{bmatrix}=\begin{bmatrix}1\\0\end{bmatrix}\)。

稳定性:这个流动过程从\(u(0)=\begin{bmatrix}1\\0\end{bmatrix}\)开始,初始值\(1\)的一部分流入初始值\(0\)中,经过无限的时间最终达到稳态\(u(\infty)=\begin{bmatrix}\frac{2}{3}\\\frac{1}{3}\end{bmatrix}\)。所以,要使得\(u(t)\to 0\),则需要负的特征值。但如果特征值为复数呢?如\(\lambda=-3+6i\),我们来计算\(\left|e^{(-3+6i)t}\right|\),其中的\(\left|e^{6it}\right|\)部分为\(\left|\cos 6t+i\sin 6t\right|=1\),因为这部分的模为\(\cos^2\alpha+\sin^2\alpha=1\),这个虚部就在单位圆上转悠。所以只有实数部分才是重要的。所以我们可以把前面的结论改为需要实部为负数的特征值。实部会决定最终结果趋近于\(0\)或\(\infty\),虚部不过是一些小杂音。

收敛态:需要其中一个特征值实部为\(0\),而其他特征值的实部皆小于\(0\)。

发散态:如果某个特征值实部大于\(0\)。上面的例子中,如果将\(A\)变为\(-A\),特征值也会变号,结果发散。

再进一步,我们想知道如何从直接判断任意二阶矩阵的特征值是否均小于零。对于二阶矩阵\(A=\begin{bmatrix}a&b\\c&d\end{bmatrix}\),矩阵的迹为\(a+d=\lambda_1+\lambda_2\),如果矩阵稳定,则迹应为负数。但是这个条件还不够,有反例迹小于\(0\)依然发散:\(\begin{bmatrix}-2&0\\0&1\end{bmatrix}\),迹为\(-1\)但是仍然发散。还需要加上一个条件,因为\(\det A=\lambda_1\cdot\lambda_2\),所以还需要行列式为正数。

总结:原方程组有两个相互耦合的未知函数,\(u_1, u_2\)相互耦合,而特征值和特征向量的作则就是解耦,也就是对角化(diagonalize)。回到原方程组\(\frac{\mathrm{d}u}{\mathrm{d}t}=Au\),将\(u\)表示为特征向量的线性组合\(u=Sv\),代入原方程有\(S\frac{\mathrm{d}v}{\mathrm{d}t}=ASv\),两边同乘以\(S^{-1}\)得\(\frac{\mathrm{d}v}{\mathrm{d}t}=S^{-1}ASv=\Lambda v\)。以特征向量为基,将\(u\)表示为\(Sv\),得到关于\(v\)的对角化方程组,新方程组不存在耦合,此时\(\begin{cases}\frac{\mathrm{d}v_1}{\mathrm{d}t}&=\lambda_1v_1\\\frac{\mathrm{d}v_2}{\mathrm{d}t}&=\lambda_2v_2\\\vdots&\vdots\\\frac{\mathrm{d}v_n}{\mathrm{d}t}&=\lambda_nv_n\end{cases}\),这是一个各未知函数间没有联系的方程组,它们的解的一般形式为\(v(t)=e^{\Lambda t}v(0)\),则原方程组的解的一般形式为\(u(t)=e^{At}u(0)=Se^{\Lambda t}S^{-1}u(0)\)。这里引入了指数部分为矩阵的形式。

指数矩阵\(e^{At}\)

在上面的结论中,我们见到了\(e^{At}\)。这种指数部分带有矩阵的情况称为指数矩阵(exponential matrix)。

理解指数矩阵的关键在于,将指数形式展开称为幂基数形式,就像\(e^x=1+\frac{x^2}{2}+\frac{x^3}{6}+\cdots\)一样,将\(e^{At}\)展开成幂级数的形式为:

\[e^{At}=I+At+\frac{(At)^2}{2}+\frac{(At)^3}{6}+\cdots+\frac{(At)^n}{n!}+\cdots \]

再说些题外话,有两个极具美感的泰勒级数:\(e^x=\sum \frac{x^n}{n!}\)与\(\frac{1}{1-x}=\sum x^n\),如果把第二个泰勒级数写成指数矩阵形式,有\((I-At)^{-1}=I+At+(At)^2+(At)^3+\cdots\),这个式子在\(t\)非常小的时候,后面的高次项近似等于零,所以可以用来近似\(I-At\)的逆矩阵,通常近似为\(I+At\),当然也可以再加几项。第一个级数对我们而言比第二个级数好,因为第一个级数总会收敛于某个值,所以\(e^x\)总会有意义,而第二个级数需要\(A\)特征值的绝对值小于\(1\)(因为涉及矩阵的幂运算)。我们看到这些泰勒级数的公式对矩阵同样适用。

回到正题,我们需要证明\(Se^{\Lambda t}S^{-1}=e^{At}\),继续使用泰勒级数:

\[e^{At}=I+At+\frac{(At)^2}{2}+\frac{(At)^3}{6}+\cdots+\frac{(At)^n}{n!}+\cdots\\ e^{At}=SS^{-1}+S\Lambda S^{-1}t+\frac{S\Lambda^2S^{-1}}{2}t^2+\frac{S\Lambda^3S^{-1}}{6}t^3+\cdots+\frac{S\Lambda^nS^{-1}}{n!}t^n+\cdots\\ e^{At}=S\left(I+\Lambda t+\frac{\Lambda^2t^2}{2}+\frac{\Lambda^3t^3}{3}+\cdots+\frac{\Lambda^nt^n}{n}+\cdots\right)S^{-1}\\ e^{At}=Se^{\Lambda t}S^{-1} \]

需要注意的是,\(e^{At}\)的泰勒级数展开是恒成立的,但我们推出的版本却需要矩阵可对角化这个前提条件。

最后,我们来看看什么是\(e^{\Lambda t}\),我们将\(e^{At}\)变为对角矩阵就是因为对角矩阵简单、没有耦合,\(e^{\Lambda t}=\begin{bmatrix}e^{\lambda_1t}&0&\cdots&0\\0&e^{\lambda_2t}&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&e^{\lambda_nt}\end{bmatrix}\)。

有了\(u(t)=Se^{\Lambda t}S^{-1}u(0)\),再来看矩阵的稳定性可知,所有特征值的实部均为负数时矩阵收敛,此时对角线上的指数收敛为\(0\)。如果我们画出复平面,则要使微分方程存在稳定解,则特征值存在于复平面的左侧(即实部为负);要使矩阵的幂收敛于\(0\),则特征值存在于单位圆内部(即模小于\(1\)),这是幂稳定区域。(上一讲的差分方程需要计算矩阵的幂。)

同差分方程一样,我们来看二阶情况如何计算,有\(y''+by'+k=0\)。我们也模仿差分方程的情形,构造方程组\(\begin{cases}y''&=-by'-ky\\y'&=y'\end{cases}\),写成矩阵形式有\(\begin{bmatrix}y''\\y'\end{bmatrix}=\begin{bmatrix}-b&-k\\1&0\end{bmatrix}\begin{bmatrix}y'\\y\end{bmatrix}\),令\(u'=\begin{bmatrix}y''\\y'\end{bmatrix}, \ u=\begin{bmatrix}y'\\y\end{bmatrix}\)。

继续推广,对于\(5\)阶微分方程\(y'''''+by''''+cy'''+dy''+ey'+f=0\),则可以写作\(\begin{bmatrix}y'''''\\y''''\\y'''\\y''\\y'\end{bmatrix}=\begin{bmatrix}-b&-c&-d&-e&-f\\1&0&0&0&0\\0&1&0&0&0\\0&0&1&0&0\\0&0&0&1&0\end{bmatrix}\begin{bmatrix}y''''\\y'''\\y''\\y'\\y\end{bmatrix}\),这样我们就把一个五阶微分方程化为\(5\times 5\)一阶方程组了,然后就是求特征值、特征向量了步骤了。

第二十四讲:马尔科夫矩阵、傅里叶级数

马尔科夫矩阵

马尔科夫矩阵(Markov matrix)是指具有以下两个特性的矩阵:

  1. 矩阵中的所有元素大于等于\(0\);(因为马尔科夫矩阵与概率有关,而概率是非负的。)
  2. 每一列的元素之和为\(1\)

对于马尔科夫矩阵,我们关心幂运算过程中的稳态(steady state)。与上一讲不同,指数矩阵关系特征值是否为\(0\),而幂运算要达到稳态需要特征值为\(1\)。

根据上面两条性质,我们可以得出两个推论:

  1. 马尔科夫矩阵必有特征值为\(1\);
  2. 其他的特征值的绝对值皆小于\(1\)。

使用第二十二讲中得到的公式进行幂运算\(u_k=A^ku_0=S\Lambda^kS^{-1}u_0=S\Lambda^kS^{-1}Sc=S\Lambda^kc=c_1\lambda_1^kx_1+c_2\lambda_2^kx_2+\cdots+c_n\lambda_n^kx_n\),从这个公式很容易看出幂运算的稳态。比如我们取\(\lambda_1=1\),其他的特征值绝对值均小于\(1\),于是在经过\(k\)次迭代,随着时间的推移,其他项都趋近于\(0\),于是在\(k\to\infty\)时,有稳态\(u_k=c_1x_1\),这也就是初始条件\(u_0\)的第\(1\)个分量。

我们来证明第一个推论,取\(A=\begin{bmatrix}0.1&0.01&0.3\\0.2&0.99&0.3\\0.7&0&0.4\end{bmatrix}\),则\(A-I=\begin{bmatrix}-0.9&0.01&0.3\\0.2&-0.01&0.3\\0.7&0&-0.6\end{bmatrix}\)。观察\(A-I\)易知其列向量中元素之和均为\(0\),因为马尔科夫矩阵的性质就是各列向量元素之和为\(1\),现在我们从每一列中减去了\(1\),所以这是很自然的结果。而如果列向量中元素和为\(0\),则矩阵的任意行都可以用“零减去其他行之和”表示出来,即该矩阵的行向量线性相关。

用以前学过的子空间的知识描述,当\(n\)阶方阵各列向量元素之和皆为\(1\)时,则有\(\begin{bmatrix}1\\1\\\vdots\\1\end{bmatrix}\)在矩阵\(A-I\)左零空间中,即\((A-I)^T\)行向量线性相关。而\(A\)特征值\(1\)所对应的特征向量将在\(A-I\)的零空间中,因为\(Ax=x\rightarrow(A-I)x=0\)。

另外,特征值具有这样一个性质:矩阵与其转置的特征值相同。因为我们在行列式一讲了解了性质10,矩阵与其转置的行列式相同,那么如果\(\det(A-\lambda I)=0\),则有\(\det(A-\lambda I)^T=0\),根据矩阵转置的性质有\(\det(A^T-\lambda I^T)=0\),即\(\det(A^T-\lambda I)=0\)。这正是\(A^T\)特征值的计算式。

然后计算特征值\(\lambda_1=1\)所对应的特征向量,\((A-I)x_1=0\),得出\(x_1=\begin{bmatrix}0.6\\33\\0.7\end{bmatrix}\),特征向量中的元素皆为正。

接下来介绍马尔科夫矩阵的应用,我们用麻省和加州这两个州的人口迁移为例:

\(\begin{bmatrix}u_{cal}\\u_{mass}\end{bmatrix}_{k+1}\begin{bmatrix}0.9&0.2\\0.1&0.8\end{bmatrix}\begin{bmatrix}u_{cal}\\u_{mass}\end{bmatrix}_k\),元素非负,列和为一。这个式子表示每年有\(10%\)的人口从加州迁往麻省,同时有\(20%\)的人口从麻省迁往加州。注意使用马尔科夫矩阵的前提条件是随着时间的推移,矩阵始终不变。

设初始情况\(\begin{bmatrix}u_{cal}\\u_{mass}\end{bmatrix}_0=\begin{bmatrix}0\\1000\end{bmatrix}\),我们先来看第一次迁徙后人口的变化情况:\(\begin{bmatrix}u_{cal}\\u_{mass}\end{bmatrix}_1=\begin{bmatrix}0.9&0.2\\0.1&0.8\end{bmatrix}\begin{bmatrix}0\\1000\end{bmatrix}=\begin{bmatrix}200\\800\end{bmatrix}\),随着时间的推移,会有越来越多的麻省人迁往加州,而同时又会有部分加州人迁往麻省。

计算特征值:我们知道马尔科夫矩阵的一个特征值为\(\lambda_1=1\),则另一个特征值可以直接从迹算出\(\lambda_2=0.7\)。

计算特征向量:带入\(\lambda_1=1\)求\(A-I\)的零空间有\(\begin{bmatrix}-0.1&0.2\\0.1&-0.2\end{bmatrix}\),则\(x_1=\begin{bmatrix}2\\1\end{bmatrix}\),此时我们已经可以得出无穷步后稳态下的结果了。\(u_{\infty}=c_1\begin{bmatrix}2\\1\end{bmatrix}\)且人口总数始终为\(1000\),则\(c_1=\frac{1000}{3}\),稳态时\(\begin{bmatrix}u_{cal}\\u_{mass}\end{bmatrix}_{\infty}=\begin{bmatrix}\frac{2000}{3}\\\frac{1000}{3}\end{bmatrix}\)。注意到特征值为\(1\)的特征向量元素皆为正。

为了求每一步的结果,我们必须解出所有特征向量。带入\(\lambda_2=0.7\)求\(A-0.7I\)的零空间有\(\begin{bmatrix}0.2&0.2\\0.1&0.1\end{bmatrix}\),则\(x_2=\begin{bmatrix}-1\\1\end{bmatrix}\)。

通过\(u_0\)解出\(c_1, c_2\),\(u_k=c_11^k\begin{bmatrix}2\\1\end{bmatrix}+c_20.7^k\begin{bmatrix}-1\\1\end{bmatrix}\),带入\(k=0\)得\(u_0=\begin{bmatrix}0\\1000\end{bmatrix}=c_1\begin{bmatrix}2\\1\end{bmatrix}+c_2\begin{bmatrix}-1\\1\end{bmatrix}\),解出\(c_1=\frac{1000}{3}, c_2=\frac{2000}{3}\)。

另外,有时人们更喜欢用行向量,此时将要使用行向量乘以矩阵,其行向量各分量之和为\(1\)。

傅里叶级数

在介绍傅里叶级数(Fourier series)之前,先来回顾一下投影。

设\(q_1,q_2,\cdots q_n\)为一组标准正交基,则向量\(v\)在该标准正交基上的展开为\(v=x_1q_1+x_2q_2+\cdots+x_nq_n\),此时我们想要得到各系数\(x_i\)的值。比如求\(x_1\)的值,我们自然想要消掉除\(x_1q_1\)外的其他项,这时只需要等式两边同乘以\(q_1^T\),因为的\(q_i\)向量相互正交且长度为\(1\),则\(q_i^Tq_j=0, q_i^2=1\)所以原式变为\(q_1^Tv=x_1\)。

写为矩阵形式有\(\Bigg[q_1\ q_2\ \cdots\ q_n\Bigg]\begin{bmatrix}x_1\\x_2\\\vdots\\x_n\end{bmatrix}=v\),即\(Qx=v\)。所以有\(x=Q^{-1}v\),而在第十七讲我们了解到标准正交基有\(Q^T=Q^{-1}\),所以我们不需要计算逆矩阵可直接得出\(x=Q^Tv\)。此时对于\(x\)的每一个分量有\(x_i=q_i^Tv\)。

接下来介绍傅里叶级数。先写出傅里叶级数的展开式:

\[f(x)=a_0+a_1\cos x+b_1\sin x+a_2\cos 2x+b_2\sin 2x+\cdots \]

傅里叶发现,如同将向量\(v\)展开(投影)到向量空间的一组标准正交基中,在函数空间中,我们也可以做类似的展开。将函数\(f(x)\)投影在一系列相互正交的函数中。函数空间中的\(f(x)\)就是向量空间中的\(v\);函数空间中的\(1,\cos x,\sin x,\cos 2x,\sin 2x,\cdots\)就是向量空间中的\(q_1,q_2,\cdots,q_n\);不同的是,函数空间是无限维的而我们以前接触到的向量空间通常是有限维的。

再来介绍何为“函数正交”。对于向量正交我们通常使用两向量内积(点乘)为零判断。我们知道对于向量\(v,w\)的内积为\(v^Tw=v_1w_1+v_2w_2+\cdots+v_nw_n=0\),也就是向量的每个分量之积再求和。而对于函数\(f(x)\cdot g(x)\)内积,同样的,我们需要计算两个函数的每个值之积而后求和,由于函数取值是连续的,所以函数内积为:

\[f^Tg=\int f(x)g(x)\mathrm{d}x \]

在本例中,由于傅里叶级数使用正余弦函数,它们的周期都可以算作\(2\pi\),所以本例的函数点积可以写作\(f^Tg=\int_0^{2\pi}f(x)g(x)\mathrm{d}x\)。我来检验一个内积\(\int_0^{2\pi}\sin{x}\cos{x}\mathrm{d}x=\left.\frac{1}{2}\sin^2x\right|_0^{2\pi}=0\),其余的三角函数族正交性结果可以参考傅里叶级数的“希尔伯特空间的解读”一节。

最后我们来看\(\cos x\)项的系数是多少(\(a_0\)是\(f(x)\)的平均值)。同向量空间中的情形一样,我们在等式两边同时做\(\cos x\)的内积,原式变为\(\int_0^{2\pi}f(x)\cos x\mathrm{d}x=a_1\int_0^{2\pi}\cos^2x\mathrm{d}x\),因为正交性等式右边仅有\(\cos x\)项不为零。进一步化简得\(a_1\pi=\int_0^{2\pi}f(x)\cos x\mathrm{d}x\rightarrow a_1=\frac{1}{\pi}\int_0^{2\pi}f(x)\cos x\mathrm{d}x\)。

于是,我们把函数\(f(x)\)展开到了函数空间的一组标准正交基上。

第二十五讲:复习二

微分方程不在本讲的范围内。下面通过往年例题复习上面的知识。

  1. 求\(a=\begin{bmatrix}2\\1\\2\end{bmatrix}\)的投影矩阵\(P\):\(\Bigg(\)由\(a\bot(b-p)\rightarrow A^T(b-A\hat x)=0\)得到\(\hat x=\left(A^TA\right)^{-1}A^Tb\),求得\(p=A\hat x=A\left(A^TA\right)^{-1}A^Tb=Pb\)最终得到\(P\Bigg)\)\(\underline{P=A\left(A^TA\right)^{-1}A^T}\stackrel{a}=\frac{aa^T}{a^Ta}=\frac{1}{9}\begin{bmatrix}4&2&4\\2&1&2\\4&2&4\end{bmatrix}\)。

    求\(P\)矩阵的特征值:观察矩阵易知矩阵奇异,且为秩一矩阵,则其零空间为\(2\)维,所以由\(Px=0x\)得出矩阵的两个特征向量为\(\lambda_1=\lambda_2=0\);而从矩阵的迹得知\(trace(P)=1=\lambda_1+\lambda_2+\lambda_3=0+0+1\),则第三个特征向量为\(\lambda_3=1\)。

    求\(\lambda_3=1\)的特征向量:由\(Px=x\)我们知道经其意义为,\(x\)过矩阵\(P\)变换后不变,又有\(P\)是向量\(a\)的投影矩阵,所以任何向量经过\(P\)变换都会落在\(a\)的列空间中,则只有已经在\(a\)的列空间中的向量经过\(P\)的变换后保持不变,即其特征向量为\(x=a=\begin{bmatrix}2\\1\\2\end{bmatrix}\),也就是\(Pa=a\)。

    有差分方程\(u_{k+1}=Pu_k,\ u_0=\begin{bmatrix}9\\9\\0\end{bmatrix}\),求解\(u_k\):我们先不急于解出特征值、特征向量,因为矩阵很特殊(投影矩阵)。首先观察\(u_1=Pu_0\),式子相当于将\(u_0\)投影在了\(a\)的列空间中,计算得\(u_1=a\frac{a^Tu_0}{a^Ta}=3a=\begin{bmatrix}6\\3\\6\end{bmatrix}\)(这里的\(3\)相当于做投影时的系数\(\hat x\)),其意义为\(u_1\)在\(a\)上且距离\(u_0\)最近。再来看看\(u_2=Pu_1\),这个式子将\(u_1\)再次投影到\(a\)的列空间中,但是此时的\(u_1\)已经在该列空间中了,再次投影仍不变,所以有\(u_k=P^ku_0=Pu_0=\begin{bmatrix}6\\3\\6\end{bmatrix}\)。

    上面的解法利用了投影矩阵的特殊性质,如果在一般情况下,我们需要使用\(AS=S\Lambda\rightarrow A=S\Lambda S^{-1} \rightarrow u_{k+1}=Au_k=A^{k+1}u_0, u_0=Sc\rightarrow u_{k+1}=S\Lambda^{k+1}S^{-1}Sc=S\Lambda^{k+1}c\),最终得到公式\(A^ku_0=c_1\lambda_1^kx_1+c_2\lambda_2^kx_2+\cdots+c_n\lambda_n^kx_n\)。题中\(P\)的特殊性在于它的两个“零特征值”及一个“一特征值”使得式子变为\(A^ku_0=c_3x_3\),所以得到了上面结构特殊的解。

  2. 将点\((1,4),\ (2,5),\ (3,8)\)拟合到一条过零点的直线上:设直线为\(y=Dt\),写成矩阵形式为\(\begin{bmatrix}1\\2\\3\end{bmatrix}D=\begin{bmatrix}4\\5\\8\end{bmatrix}\),即\(AD=b\),很明显\(D\)不存在。利用公式\(A^TA\hat D=A^Tb\)得到\(14D=38,\ \hat D=\frac{38}{14}\),即最佳直线为\(y=\frac{38}{14}t\)。这个近似的意义是将\(b\)投影在了\(A\)的列空间中。

  3. 求\(a_1=\begin{bmatrix}1\\2\\3\end{bmatrix}\ a_2=\begin{bmatrix}1\\1\\1\end{bmatrix}\)的正交向量:找到平面\(A=\Bigg[a_1,a_2\Bigg]\)的正交基,使用Gram-Schmidt法,以\(a_1\)为基准,正交化\(a_2\),也就是将\(a_2\)中平行于\(a_1\)的分量去除,即\(a_2-xa_1=a_2-\frac{a_1^Ta_2}{a_1^Ta_1}a_1=\begin{bmatrix}1\\1\\1\end{bmatrix}-\frac{6}{14}\begin{bmatrix}1\\2\\3\end{bmatrix}\)。

  4. 有\(4\times 4\)矩阵\(A\),其特征值为\(\lambda_1,\lambda_2,\lambda_3,\lambda_4\),则矩阵可逆的条件是什么:矩阵可逆,则零空间中只有零向量,即\(Ax=0x\)没有非零解,则零不是矩阵的特征值。

    \(\det A^{-1}\)是什么:\(\det A^{-1}=\frac{1}{\det A}\),而\(\det A=\lambda_1\lambda_2\lambda_3\lambda_4\),所以有\(\det A^{-1}=\frac{1}{\lambda_1\lambda_2\lambda_3\lambda_4}\)。

    \(trace(A+I)\)的迹是什么:我们知道\(trace(A)=a_{11}+a_{22}+a_{33}+a_{44}=\lambda_1+\lambda_2+\lambda_3+\lambda_4\),所以有\(trace(A+I)=a_{11}+1+a_{22}+1+a_{33}+1+a_{44}+1=\lambda_1+\lambda_2+\lambda_3+\lambda_4+4\)。

  5. 有矩阵\(A_4=\begin{bmatrix}1&1&0&0\\1&1&1&0\\0&1&1&1\\0&0&1&1\end{bmatrix}\),求\(D_n=?D_{n-1}+?D_{n-2}\):求递归式的系数,使用代数余子式将矩阵安第一行展开得\(\det A_4=1\cdot\begin{vmatrix}1&1&0\\1&1&1\\0&1&1\end{vmatrix}-1\cdot\begin{vmatrix}1&1&0\\0&1&1\\0&1&1\end{vmatrix}=1\cdot\begin{vmatrix}1&1&0\\1&1&1\\0&1&1\end{vmatrix}-1\cdot\begin{vmatrix}1&1\\1&1\end{vmatrix}=\det A_3-\det A_2\)。则可以看出有规律\(D_n=D_{n-1}-D_{n-2}, D_1=1, D_2=0\)。

    使用我们在差分方程中的知识构建方程组\(\begin{cases}D_n&=D_{n-1}-D_{n-2}\\D_{n-1}&=D_{n-1}\end{cases}\),用矩阵表达有\(\begin{bmatrix}D_n\\D_{n-1}\end{bmatrix}=\begin{bmatrix}1&-1\\1&0\end{bmatrix}\begin{bmatrix}D_{n-1}\\D_{n-2}\end{bmatrix}\)。计算系数矩阵\(A_c\)的特征值,\(\begin{vmatrix}1-\lambda&1\\1&-\lambda\end{vmatrix}=\lambda^2-\lambda+1=0\),解得\(\lambda_1=\frac{1+\sqrt{3}i}{2},\lambda_2=\frac{1-\sqrt{3}i}{2}\),特征值为一对共轭复数。

    要判断递归式是否收敛,需要计算特征值的模,即实部平方与虚部平方之和\(\frac{1}{4}+\frac{3}{4}=1\)。它们是位于单位圆\(e^{i\theta}\)上的点,即\(\cos\theta+i\sin\theta\),从本例中可以计算出\(\theta=60^\circ\),也就是可以将特征值写作\(\lambda_1=e^{i\pi/3},\lambda_2=e^{-i\pi/3}\)。注意,从复平面单位圆上可以看出,这些特征值的六次方将等于一:\(e^{2\pi i}=e^{2\pi i}=1\)。继续深入观察这一特性对矩阵的影响,\(\lambda_1^6=\lambda^6=1\),则对系数矩阵有\(A_c^6=I\)。则系数矩阵\(A_c\)服从周期变化,既不发散也不收敛。

  6. 有这样一类矩阵\(A_4=\begin{bmatrix}0&1&0&0\\1&0&2&0\\0&2&0&3\\0&0&3&0\end{bmatrix}\),求投影到\(A_3\)列空间的投影矩阵:有\(A_3=\begin{bmatrix}0&1&0\\1&0&2\\0&2&0\end{bmatrix}\),按照通常的方法求\(P=A\left(A^TA\right)A^T\)即可,但是这样很麻烦。我们可以考察这个矩阵是否可逆,因为如果可逆的话,\(\mathbb{R}^4\)空间中的任何向量都会位于\(A_4\)的列空间,其投影不变,则投影矩阵为单位矩阵\(I\)。所以按行展开求行列式\(\det A_4=-1\cdot-1\cdot-3\cdot-3=9\),所以矩阵可逆,则\(P=I\)。

    求\(A_3\)的特征值及特征向量:\(\left|A_3-\lambda I\right|=\begin{vmatrix}-\lambda&1&0\\1&-\lambda&2\\0&2&-\lambda\end{vmatrix}=-\lambda^3+5\lambda=0\),解得\(\lambda_1=0,\lambda_2=\sqrt 5,\lambda_3=-\sqrt 5\)。

    我们可以猜测这一类矩阵的规律:奇数阶奇异,偶数阶可逆。

第二十六讲:对称矩阵及正定性

对称矩阵

前面我们学习了矩阵的特征值与特征向量,也了解了一些特殊的矩阵及其特征值、特征向量,特殊矩阵的特殊性应该会反映在其特征值、特征向量中。如马尔科夫矩阵,有一特征值为\(1\),本讲介绍(实)对称矩阵。

先提前介绍两个对称矩阵的特性:

  1. 特征值为实数;(对比第二十一讲介绍的旋转矩阵,其特征值为纯虚数。)
  2. 特征向量相互正交。(当特征值重复时,特征向量也可以从子空间中选出相互正交正交的向量。)

典型的状况是,特征值不重复,特征向量相互正交。

观察\((1)\)式,我们发现这个分解本身就代表着对称,\(\left(Q\varLambda Q^T\right)^T=\left(Q^T\right)^T\varLambda^TQ^T=Q\varLambda Q^T\)。\((1)\)式在数学上叫做谱定理(spectral theorem),谱就是指矩阵特征值的集合。(该名称来自光谱,指一些纯事物的集合,就像将特征值分解成为特征值与特征向量。)在力学上称之为主轴定理(principle axis theorem),从几何上看,它意味着如果给定某种材料,在合适的轴上来看,它就变成对角化的,方向就不会重复。

继续研究\(A=Q\varLambda Q^T=\Bigg[q_1\ q_2\ \cdots\ q_n\Bigg]\begin{bmatrix}\lambda_1& &\cdots& \\&\lambda_2&\cdots&\\\vdots&\vdots&\ddots&\vdots\\& &\cdots&\lambda_n\end{bmatrix}\begin{bmatrix}\quad q_1^T\quad\\\quad q_1^T\quad\\\quad \vdots \quad\\\quad q_1^T\quad\end{bmatrix}=\lambda_1q_1q_1^T+\lambda_2q_2q_2^T+\cdots+\lambda_nq_nq_n^T\),注意这个展开式中的\(qq^T\),\(q\)是单位列向量所以\(q^Tq=1\),结合我们在第十五讲所学的投影矩阵的知识有\(\frac{qq^T}{q^Tq}=qq^T\)是一个投影矩阵,很容易验证其性质,比如平方它会得到\(qq^Tqq^T=qq^T\)于是多次投影不变等。

每一个对称矩阵都可以分解为一系列相互正交的投影矩阵。

在知道对称矩阵的特征值皆为实数后,我们再来讨论这些实数的符号,因为特征值的正负号会影响微分方程的收敛情况(第二十三讲,需要实部为负的特征值保证收敛)。用消元法取得矩阵的主元,观察主元的符号,主元符号的正负数量与特征向量的正负数量相同

正定性

如果对称矩阵是“好矩阵”,则正定矩阵(positive definite)是其一个更好的子类。正定矩阵指特征值均为正数的矩阵(根据上面的性质有矩阵的主元均为正)。

举个例子,\(\begin{bmatrix}5&2\\2&3\end{bmatrix}\),由行列式消元知其主元为\(5,\frac{11}{5}\),按一般的方法求特征值有\(\begin{vmatrix}5-\lambda&2\\2&3-lambda\end{vmatrix}=\lambda^2-8\lambda+11=0, \lambda=4\pm\sqrt 5\)。

正定矩阵的另一个性质是,所有子行列式为正。对上面的例子有\(\begin{vmatrix}5\end{vmatrix}=5, \begin{vmatrix}5&2\\2&3\end{vmatrix}=11\)。

我们看到正定矩阵将早期学习的的消元主元、中期学习的的行列式、后期学习的特征值结合在了一起。

第二十七讲:复数矩阵和快速傅里叶变换

本讲主要介绍复数向量、复数矩阵的相关知识(包括如何做复数向量的点积运算、什么是复数对称矩阵等),以及傅里叶矩阵(最重要的复数矩阵)和快速傅里叶变换。

复数矩阵运算

先介绍复数向量,我们不妨换一个字母符号来表示:\(z=\begin{bmatrix}z_1\\z_2\\\vdots\\z_n\end{bmatrix}\),向量的每一个分量都是复数。此时\(z\)不再属于\(\mathbb{R}^n\)实向量空间,它现在处于\(\mathbb{C}^n\)复向量空间。

计算复向量的模

对比实向量,我们计算模只需要计算\(\left|v\right|=\sqrt{v^Tv}\)即可,而如果对复向量使用\(z^Tz\)则有\(z^Tz=\begin{bmatrix}z_1&z_2&\cdots&z_n\end{bmatrix}\begin{bmatrix}z_1\\z_2\\\vdots\\z_n\end{bmatrix}=z_1^2+z_2^2+\cdots+z_n^2\),这里\(z_i\)是复数,平方后虚部为负,求模时本应相加的运算变成了减法。(如向量\(\begin{bmatrix}1&i\end{bmatrix}\),右乘其转置后结果为\(0\),但此向量的长度显然不是零。)

根据上一讲我们知道,应使用\(\left|z\right|=\sqrt{\bar{z}^Tz}\),即\(\begin{bmatrix}\bar z_1&\bar z_2&\cdots&\bar z_n\end{bmatrix}\begin{bmatrix}z_1\\z_2\\\vdots\\z_n\end{bmatrix}\),即使用向量共轭的转置乘以原向量即可。(如向量\(\begin{bmatrix}1&i\end{bmatrix}\),右乘其共轭转置后结果为\(\begin{bmatrix}1&-i\end{bmatrix}\begin{bmatrix}1\\i\end{bmatrix}=2\)。)

我们把共轭转置乘以原向量记为\(z^Hz\),\(H\)读作埃尔米特(人名为Hermite,形容词为Hermitian)

计算向量的内积

有了复向量模的计算公式,同理可得,对于复向量,内积不再是实向量的\(y^Tx\)形式,复向量内积应为\(y^Hx\)。

对称性

对于实矩阵,\(A^T=A\)即可表达矩阵的对称性。而对于复矩阵,我们同样需要求一次共轭\(\bar{A}^T=A\)。举个例子\(\begin{bmatrix}2&3+i\\3-i&5\end{bmatrix}\)是一个复数情况下的对称矩阵。这叫做埃尔米特矩阵,有性质\(A^H=A\)。

正交性

在第十七讲中,我们这样定义标准正交向量:\(q_i^Tq_j=\begin{cases}0\quad i\neq j\\1\quad i=j\end{cases}\)。现在,对于复向量我们需要求共轭:\(\bar{q}_i^Tq_j=q_i^Hq_j=\begin{cases}0\quad i\neq j\\1\quad i=j\end{cases}\)。

第十七讲中的标准正交矩阵:\(Q=\Bigg[q_1\ q_2\ \cdots\ q_n\Bigg]\)有\(Q^TQ=I\)。现在对于复矩阵则有\(Q^HQ=I\)。

就像人们给共轭转置起了个“埃尔米特”这个名字一样,正交性(orthogonal)在复数情况下也有了新名字,酉(unitary),酉矩阵(unitary matrix)与正交矩阵类似,满足\(Q^HQ=I\)的性质。而前面提到的傅里叶矩阵就是一个酉矩阵。

傅里叶矩阵

\(n\)阶傅里叶矩阵\(F_n=\begin{bmatrix}1&1&1&\cdots&1\\1&w&w^2&\cdots&w^{n-1}\\1&w^2&w^4&\cdots&w^{2(n-1)}\\\vdots&\vdots&\vdots&\ddots&\vdots\\1&w^{n-1}&w^{2(n-1)}&\cdots&w^{(n-1)^2}\end{bmatrix}\),对于每一个元素有\((F_n)_{ij}=w^{ij}\quad i,j=0,1,2,\cdots,n-1\)。矩阵中的\(w\)是一个非常特殊的值,满足\(w^n=1\),其公式为\(w=e^{i2\pi/n}\)。易知\(w\)在复平面的单位圆上,\(w=\cos\frac{2\pi}{n}+i\sin\frac{2\pi}{n}\)。

在傅里叶矩阵中,当我们计算\(w\)的幂时,\(w\)在单位圆上的角度翻倍。比如在\(6\)阶情形下,\(w=e^{2\pi/6}\),即位于单位圆上\(60^\circ\)角处,其平方位于单位圆上\(120^\circ\)角处,而\(w^6\)位于\(1\)处。从开方的角度看,它们是\(1\)的\(6\)个六次方根,而一次的\(w\)称为原根。

快速傅里叶变换(Fast Fourier transform/FFT)

对于傅里叶矩阵,\(F_6,\ F_3\)、\(F_8,\ F_4\)、\(F_{64},\ F_{32}\)之间有着特殊的关系。

举例,有傅里叶矩阵\(F_64\),一般情况下,用一个列向量右乘\(F_{64}\)需要约\(64^2\)次计算,显然这个计算量是比较大的。我们想要减少计算量,于是想要分解\(F_{64}\),联系到\(F_{32}\),有\(\Bigg[F_{64}\Bigg]=\begin{bmatrix}I&D\\I&-D\end{bmatrix}\begin{bmatrix}F_{32}&0\\0&F_{32}\end{bmatrix}\begin{bmatrix}1&&\cdots&&&0&&\cdots&&\\0&&\cdots&&&1&&\cdots&&\\&1&\cdots&&&&0&\cdots&&\\&0&\cdots&&&&1&\cdots&&\\&&&\ddots&&&&&\ddots&&\\&&&\ddots&&&&&\ddots&&\\&&&\cdots&1&&&&\cdots&0\\&&&\cdots&0&&&&\cdots&1\end{bmatrix}\)。

我们分开来看等式右侧的这三个矩阵:

所以我们把\(64^2\)复杂度的计算化简为\(2\times 32^2+32\)复杂度的计算,我们可以进一步化简\(F_{32}\)得到与\(F_{16}\)有关的式子\(\begin{bmatrix}I_{32}&D_{32}\\I_{32}&-D_{32}\end{bmatrix}\begin{bmatrix}I_{16}&D_{16}&&\\I_{16}&-D_{16}&&\\&&I_{16}&D_{16}\\&&I_{16}&-D_{16}\end{bmatrix}\begin{bmatrix}F_{16}&&&\\&F_{16}&&\\&&F_{16}&\\&&&F_{16}\end{bmatrix}\begin{bmatrix}P_{16}&\\&P_{16}\end{bmatrix}\Bigg[\ P_{32}\ \Bigg]\)。而\(32^2\)的计算量进一步分解为\(2\times 16^2+16\)的计算量,如此递归下去我们最终得到含有一阶傅里叶矩阵的式子。

来看化简后计算量,\(2\left(2\left(2\left(2\left(2\left(2\left(1\right)^2+1\right)+2\right)+4\right)+8\right)+16\right)+32\),约为\(6\times 32=\log_264\times \frac{64}{2}\),算法复杂度为\(\frac{n}{2}\log_2n\)。

于是原来需要\(n^2\)的运算现在只需要\(\frac{n}{2}\log_2n\)就可以实现了。不妨看看\(n=10\)的情况,不使用FFT时需要\(n^2=1024\times 1024\)次运算,使用FFT时只需要\(\frac{n}{2}\log_2n=5\times 1024\)次运算,运算量大约是原来的\(\frac{1}{200}\)。

下一讲将继续介绍特征值、特征向量及正定矩阵。

第二十八讲:正定矩阵和最小值

本讲我们会了解如何完整的测试一个矩阵是否正定,测试\(x^TAx\)是否具有最小值,最后了解正定的几何意义——椭圆(ellipse)和正定性有关,双曲线(hyperbola)与正定无关。另外,本讲涉及的矩阵均为实对称矩阵。

正定性的判断

我们仍然从二阶说起,有矩阵\(A=\begin{bmatrix}a&b\\b&d\end{bmatrix}\),判断其正定性有以下方法:

  1. 矩阵的所有特征值大于零则矩阵正定:\(\lambda_1>0,\ \lambda_2>0\);
  2. 矩阵的所有顺序主子阵(leading principal submatrix)的行列式(即顺序主子式,leading principal minor)大于零则矩阵正定:\(a>0,\ ac-b^2>0\);
  3. 矩阵消元后主元均大于零:\(a>0,\ \frac{ac-b^2}{a}>0\);
  4. \(x^TAx>0\);

大多数情况下使用4来定义正定性,而用前三条来验证正定性。

来计算一个例子:\(A=\begin{bmatrix}2&6\\6&?\end{bmatrix}\),在\(?\)处填入多少才能使矩阵正定?

接下来计算一个三阶矩阵,\(A=\begin{bmatrix}2&-1&0\\-1&2&-1\\0&-1&2\end{bmatrix}\),它是正定的吗?函数\(x^TAx\)是多少?函数在原点去最小值吗?图像是什么样的?

现在我们将矩阵\(A\)分解为\(A=Q\Lambda Q^T\),可以发现上面说到的各种元素都可以表示在这个分解的矩阵中,我们称之为主轴定理(principal axis theorem),即特征向量说明主轴的方向、特征值说明主轴的长度。

\(A=Q\Lambda Q^T\)是特征值相关章节中最重要的公式。

第二十九讲:相似矩阵和若尔当形

在本讲的开始,先接着上一讲来继续说一说正定矩阵。

接下来进入本讲的正题。

相似矩阵

先列出定义:矩阵\(A,\ B\)对于某矩阵\(M\)满足\(B=M^{-1}AM\)时,成\(A,\ B\)互为相似矩阵。

对于在对角化一讲(第二十二讲)中学过的式子\(S^{-1}AS=\Lambda\),则有\(A\)相似于\(\Lambda\)。

所以,相似矩阵有相同的特征值。

现在我们来证明这个性质,有\(Ax=\lambda x,\ B=M^{-1}AM\),第一个式子化为\(AMM^{-1}x=\lambda x\),接着两边同时左乘\(M^{-1}\)得\(M^{-1}AMM^{-1}x=\lambda M^{-1}x\),进行适当的分组得\(\left(M^{-1}AM\right)M^{-1}x=\lambda M^{-1}x\)即\(BM^{-1}x=\lambda M^{-1}x\)。

\(BM^{-1}=\lambda M^{-1}x\)可以解读成矩阵\(B\)与向量\(M^{-1}x\)之积等于\(\lambda\)与向量\(M^{-1}x\)之积,也就是\(B\)的仍为\(\lambda\),而特征向量变为\(M^{-1}x\)。

以上就是我们得到的一族特征值为\(3,\ 1\)的矩阵,它们具有相同的特征值。接下来看特征值重复时的情形。

若尔当形在过去是线性代数的核心知识,但现在不是了(现在是下一讲的奇异值分解),因为它并不容易计算。

若尔当形

再来看一个更加“糟糕”的矩阵:

若尔当块的定义型为\(J_i=\begin{bmatrix}\lambda_i&1&&\cdots&\\&\lambda_i&1&\cdots&\\&&\lambda_i&\cdots&\\\vdots&\vdots&\vdots&\ddots&\\&&&&\lambda_i\end{bmatrix}\),它的对角线上只为同一个数,仅有一个特征向量。

所有有,每一个矩阵\(A\)都相似于一个若尔当矩阵,型为\(J=\left[\begin{array}{c|c|c|c}J_1&&&\\\hline&J_2&&\\\hline&&\ddots&\\\hline&&&J_d\end{array}\right]\)。注意,对角线上方还有\(1\)。若尔当块的个数即为矩阵特征值的个数。

在矩阵为“好矩阵”的情况下,\(n\)阶矩阵将有\(n\)个不同的特征值,那么它可以对角化,所以它的若尔当矩阵就是\(\Lambda\),共\(n\)个特征向量,有\(n\)个若尔当块。

第三十讲:奇异值分解

本讲我们介绍将一个矩阵写为\(A=U\varSigma V^T\),分解的因子分别为正交矩阵、对角矩阵、正交矩阵,与前面几讲的分解不同的是,这两个正交矩阵通常是不同的,而且这个式子可以对任意矩阵使用,不仅限于方阵、可对角化的方阵等。

我们现在要做的是,在\(A\)的列空间中找到一组特殊的正交基\(v_1,v_2,\cdots,v_r\),这组基在\(A\)的作用下可以转换为\(A\)的行空间中的一组正交基\(u_1,u_2,\cdots,u_r\)。

用矩阵语言描述为\(A\Bigg[v_1\ v_2\ \cdots\ v_r\Bigg]=\Bigg[\sigma_1u_1\ \sigma_2u_2\ \cdots\ \sigma_ru_r\Bigg]=\Bigg[u_1\ u_2\ \cdots\ u_r\Bigg]\begin{bmatrix}\sigma_1&&&\\&\sigma_2&&\\&&\ddots&\\&&&\sigma_n\end{bmatrix}\),即\(Av_1=\sigma_1u_1,\ Av_2=\sigma_2u_2,\cdots,Av_r=\sigma_ru_r\),这些\(\sigma\)是缩放因子,表示在转换过程中有拉伸或压缩。而\(A\)的左零空间和零空间将体现在\(\sigma\)的零值中。

另外,如果算上左零、零空间,我们同样可以对左零、零空间取标准正交基,然后写为\(A\Bigg[v_1\ v_2\ \cdots\ v_r\ v_{r+1}\ \cdots\ v_m\Bigg]=\Bigg[u_1\ u_2\ \cdots\ u_r\ u_{r+1}\ \cdots \ u_n\Bigg]\left[\begin{array}{c c c|c}\sigma_1&&&\\&\ddots&&\\&&\sigma_r&\\\hline&&&\begin{bmatrix}0\end{bmatrix}\end{array}\right]\),此时\(U\)是\(m\times m\)正交矩阵,\(\varSigma\)是\(m\times n\)对角矩阵,\(V^T\)是\(n\times n\)正交矩阵。

最终可以写为\(AV=U\varSigma\),可以看出这十分类似对角化的公式,矩阵\(A\)被转化为对角矩阵\(\varSigma\),我们也注意到\(U,\ V\)是两组不同的正交基。(在正定的情况下,\(U,\ V\)都变成了\(Q\)。)。进一步可以写作\(A=U\varSigma V^{-1}\),因为\(V\)是标准正交矩阵所以可以写为\(A=U\varSigma V^T\)

计算一个例子,\(A=\begin{bmatrix}4&4\\-3&3\end{bmatrix}\),我们需要找到:

在\(A=U\varSigma V^T\)中有两个标准正交矩阵需要求解,我们希望一次只解一个,如何先将\(U\)消去来求\(V\)?

这个技巧会经常出现在长方形矩阵中:求\(A^TA\),这是一个对称正定矩阵(至少是半正定矩阵),于是有\(A^TA=V\varSigma^TU^TU\varSigma V^T\),由于\(U\)是标准正交矩阵,所以\(U^TU=I\),而\(\varSigma^T\varSigma\)是对角线元素为\(\sigma^2\)的对角矩阵。

现在有\(A^TA=V\begin{bmatrix}\sigma_1&&&\\&\sigma_2&&\\&&\ddots&\\&&&\sigma_n\end{bmatrix}V^T\),这个式子中\(V\)即是\(A^TA\)的特征向量矩阵而\(\varSigma^2\)是其特征值矩阵。

同理,我们只想求\(U\)时,用\(AA^T\)消掉\(V\)即可。

我们来计算\(A^TA=\begin{bmatrix}4&-3\\4&3\end{bmatrix}\begin{bmatrix}4&4\\-3&3\end{bmatrix}=\begin{bmatrix}25&7\\7&25\end{bmatrix}\),对于简单的矩阵可以直接观察得到特征向量\(A^TA\begin{bmatrix}1\\1\end{bmatrix}=32\begin{bmatrix}1\\1\end{bmatrix},\ A^TA\begin{bmatrix}1\\-1\end{bmatrix}=18\begin{bmatrix}1\\-1\end{bmatrix}\),化为单位向量有\(\sigma_1=32,\ v_1=\begin{bmatrix}\frac{1}{\sqrt{2}}\\\frac{1}{\sqrt{2}}\end{bmatrix},\ \sigma_2=18,\ v_2=\begin{bmatrix}\frac{1}{\sqrt{2}}\\-\frac{1}{\sqrt{2}}\end{bmatrix}\)。

到目前为止,我们得到\(\begin{bmatrix}4&4\\-3&3\end{bmatrix}=\begin{bmatrix}u_?&u_?\\u_?&u_?\end{bmatrix}\begin{bmatrix}\sqrt{32}&0\\0&\sqrt{18}\end{bmatrix}\begin{bmatrix}\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}\\\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}\end{bmatrix}\),接下来继续求解\(U\)。

\(AA^T=U\varSigma V^TV\varSigma^TU^T=U\varSigma^2U^T\),求出\(AA^T\)的特征向量即可得到\(U\),\(\begin{bmatrix}4&4\\-3&3\end{bmatrix}\begin{bmatrix}4&-3\\4&3\end{bmatrix}=\begin{bmatrix}32&0\\0&18\end{bmatrix}\),观察得\(AA^T\begin{bmatrix}1\\0\end{bmatrix}=32\begin{bmatrix}1\\0\end{bmatrix},\ AA^T\begin{bmatrix}0\\1\end{bmatrix}=18\begin{bmatrix}0\\1\end{bmatrix}\)。但是我们不能直接使用这一组特征向量,因为式子\(AV=U\varSigma\)明确告诉我们,一旦\(V\)确定下来,\(U\)也必须取能够满足该式的向量,所以此处\(Av_2=\begin{bmatrix}0\\-\sqrt{18}\end{bmatrix}=u_2\sigma_2=\begin{bmatrix}0\\-1\end{bmatrix}\sqrt{18}\),则\(u_1=\begin{bmatrix}1\\0\end{bmatrix},\ u_2=\begin{bmatrix}0\\-1\end{bmatrix}\)。(这个问题在本讲的官方笔记中有详细说明。)

最终,我们得到\(\begin{bmatrix}4&4\\-3&3\end{bmatrix}=\begin{bmatrix}1&0\\0&-1\end{bmatrix}\begin{bmatrix}\sqrt{32}&0\\0&\sqrt{18}\end{bmatrix}\begin{bmatrix}\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}\\\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}\end{bmatrix}\)。

再做一个例子,\(A=\begin{bmatrix}4&3\\8&6\end{bmatrix}\),这是个秩一矩阵,有零空间。\(A\)的行空间为\(\begin{bmatrix}4\\3\end{bmatrix}\)的倍数,\(A\)的列空间为\(\begin{bmatrix}4\\8\end{bmatrix}\)的倍数。

最终得到\(\begin{bmatrix}4&3\\8&6\end{bmatrix}=\begin{bmatrix}1&\underline {2}\\2&\underline{-1}\end{bmatrix}\begin{bmatrix}\sqrt{125}&0\\0&\underline{0}\end{bmatrix}\begin{bmatrix}0.8&0.6\\\underline{0.6}&\underline{-0.8}\end{bmatrix}\),其中下划线部分都是与零空间相关的部分。

通过将矩阵写为\(Av_i=\sigma_iu_i\)形式,将矩阵对角化,向量\(u,\ v\)之间没有耦合,\(A\)乘以每个\(v\)都能得到一个相应的\(u\)。

第三十一讲:线性变换及对应矩阵

如何判断一个操作是不是线性变换?线性变换需满足以下两个要求:

\[T(v+w)=T(v)+T(w)\\ T(cv)=cT(v) \]

即变换\(T\)需要同时满足加法和数乘不变的性质。将两个性质合成一个式子为:\(T(cv+dw)=cT(v)+dT(w)\)

例1,二维空间中的投影操作,\(T: \mathbb{R}^2\to\mathbb{R}^2\),它可以将某向量投影在一条特定直线上。检查一下投影操作,如果我们将向量长度翻倍,则其投影也翻倍;两向量相加后做投影与两向量做投影再相加结果一致。所以投影操作是线性变换。

“坏”例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()

png

plt.close(fig)

比如,上图中向量长度翻倍,再做平移,明显与向量平移后再翻倍的结果不一致。

有时我们也可以用一个简单的特例判断线性变换,检查\(T(0)\stackrel{?}{=}0\)。零向量平移后结果并不为零。

所以平面平移操作并不是线性变换。

“坏”例2,求模运算,\(T(v)=\|v\|,\ T:\mathbb{R}^3\to\mathbb{R}^1\),这显然不是线性变换,比如如果我们将向量翻倍则其模翻倍,但如果我将向量翻倍取负,则其模依然翻倍。所以\(T(-v)\neq -T(v)\)

例2,旋转\(45^\circ\)操作,\(T:\mathbb{R}^2\to\mathbb{R}^2\),也就是将平面内一个向量映射为平面内另一个向量。检查可知,如果向量翻倍,则旋转后同样翻倍;两个向量先旋转后相加,与这两个向量先相加后旋转得到的结果一样。

所以从上面的例子我们知道,投影与旋转都是线性变换。

例3,矩阵乘以向量,\(T(v)=Av\),这也是一个(一系列)线性变换,不同的矩阵代表不同的线性变换。根据矩阵的运算法则有\(A(v+w)=A(v)+A(w),\ A(cv)=cAv\)。比如取\(A=\begin{bmatrix}1&0\\0&-1\end{bmatrix}\),作用于平面上的向量\(v\),会导致\(v\)的\(x\)分量不变,而\(y\)分量取反,也就是图像沿\(x\)轴翻转。

线性变换的核心,就是该变换使用的相应的矩阵。

比如我们需要做一个线性变换,将一个三维向量降至二维,\(T:\mathbb{R}^3\to\mathbb{R}^2\),则在\(T(v)=Av\)中,\(v\in\mathbb{R}^3,\ T(v)\in\mathbb{R}^2\),所以\(A\)应当是一个\(2\times 3\)矩阵。

如果我们希望知道线性变换\(T\)对整个输入空间\(\mathbb{R}^n\)的影响,我们可以找到空间的一组基\(v_1,\ v_2,\ \cdots,\ v_n\),检查\(T\)对每一个基的影响\(T(v_1),\ T(v_2),\ \cdots,\ T(v_n)\),由于输入空间中的任意向量都满足:

\[v=c_1v_1+c_2v_2+\cdots+c_nv_n\tag{1} \]

所以我们可以根据\(T(v)\)推出线性变换\(T\)对空间内任意向量的影响,得到:

\[T(v)=c_1T(v_1)+c_2T(v_2)+\cdots+c_nT(v_n)\tag{2} \]

现在我们需要考虑,如何把一个与坐标无关的线性变换变成一个与坐标有关的矩阵呢?

在\(1\)式中,\(c_1,c_2,\cdots,c_n\)就是向量\(v\)在基\(v_1,v_2,\cdots,v_n\)上的坐标,比如分解向量\(v=\begin{bmatrix}3\\2\\4\end{bmatrix}=3\begin{bmatrix}1\\0\\0\end{bmatrix}+2\begin{bmatrix}0\\1\\0\end{bmatrix}+4\begin{bmatrix}0\\0\\1\end{bmatrix}\),式子将向量\(v\)分解在一组标准正交基\(\begin{bmatrix}1\\0\\0\end{bmatrix},\begin{bmatrix}0\\1\\0\end{bmatrix},\begin{bmatrix}0\\0\\1\end{bmatrix}\)上。当然,我们也可以选用矩阵的特征向量作为基向量,基的选择是多种多样的。

我们打算构造一个矩阵\(A\)用以表示线性变换\(T:\mathbb{R}^n\to\mathbb{R}^m\)。我们需要两组基,一组用以表示输入向量,一组用以表示输出向量。令\(v_1,v_2,\cdots,v_n\)为输入向量的基,这些向量来自\(\mathbb{R}^n\);\(w_1,w_2,\cdots,w_m\)作为输出向量的基,这些向量来自\(\mathbb{R}^m\)。

我们用二维空间的投影矩阵作为例子:

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()

png

plt.close(fig)

从图中可以看到,设输入向量的基为\(v_1,v_2\),\(v_1\)就在投影上,而\(v_2\)垂直于投影方向,输出向量的基为\(w_1,w_2\),而\(v_1=w_1,v_2=w_2\)。那么如果输入向量为\(v=c_1v_1+c_2v_2\),则输出向量为\(T(v)=c_1v_1\),也就是线性变换去掉了法线方向的分量,输入坐标为\((c_1,c_2)\),输出坐标变为\((c_1,0)\)。

找出这个矩阵并不困难,\(Av=w\),则有\(\begin{bmatrix}1&0\\0&0\end{bmatrix}\begin{bmatrix}c_1\\c_2\end{bmatrix}=\begin{bmatrix}c_1\\0\end{bmatrix}\)。

本例中我们选取的基极为特殊,一个沿投影方向,另一个沿投影法线方向,其实这两个向量都是投影矩阵的特征向量,所以我们得到的线性变换矩阵是一个对角矩阵,这是一组很好的基。

所以,如果我们选取投影矩阵的特征向量作为基,则得到的线性变换矩阵将是一个包含投影矩阵特征值的对角矩阵。

继续这个例子,我们不再选取特征向量作为基,而使用标准基\(v_1=\begin{bmatrix}1\\0\end{bmatrix},v_2=\begin{bmatrix}0\\1\end{bmatrix}\),我们继续使用相同的基作为输出空间的基,即\(v_1=w_1,v_2=w_2\)。此时投影矩阵为\(P=\frac{aa^T}{a^Ta}=\begin{bmatrix}\frac{1}{2}&\frac{1}{2}\\\frac{1}{2}&\frac{1}{2}\end{bmatrix}\),这个矩阵明显没有上一个矩阵“好”,不过这个矩阵也是一个不错的对称矩阵。

总结通用的计算线性变换矩阵\(A\)的方法:

最后我们介绍一种不一样的线性变换,\(T=\frac{\mathrm{d}}{\mathrm{d}x}\):

最后,矩阵的逆相当于对应线性变换的逆运算,矩阵的乘积相当于线性变换的乘积,实际上矩阵乘法也源于线性变换。

第三十二讲:基变换和图像压缩

图像压缩

本讲我们介绍一种图片有损压缩的一种方法:JPEG。

假设我们有一张图片,长宽皆为\(512\)个像素,我们用\(x_i\)来表示第\(i\)个像素,如果是灰度照片,通常\(x_i\)可以在\([0,255]\)上取值,也就是8 bits。对于这承载这张图片信息的向量\(x\)来说,有\(x\in\mathbb{R}^n, n=512^2\)。而如果是彩色照片,通常需要三个量来表示一个像素,则向量长度也会变为现在的三倍。

如此大的数据不经过压缩很难广泛传播。教学录像采用的压缩方法就是JPEG(Joint Photographic Expert Group,联合图像专家组),该方法采用的就是基变换的方式压缩图像。比如说一块干净的黑白,其附近的像素值应该非常接近,此时如果一个像素一个像素的描述黑白灰度值就太浪费空间了,所以标准基在这种情况下并不能很好的利用图片的特性。

我们知道,标准基是 \(\begin{bmatrix}1\\0\\\vdots\\0\end{bmatrix}\begin{bmatrix}0\\1\\\vdots\\0\end{bmatrix}\cdots\begin{bmatrix}0\\0\\\vdots\\1\end{bmatrix}\),我们想寻找一个更好的基。

我们试试使用别的基描述图片,比如:

傅里叶基

现在我们来介绍傅里叶基,以\(8\times 8\)傅里叶基为例(这表示我们每次只处理\(8\times 8\)像素的一小块图像):

\(F_n=\begin{bmatrix}1&1&1&\cdots&1\\1&w&w^2&\cdots&w^{n-1}\\1&w^2&w^4&\cdots&w^{2(n-1)}\\\vdots&\vdots&\vdots&\ddots&\vdots\\1&w^{n-1}&w^{2(n-1)}&\cdots&w^{(n-1)^2}\end{bmatrix},\ w=e^{i2\pi/n},\ n=8\),我们不需要深入\(8\)阶傅里叶基的细节,先看看使用傅里叶基的思路是怎样的。

每次处理\(8\times 8\)的一小块时,会遇到\(64\)个像素,也就是\(64\)个基向量,\(64\)个系数,在\(64\)维空间中利用傅里叶向量做基变换:

我们再来提一下视频压缩:视频是一系列连续图像,且相近的帧非常接近,而我们的压缩算法就需要利用这个相近性质。在实际生活中,从时间与空间的角度讲,事物不会瞬间改变。

小波基

接下来介绍另一组基,它是傅里叶基的竞争对手,名为小波(wavelets),同样以\(8\times 8\)为例:
\(\begin{bmatrix}1\\1\\1\\1\\1\\1\\1\\1\end{bmatrix} \begin{bmatrix}1\\1\\1\\1\\-1\\-1\\-1\\-1\end{bmatrix} \begin{bmatrix}1\\1\\-1\\-1\\0\\0\\0\\0\end{bmatrix} \begin{bmatrix}0\\0\\0\\0\\1\\1\\-1\\-1\end{bmatrix} \begin{bmatrix}1\\-1\\0\\0\\0\\0\\0\\0\end{bmatrix} \begin{bmatrix}0\\0\\1\\-1\\0\\0\\0\\0\end{bmatrix} \begin{bmatrix}0\\0\\0\\0\\1\\-1\\0\\0\end{bmatrix} \begin{bmatrix}0\\0\\0\\0\\0\\0\\1\\-1\end{bmatrix}\)。

可以看出傅里叶基中频率最高的向量为小波后四个基向量之和。

在标准基下的一组(按八个一组计算,\(P\in\mathbb{R}^8\))像素\(P=\begin{bmatrix}p_1\\p_2\\\vdots\\p_8\end{bmatrix}=c_1w_1+c_2w_2+\cdots+c_nw_n=\Bigg[w_1\ w_2\ \cdots\ w_n\Bigg]\begin{bmatrix}c_1\\c_2\\\vdots\\c_n\end{bmatrix}\),即\(P=WC\),我们需要计算像素向量在另一组基下系数,所以有\(C=W^{-1}P\)。

此时我们发现,如果选取“好的基”会使得逆矩阵的求解过程变简单,所谓“好的基”:

要想继续了解小波基,可以参考一篇非常精彩的文章能不能通俗的讲解下傅立叶分析和小波分析之间的关系?——“咚懂咚懂咚“的答案

基变换

前面介绍小波基的时候我们就已经做了一次基变换。

将目标基的向量按列组成矩阵\(W\),则基变换就是\(\Bigg[x\Bigg]\xrightarrow{x=Wc}\Bigg[c\Bigg]\)。

看一个例子,有线性变换\(T:\mathbb{R}^8\to\mathbb{R}^8\),在第一组基\(v_1,v_2,\cdots,v_8\)上计算得到矩阵\(A\),在第二组基\(w_1,w_2,\cdots,w_n\)上计算得到矩阵\(B\)。先说结论,矩阵\(A,B\)是相似的,也就是有\(B=M^{-1}AM\),而\(M\)就是基变换矩阵。

进行基变换时会发生两件事:

  1. 每个向量都会有一组新的坐标,而\(x=Wc\)就是新旧坐标的关系;

  2. 每个线性变换都会有一个新的矩阵,而\(B=M^{-1}AM\)就是新旧矩阵的关系。

    再来看什么是\(A\)矩阵?

    对于第一组基\(v_1,v_2,\cdots,v_8\),要完全了解线性变换\(T\),只需要知道\(T\)作用在基的每一个向量上会产生什么结果即可。因为在这个基下的每一个向量都可以写成\(x=c_1v_1+c_2v_2+\cdots+c_8v_8\)的形式,所以\(T(x)=c_1T(v_1)+c_2T(v_2)+\cdots+c_8T(v_8)\)。

    而且\(T(v_1)=a_{11}v_1+a_{21}v_2+\cdots+a_{81}v_8,\ T(v_2)=a_{12}v_1+a_{22}v_2+\cdots+a_{82}v_8,\ \cdots\),则矩阵\(\begin{bmatrix}A\end{bmatrix}=\left[\begin{array}{c|c|c|c}a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{m1}&a_{m2}&\cdots&a_{mn}\\\end{array}\right]\)

    这些都是上一讲结尾所涉及的知识。

最后我们以一个更加特殊的基收场,设\(v_1,v_2,\cdots,v_n\)是一组特征向量,也就是\(T(v_i)=\lambda_1v_i\),那么问题就是矩阵\(A\)是什么?

继续使用线性变换中学到的,输入的第一个向量\(v_1\)经由\(T\)加工后得到\(\lambda_1v_1\),第二个向量\(v_2\xrightarrow{T}\lambda_2v_2\),继续做下去,最终有\(v_n=v_n\xrightarrow{T}\lambda_nv_n\)。除了\(\lambda_iv_i\)外的其他基向量都变为\(0\),那么矩阵\(A=\begin{bmatrix}\lambda_1&&&\\&\lambda_2&&\\&&\ddots&\\&&&\lambda_n\end{bmatrix}\)。

这是一个非常完美的基,我们在图像处理中最想要的就是这种基,但是找出像素矩阵的特征向量代价太大,所以我们找了一些代价小同时效果也不错的基,比如小波基、傅里叶基等等。

第三十三讲:单元检测3复习

在上一次复习中,我们已经涉及了求特征值与特征向量(通过解方程\(\det(A-\lambda I)=0\)得出\(\lambda\),再将\(\lambda\)带入\(A-\lambda I\)求其零空间得到\(x\))。

接下的章节来我们学习了:

现在,我们继续通过例题复习这些知识点。

  1. 解方程\(\frac{\mathrm{d}u}{\mathrm{d}t}=Au=\begin{bmatrix}0&-1&0\\1&0&-1\\0&1&0\end{bmatrix}u\)

    首先通过\(A\)的特征值/向量求通解\(u(t)=c_1e^{\lambda_1t}x_1+c_2e^{\lambda_2t}x_2+c_3e^{\lambda_3t}x_3\),很明显矩阵是奇异的,所以有\(\lambda_1=0\);

    继续观察矩阵会发现\(A^T=-A\),这是一个反对称矩阵(anti-symmetric)或斜对陈矩阵(skew-symmetric),这与我们在第二十一讲介绍过的旋转矩阵类似,它的特征值应该为纯虚数(特征值在虚轴上),所以我们猜测其特征值应为\(0\cdot i,\ b\cdot i,\ -b\cdot i\)。通过解\(\det(A-\lambda I)=0\)验证一下:\(\begin{bmatrix}-\lambda&-1&0\\1&-\lambda&-1\\0&1&\lambda\end{bmatrix}=\lambda^3+2\lambda=0, \lambda_2=\sqrt 2i, \lambda_3=-\sqrt 2i\)。

    此时\(u(t)=c_1+c_2e^{\sqrt 2it}x_2+c_3e^{-\sqrt 2it}x_3\),\(e^{\sqrt 2it}\)始终在复平面单位圆上,所以\(u(t)\)及不发散也不收敛,它只是具有周期性。当\(t=0\)时有\(u(0)=c_1+c_2+c_3\),如果使\(e^{\sqrt 2iT}=1\)即\(\sqrt 2iT=2\pi i\)则也能得到\(u(T)=c_1+c_2+c_3\),周期\(T=\pi\sqrt 2\)。

    另外,反对称矩阵同对称矩阵一样,具有正交的特征向量。当矩阵满足什么条件时,其特征向量相互正交?答案是必须满足\(AA^T=A^TA\)。所以对称矩阵\(A=A^T\)满足此条件,同时反对称矩阵\(A=-A^T\)也满足此条件,而正交矩阵\(Q^{-1}=Q^T\)同样满足此条件,这三种矩阵的特征向量都是相互正交的。

    上面的解法并没有求特征向量,进而通过\(u(t)=e^{At}u(0)\)得到通解,现在我们就来使用指数矩阵来接方程。如果矩阵可以对角化(在本例中显然可以),则\(A=S\Lambda S^{-1}, e^{At}=Se^{\Lambda t}S^{-1}=S\begin{bmatrix}e^{\lambda_1t}&&&\\&e^{\lambda_1t}&&\\&&\ddots&\\&&&e^{\lambda_1t}\end{bmatrix}S^{-1}\),这个公式在能够快速计算\(S,\lambda\)时很方便求解。

  2. 已知矩阵的特征值\(\lambda_1=0,\lambda_2=c,\lambda_3=2\),特征向量\(x_1=\begin{bmatrix}1\\1\\1\end{bmatrix},x_2=\begin{bmatrix}1&-1&0\end{bmatrix},x_3=\begin{bmatrix}1\\1\\-2\end{bmatrix}\):

    \(c\)如何取值才能保证矩阵可以对角化?其实可对角化只需要有足够的特征向量即可,而现在特征向量已经足够,所以\(c\)可以取任意值。

    \(c\)如何取值才能保证矩阵对称?我们知道,对称矩阵的特征值均为实数,且注意到给出的特征向量是正交的,有了实特征值及正交特征向量,我们就可以得到对称矩阵。

    \(c\)如何取值才能使得矩阵正定?已经有一个零特征值了,所以矩阵不可能是正定的,但可以是半正定的,如果\(c\)去非负实数。

    \(c\)如何取值才能使得矩阵是一个马尔科夫矩阵?在第二十四讲我们知道马尔科夫矩阵的性质:必有特征值等于\(1\),其余特征值均小于\(1\),所以\(A\)不可能是马尔科夫矩阵。

    \(c\)取何值才能使得\(P=\frac{A}{2}\)是一个投影矩阵?我们知道投影矩阵的一个重要性质是\(P^2=P\),所以有对其特征值有\(\lambda^2=\lambda\),则\(c=0,2\)。

    题设中的正交特征向量意义重大,如果没有正交这个条件,则矩阵\(A\)不会是对称、正定、投影矩阵。因为特征向量的正交性我们才能直接去看特征值的性质。

  3. 复习奇异值分解,\(A=U\varSigma V^T\):

    先求正交矩阵\(V\):\(A^TA=V\varSigma^TU^TU\varSigma V^T=V\left(\varSigma^T\varSigma\right)V^T\),所以\(V\)是矩阵\(A^TA\)的特征向量矩阵,而矩阵\(\varSigma^T\varSigma\)是矩阵\(A^TA\)的特征值矩阵,即\(A^TA\)的特征值为\(\sigma^2\)。

    接下来应该求正交矩阵\(U\):\(AA^T=U\varSigma^TV^TV\varSigma U^T=U\left(\varSigma^T\varSigma\right)U^T\),但是请注意,我们在这个式子中无法确定特征向量的符号,我们需要使用\(Av_i=\sigma_iu_i\),通过已经求出的\(v_i\)来确定\(u_i\)的符号(因为\(AV=U\varSigma\)),进而求出\(U\)。

    已知\(A=\bigg[u_1\ u_2\bigg]\begin{bmatrix}3&0\\0&2\end{bmatrix}\bigg[v_1\ v_2\bigg]^T\)

    从已知的\(\varSigma\)矩阵可以看出,\(A\)矩阵是非奇异矩阵,因为它没有零奇异值。另外,如果把\(\varSigma\)矩阵中的\(2\)改成\(-5\),则题目就不再是奇异值分解了,因为奇异值不可能为负;如果将\(2\)变为\(0\),则\(A\)是奇异矩阵,它的秩为\(1\),零空间为\(1\)维,\(v_2\)在其零空间中。

  4. \(A\)是正交对称矩阵,那么它的特征值具有什么特点

    首先,对于对称矩阵,有特征值均为实数;然后是正交矩阵,直觉告诉我们\(|\lambda|=1\)。来证明一下,对于\(Qx=\lambda x\),我们两边同时取模有\(\|Qx\|=|\lambda|\|x\|\),而正交矩阵不会改变向量长度,所以有\(\|x\|=|\lambda|\|x\|\),因此\(\lambda=\pm1\)。

    \(A\)是正定的吗?并不一定,因为特征向量可以取\(-1\)。

    \(A\)的特征值没有重复吗?不是,如果矩阵大于\(2\)阶则必定有重复特征值,因为只能取\(\pm1\)。

    \(A\)可以被对角化吗?是的,任何对称矩阵、任何正交矩阵都可以被对角化。

    \(A\)是非奇异矩阵吗?是的,正交矩阵都是非奇异矩阵。很明显它的特征值都不为零。

    证明\(P=\frac{1}{2}(A+I)\)是投影矩阵

    我们使用投影矩阵的性质验证,首先由于\(A\)是对称矩阵,则\(P\)一定是对称矩阵;接下来需要验证\(P^2=P\),也就是\(\frac{1}{4}\left(A^2+2A+I\right)=\frac{1}{2}(A+I)\)。来看看\(A^2\)是什么,\(A\)是正交矩阵则\(A^T=A^{-1}\),而\(A\)又是对称矩阵则\(A=A^T=A^{-1}\),所以\(A^2=I\)。带入原式有\(\frac{1}{4}(2A+2I)=\frac{1}{2}(A+I)\),得证。

    我们可以使用特征值验证,\(A\)的特征值可以取\(\pm1\),则\(A+I\)的特征值可以取\(0,2\),\(\frac{1}{2}(A+I)\)的特征值为\(0,1\),特征值满足投影矩阵且它又是对称矩阵,得证。

第三十四讲:左右逆和伪逆

前面我们涉及到的逆(inverse)都是指左、右乘均成立的逆矩阵,即\(A^{-1}A=I=AA^{-1}\)。在这种情况下,\(m\times n\)矩阵\(A\)满足\(m=n=rank(A)\),也就是满秩方阵。

左逆(left inserve)

记得我们在最小二乘一讲(第十六讲)介绍过列满秩的情况,也就是列向量线性无关,但行向量通常不是线性无关的。常见的列满秩矩阵\(A\)满足\(m>n=rank(A)\)。

列满秩时,列向量线性无关,所以其零空间中只有零解,方程\(Ax=b\)可能有一个唯一解(\(b\)在\(A\)的列空间中,此特解就是全部解,因为通常的特解可以通过零空间中的向量扩展出一组解集,而此时零空间只有列向量),也可能无解(\(b\)不在\(A\)的列空间中)。

另外,此时行空间为\(\mathbb{R}^n\),也正印证了与行空间互为正交补的零空间中只有列向量。

现在来观察\(A^TA\),也就是在\(m>n=rank(A)\)的情况下,\(n\times m\)矩阵乘以\(m\times n\)矩阵,结果为一个满秩的\(n\times n\)矩阵,所以\(A^TA\)是一个可逆矩阵。也就是说\(\underbrace{\left(A^TA\right)^{-1}A^T}A=I\)成立,而大括号部分的\(\left(A^TA\right)^{-1}A^T\)称为长方形矩阵\(A\)的左逆

\[A^{-1}_{left}=\left(A^TA\right)^{-1}A^T \]

顺便复习一下最小二乘一讲,通过关键方程\(A^TA\hat x=A^Tb\),\(A^{-1}_{left}\)被当做一个系数矩阵乘在\(b\)向量上,求得\(b\)向量投影在\(A\)的列空间之后的解\(\hat x=\left(A^TA\right)^{-1}A^Tb\)。如果我们强行给左逆左乘矩阵\(A\),得到的矩阵就是投影矩阵\(P=A\left(A^TA\right)^{-1}A^T\),来自\(p=A\hat x=A\left(A^TA\right)^{-1}A^T\),它将右乘的向量\(b\)投影在矩阵\(A\)的列空间中。

再来观察\(AA^T\)矩阵,这是一个\(m\times m\)矩阵,秩为\(rank(AA^T)=n<m\),也就是说\(AA^T\)是不可逆的,那么接下来我们看看右逆。

右逆(right inverse)

可以与左逆对称的看,右逆也就是研究\(m\times n\)矩阵\(A\)行满秩的情况,此时\(n>m=rank(A)\)。对称的,其左零空间中仅有零向量,即没有行向量的线性组合能够得到零向量。

行满秩时,矩阵的列空间将充满向量空间\(C(A)=\mathbb{R}^m\),所以方程\(Ax=b\)总是有解集,由于消元后有\(n-m\)个自由变量,所以方程的零空间为\(n-m\)维。

与左逆对称,再来观察\(AA^T\),在\(n>m=rank(A)\)的情况下,\(m\times n\)矩阵乘以\(n\times m\)矩阵,结果为一个满秩的\(m\times m\)矩阵,所以此时\(AA^T\)是一个满秩矩阵,也就是\(AA^T\)可逆。所以\(A\underbrace{A^T\left(AA^T\right)}=I\),大括号部分的\(A^T\left(AA^T\right)\)称为长方形矩阵的右逆

\[A^{-1}_{right}=A^T\left(AA^T\right) \]

同样的,如果我们强行给右逆右乘矩阵\(A\),将得到另一个投影矩阵\(P=A^T\left(AA^T\right)A\),与上一个投影矩阵不同的是,这个矩阵的\(A\)全部变为\(A^T\)了。所以这是一个能够将右乘的向量\(b\)投影在\(A\)的行空间中。

前面我们提及了逆(方阵满秩),并讨论了左逆(矩阵列满秩)、右逆(矩阵行满秩),现在看一下第四种情况,\(m\times n\)矩阵\(A\)不满秩的情况。

伪逆(pseudo inverse)

有\(m\times n\)矩阵\(A\),满足\(rank(A)\lt min(m,\ n)\),则

现在任取一个向量\(x\),乘上\(A\)后结果\(Ax\)一定落在矩阵\(A\)的列空间\(C(A)\)中。而根据维数,\(x\in\mathbb{R}^n,\ Ax\in\mathbb{R}^m\),那么我们现在猜测,输入向量\(x\)全部来自矩阵的行空间,而输出向量\(Ax\)全部来自矩阵的列空间,并且是一一对应的关系,也就是\(\mathbb{R}^n\)的\(r\)维子空间到\(\mathbb{R}^m\)的\(r\)维子空间的映射。

而矩阵\(A\)现在有这些零空间存在,其作用是将某些向量变为零向量,这样\(\mathbb{R}^n\)空间的所有向量都包含在行空间与零空间中,所有向量都能由行空间的分量和零空间的分量构成,变换将零空间的分量消除。但如果我们只看行空间中的向量,那就全部变换到列空间中了。

那么,我们现在只看行空间与列空间,在行空间中任取两个向量\(x,\ y\in C(A^T)\),则有\(Ax\neq Ay\)。所以从行空间到列空间,变换\(A\)是个不错的映射,如果限制在这两个空间上,\(A\)可以说“是个可逆矩阵”,那么它的逆就称作伪逆,而这个伪逆的作用就是将列空间的向量一一映射到行空间中。通常,伪逆记作\(A^+\),因此\(Ax=(Ax),\ y=A^+(Ay)\)。

现在我们来证明对于\(x,y\in C\left(A^T\right),\ x\neq y\),有\(Ax,Ay\in C(A),\ Ax\neq Ay\):

伪逆在统计学中非常有用,以前我们做最小二乘需要矩阵列满秩这一条件,只有矩阵列满秩才能保证\(A^TA\)是可逆矩阵,而统计中经常出现重复测试,会导致列向量线性相关,在这种情况下\(A^TA\)就成了奇异矩阵,这时候就需要伪逆。

接下来我们介绍如何计算伪逆\(A^+\):

其中一种方法是使用奇异值分解,\(A=U\varSigma V^T\),其中的对角矩阵型为\(\varSigma=\left[\begin{array}{c c c|c}\sigma_1&&&\\&\ddots&&\\&&\sigma_2&\\\hline&&&\begin{bmatrix}0\end{bmatrix}\end{array}\right]\),对角线非零的部分来自\(A^TA,\ AA^T\)比较好的部分,剩下的来自左/零空间。

我们先来看一下\(\varSigma\)矩阵的伪逆是多少,这是一个\(m\times n\)矩阵,\(rank(\varSigma)=r\),\(\varSigma^+=\left[\begin{array}{c c c|c}\frac{1}{\sigma_1}&&&\\&\ddots&&\\&&\frac{1}{\sigma_r}&\\\hline&&&\begin{bmatrix}0\end{bmatrix}\end{array}\right]\),伪逆与原矩阵有个小区别:这是一个\(n\times m\)矩阵。则有\(\varSigma\varSigma^+=\left[\begin{array}{c c c|c}1&&&\\&\ddots&&\\&&1&\\\hline&&&\begin{bmatrix}0\end{bmatrix}\end{array}\right]_{m\times m}\),\(\varSigma^+\varSigma=\left[\begin{array}{c c c|c}1&&&\\&\ddots&&\\&&1&\\\hline&&&\begin{bmatrix}0\end{bmatrix}\end{array}\right]_{n\times n}\)。

观察\(\varSigma\varSigma^+\)和\(\varSigma^+\varSigma\)不难发现,\(\varSigma\varSigma^+\)是将向量投影到列空间上的投影矩阵,而\(\varSigma^+\varSigma\)是将向量投影到行空间上的投影矩阵。我们不论是左乘还是右乘伪逆,得到的不是单位矩阵,而是投影矩阵,该投影将向量带入比较好的空间(行空间和列空间,而不是左/零空间)。

接下来我们来求\(A\)的伪逆:

\[A^+=V\varSigma^+U^T \]

第三十五讲:期末复习

依然是从以往的试题入手复习知识点。

  1. 已知\(m\times n\)矩阵\(A\),有\(Ax=\begin{bmatrix}1\\0\\0\end{bmatrix}\)无解;\(Ax=\begin{bmatrix}0\\1\\0\end{bmatrix}\)仅有唯一解,求关于\(m,n,rank(A)\)的信息。

    首先,最容易判断的是\(m=3\);而根据第一个条件可知,矩阵不满秩,有\(r<m\);根据第二个条件可知,零空间仅有零向量,也就是矩阵消元后没有自由变量,列向量线性无关,所以有\(r=n\)。

    综上,有\(m=3>n=r\)。

    根据所求写出一个矩阵\(A\)的特例:\(A=\begin{bmatrix}0&0\\1&0\\0&1\end{bmatrix}\)。

    \(\det A^TA\stackrel{?}{=}\det AA^T\):不相等,因为\(A^TA\)可逆而\(AA^T\)不可逆,所以行列式不相等。(但是对于方阵,\(\det AB=\det BA\)恒成立。)

    \(A^TA\)可逆吗?是,因为\(r=n\),矩阵列向量线性无关,即列满秩。

    \(AA^T\)正定吗?否,因为\(AA^T\)是\(3\times n\)矩阵与\(n\times 3\)矩阵之积,是一个三阶方阵,而\(AA^T\)秩为\(2\),所以不是正定矩阵。(不过\(AA^T\)一定是半正定矩阵。)

    求证\(A^Ty=c\)至少有一个解:因为\(A\)的列向量线性无关,所以\(A^T\)的行向量线性无关,消元后每行都有主元,且总有自由变量,所以零空间中有非零向量,零空间维数是\(m-r\)(可以直接从\(\dim N\left(A^T\right)=m-r\)得到结论)。

  2. 设\(A=\Bigg[v_1\ v_2\ v_3\Bigg]\),对于\(Ax=v_1-v_2+v_3\),求\(x\)。

    按列计算矩阵相乘,有\(x=\begin{bmatrix}1\\-1\\1\end{bmatrix}\)。

    若Ax=v_1-v_2+v_3=0,则解是唯一的吗?为什么。如果解释唯一的,则零空间中只有零向量,而在此例中\(x=\begin{bmatrix}1\\-1\\1\end{bmatrix}\)就在零空间中,所以解不唯一。

    若\(v_1,v_2,v_3\)是标准正交向量,那么怎样的线性组合\(c_1v_1+c_2v_2\)能够最接近\(v_3\)?此问是考察投影概念,由于是正交向量,所以只有\(0\)向量最接近\(v_3\)。

  3. 矩阵\(A=\begin{bmatrix}.2&.4&.3\\.4&.2&.3\\.4&.4&.4\end{bmatrix}\),求稳态。

    这是个马尔科夫矩阵,前两之和为第三列的两倍,奇异矩阵总有一个特征值为\(0\),而马尔科夫矩阵总有一个特征值为\(1\),剩下一个特征值从矩阵的迹得知为\(-.2\)。

    再看马尔科夫过程,设从\(u(0)\)开始,\(u_k=A^ku_0, u_0=\begin{bmatrix}0\\10\\0\end{bmatrix}\)。先代入特征值\(\lambda_1=0,\ \lambda_2=1,\ \lambda_3=-.2\)查看稳态\(u_k=c_1\lambda_1^kx_1+c_2\lambda_2^kx_2+c_3\lambda_3^kx_3\),当\(k\to\infty\),第一项与第三项都会消失,剩下\(u_\infty=c_2x_2\)。

    到这里我们只需求出\(\lambda_2\)对应的特征向量即可,带入特征值求解\((A-I)x=0\),有\(\begin{bmatrix}-.8&.4&.3\\.4&-.8&.3\\.4&.4&-.6\end{bmatrix}\begin{bmatrix}?\\?\\?\end{bmatrix}=\begin{bmatrix}0\\0\\0\end{bmatrix}\),可以消元得,也可以直接观察得到\(x_2=\begin{bmatrix}3\\3\\4\end{bmatrix}\)。

    剩下就是求\(c_2\)了,可以通过\(u_0\)一一解出每个系数,但是这就需要解出每一个特征值。另一种方法,我们可以通过马尔科夫矩阵的特性知道,对于马尔科夫过程的每一个\(u_k\)都有其分量之和与初始值分量之和相等,所以对于\(x_2=\begin{bmatrix}3\\3\\4\end{bmatrix}\),有\(c_2=1\)。所以最终结果是\(u_\infty=\begin{bmatrix}3\\3\\4\end{bmatrix}\)。

  4. 对于二阶方阵,回答以下问题:

    求投影在直线\(a=\begin{bmatrix}4\\-3\end{bmatrix}\)上的投影矩阵:应为\(P=\frac{aa^T}{a^Ta}\)。

    已知特征值\(\lambda_1=2,\ x_1=\begin{bmatrix}1\\2\end{bmatrix}\quad \lambda_2=3,\ x_2=\begin{bmatrix}2\\1\end{bmatrix}\)求原矩阵\(A\):从对角化公式得\(A=S\Lambda S^{-1}=\begin{bmatrix}1&2\\2&1\end{bmatrix}\begin{bmatrix}0&0\\0&3\end{bmatrix}\begin{bmatrix}1&2\\2&1\end{bmatrix}^{-1}\),解之即可。

    \(A\)是一个实矩阵,且对任意矩阵\(B\),\(A\)都不能分解成\(A=B^TB\),给出\(A\)的一个例子:我们知道\(B^TB\)是对称的,所以给出一个非对称矩阵即可。
    矩阵\(A\)有正交的特征向量,但不是对称的,给出一个\(A\)的例子:我们在三十三讲提到过,反对称矩阵,因为满足\(AA^T=A^TA\)而同样具有正交的特征向量,所以有\(A=\begin{bmatrix}0&1\\-1&0\end{bmatrix}\)或旋转矩阵\(\begin{bmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{bmatrix}\),这些矩阵都具有复数域上的正交特征向量组。

  5. 最小二乘问题,因为时间的关系直接写出计算式和答案,\(\begin{bmatrix}1&0\\1&1\\1&2\end{bmatrix}\begin{bmatrix}C\\D\end{bmatrix}=\begin{bmatrix}3\\4\\1\end{bmatrix}(Ax=b)\),解得\(\begin{bmatrix}\hat C\\\hat D\end{bmatrix}=\begin{bmatrix}\frac{11}{3}\\-1\end{bmatrix}\)。

    求投影后的向量\(p\):向量\(p\)就是向量\(b\)在矩阵\(A\)列空间中的投影,所以\(p=\begin{bmatrix}p_1\\p_2\\p_3\end{bmatrix}=\begin{bmatrix}1&0\\1&1\\1&2\end{bmatrix}\begin{bmatrix}\hat C\\\hat D\end{bmatrix}\)。

    求拟合直线的图像:\(x=0,1,2\)时\(y=p_1,p_2,p_2\)所在的直线的图像,\(y=\hat C+\hat Dx\)即\(y=\frac{11}{3}-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()

png

plt.close(fig)

标签:特征值,begin,end,灰飞烟灭,矩阵,bmatrix,06,MIT18,lambda
来源: https://www.cnblogs.com/limbercode/p/16154859.html