线性代数杂记

部分内容来自 [MIT 18.06 Linear Algebra, Spring 2005, Gilbert Strang]。

估计 \(n\) 阶方阵 \(\mathbf{A}\) 的 \(\mathbf{LU}\) 分解的计算量

1. 第一步，将 \(a_{11}\) 作为主元，需要的运算量约为 \(n^2\)。

\[ \begin{pmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \dots & a_{nn} \end{pmatrix} \underrightarrow{消元} \begin{pmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ 0 & a_{22} & \dots & a_{2n} \\ 0 & \vdots & \ddots & \vdots \\ 0 & a_{n2} & \dots & a_{nn} \end{pmatrix} \]

2. 以此类推，接下来每一步计算量约为 \((n-1)^2, (n-2)^2, \dots, 2^2, 1^2\)。

3. 则将 \(\mathbf{A}\) 变换为 \(\mathbf{LU}\) 的总运算量应为 \(O(n^2+(n-1)^2+\dots+2^2+1^2)\)，即 \(O(\frac{n^3}{3})\)。

四个基本子空间

对于 \(m\times n\) 矩阵 \(\mathbf{A}\)，\(\mathrm{rank\ }\mathbf{A}=r\) 有：

• 行空间 \(C(\mathbf{A}^\top)\in\mathbb{R}^n,\ \mathrm{dim\ }C(\mathbf{A}^\top)=r\)，等价阶梯形的非零行可组成行空间的一组基。

• 零空间 \(N(\mathbf{A})\in\mathbb{R}^n,\ \mathrm{dim\ }N(\mathbf{A})=n-r\)，自由元所在的列即可组成零空间的一组基。

• 列空间 \(C(\mathbf{A})\in\mathbb{R}^m,\ \mathrm{dim\ }C(\mathbf{A})=r\)，主元所在的列即可组成列空间的一组基。

• 左零空间 \(N(\mathbf{A}^\top)\in\mathbb{R}^m,\ \mathrm{dim\ }N(\mathbf{A}^\top)=m-r\)，采用 Gauss-Jordan 消元，将增广矩阵 \([\mathbf{A}_{m \times n}\mid \mathbf{I}_{m \times m}]\) 中 \(\mathbf{A}\) 的部分划为简化行阶梯形式 \([\mathbf{r}_{m \times n}\mid \mathbf{E}_{m \times m}]\)，此时矩阵 \(\mathbf{E}\) 会将所有的行变换记录下来，\(\mathbf{EA}=\mathbf{r}\)，\(\mathbf{r}\) 为简化阶梯形。

对于左零空间，有 \(\mathbf{A}^\top\mathbf{y}=\mathbf{0} \to (\mathbf{A}^\top\mathbf{y})^\top=\mathbf{0}^\top\to \mathbf{y}^\top\mathbf{A}=\mathbf{0}^\top\)，因此得名。

子空间投影

有 \(n\) 维列向量 \(\mathbf{A}, \mathbf{b}\)，作 \(\mathbf{b}\) 在 \(\mathbf{A}\) 上的投影 \(\mathbf{p}\)，得 \(\mathbf{E}=\mathbf{b}-\mathbf{p}\)，记 \(\mathbf{p}=x\mathbf{A}=\mathbf{A}\cdot x\)。由垂直关系得

\[ \begin{gather} \mathbf{A}^\top\mathbf{E}=\mathbf{A}^\top(\mathbf{b}-\mathbf{p})=\mathbf{A}^\top(\mathbf{b}-x\mathbf{A})=0 \\ x\mathbf{A}^\top\mathbf{A}=\mathbf{A}^\top\mathbf{b} \\ x=\frac{\mathbf{A}^\top\mathbf{b}}{\mathbf{A}^\top\mathbf{A}} \\ \mathbf{p}=\mathbf{A}\frac{\mathbf{A}^\top\mathbf{b}}{\mathbf{A}^\top\mathbf{A}} \end{gather} \]

从上面的式子可以看出，如果将 \(\mathbf{b}\) 变为 \(2\mathbf{b}\) 则 \(\mathbf{p}\) 也会翻倍，如果将 \(\mathbf{A}\) 变为 \(2\mathbf{A}\) 则 \(\mathbf{p}\) 不变。

设投影矩阵 \(\mathbf{p}\) 是一个 \(n\times n\) 的矩阵，\(\mathbf{p}=\frac{\mathbf{aa}^{\mathrm{T}}}{\mathbf{A}^\top\mathbf{A}}\)，则可以说投影矩阵作用于某个向量后，得到其投影向量 \(\mathbf{p}=\mathbf{Pb}\)。

投影矩阵 \(\mathbf{p}\) 的列空间 \(C(\mathbf{p})\) 是一条通过 \(\mathbf{A}\) 的直线，而 \(\mathrm{rank\ }\mathbf{p}=1\)(一列乘以一行：\(\mathbf{aa}^{\mathrm{T}}\)，而这一列向量 \(\mathbf{A}\) 是该矩阵的基)。

投影矩阵的性质：

• \(\mathbf{p}=\mathbf{p}^{\mathrm{T}}\)，即投影矩阵是一个对称矩阵。
• 如果对一个向量做两次投影，即 \(\mathbf{PPb}\)，则其结果仍然与 \(\mathbf{Pb}\) 相同，也就是\(\mathbf{p}^2=\mathbf{p}\)。

为什么我们需要投影？因为有时候 \(\mathbf{Ax}=\mathbf{b}\) 无解，我们只能求出最接近的那个解。

\(\mathbf{Ax}\) 总是在 \(\mathbf{A}\) 的列空间中，而 \(\mathbf{b}\) 却不一定，这是问题所在，所以我们可以将 \(\mathbf{b}\) 变为 \(\mathbf{A}\) 的列空间中最接近的那个向量，即将无解的 \(\mathbf{Ax}=\mathbf{b}\) 变为求有解的\(\mathbf{A}\hat{\mathbf{x}}=\mathbf{p}\)(\(\mathbf{p}\) 是 \(\mathbf{b}\) 在 \(\mathbf{A}\) 的列空间中的投影，\(\hat{\mathbf{x}}\) 不再是那个不存在的 \(\mathbf{x}\)，而是最接近的解)。

现在来看 \(\mathbb{R}^3\) 中的情形，将向量 \(\mathbf{b}\) 投影在平面 \(\mathbf{A}\) 上。同样的，\(\mathbf{p}\) 是向量 \(\mathbf{b}\) 在平面 \(\mathbf{A}\) 上的投影，\(\mathbf{E}\) 是垂直于平面 \(\mathbf{A}\) 的向量，即 \(\mathbf{b}\) 在平面 \(\mathbf{A}\) 法方向的分量。

设平面 \(\mathbf{A}\) 的一组基为 \(\mathbf{A}_1, \mathbf{A}_2\)，则投影向量 \(\mathbf{p}=\hat{\mathbf{x}_1}\mathbf{A}_1+\hat{\mathbf{x}_2}\mathbf{A}_2\)，我们更倾向于写作 \(\mathbf{p}=\mathbf{A}\hat{\mathbf{x}}\)，这里如果我们求出 \(\hat{\mathbf{x}}\)，则该解就是无解方程组最近似的解。

现在问题的关键在于找 \(\mathbf{E}=\mathbf{b}-\mathbf{A}\hat{\mathbf{x}}\)，使它垂直于平面，因此我们得到两个方程

\[ \begin{cases} \mathbf{A}_1^\top(\mathbf{b}-\mathbf{A}\hat{\mathbf{x}})=0 \\ \mathbf{A}_2^\top(\mathbf{b}-\mathbf{A}\hat{\mathbf{x}})=0 \end{cases}, \]

将方程组写成矩阵形式

\[ \begin{pmatrix} \mathbf{A}_1^\top\\ \mathbf{A}_2^\top \end{pmatrix} \left(\mathbf{b}-\mathbf{A}\hat{\mathbf{x}}\right)=\begin{pmatrix} 0 \\ 0 \end{pmatrix}, \]

即

\[ \mathbf{A}^\top\left(\mathbf{b}-\mathbf{A}\hat{\mathbf{x}}\right)=0. \]

比较该方程与 \(\mathbb{R}^2\) 中的投影方程，发现只是向量 \(\mathbf{A}\) 变为矩阵 \(\mathbf{A}\) 而已，本质上就是 \(\mathbf{A}^\top\mathbf{E}=0\)。所以，\(\mathbf{E}\) 在 \(\mathbf{A}^\top\)的零空间中(\(\mathbf{E}\in N(\mathbf{A}^\top)\))。我们知道，左零空间 \(\bot\) 列空间，则有 \(\mathbf{E}\bot C(\mathbf{A})\)，与我们设想的一致。

再化简方程得 \(\mathbf{A}^\top\mathbf{A}\hat{\mathbf{x}}=\mathbf{A}^\top\mathbf{b}\)，比较在 \(\mathbb{R}^2\) 中的情形，\(\mathbf{A}^\top\mathbf{A}\) 是一个数字，而 \(\mathbf{A}^\top\mathbf{A}\) 是一个 \(n\) 阶方阵，而解出的 \(\hat{\mathbf{x}}\) 可以看做两个数字的比值。现在在 \(\mathbb{R}^3\) 中，我们需要再次考虑：什么是 \(\hat{\mathbf{x}}\)？投影是什么？投影矩阵又是什么？

• 第一个问题：\(\hat{\mathbf{x}}=(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top\mathbf{b}\)；
• 第二个问题：\(\mathbf{p}=\mathbf{A}\hat{\mathbf{x}}=\underline{\mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top}\mathbf{b}\)，回忆在 \(\mathbb{R}^2\) 中的情形，下划线部分就是原来的 \(\frac{\mathbf{aa}^\top}{\mathbf{A}^\top\mathbf{A}}\)；
• 第三个问题：易看出投影矩阵就是下划线部分 \(\mathbf{p}=\mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top\)。

这里还需要注意一个问题，\(\mathbf{p}=\mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top\) 是不能继续化简的，因为这里的 \(\mathbf{A}\) 并不是一个可逆方阵。 也可以换一种思路，如果 \(\mathbf{A}\) 是一个 \(n\) 阶可逆方阵，则 \(\mathbf{A}\) 的列空间是整个 \(\mathbb{R}^n\) 空间，于是 \(\mathbf{b}\) 在 \(\mathbb{R}^n\) 上的投影矩阵确实变为了 \(\mathbf{I}\)，因为 \(\mathbf{b}\) 已经在空间中了，其投影不再改变。

证明投影矩阵 \(\mathbf{p}\) 的性质

• 证明 \(\mathbf{p}=\mathbf{p}^\top\)：

\[ \mathbf{p}^\top=\left[\mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top\right]^\top=\mathbf{A}\left[(\mathbf{A}^\top\mathbf{A})^{-1}\right]^\top\mathbf{A}^\top, \]

而 \(\mathbf{A}^\top\mathbf{A}\) 是对称的，所以其逆也是对称的，于是有

\[ \mathbf{A}\left[(\mathbf{A}^\top\mathbf{A})^{-1}\right]^\top\mathbf{A}^\top=\mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top=\mathbf{p}, \]

得证。

• 证明 \(\mathbf{p}^2=\mathbf{p}\)：

\[ \begin{align} \mathbf{p}^2 & = \left[\mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top\right]\left[\mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top\right] \\ & = \mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\left[(\mathbf{A}^\top\mathbf{A})(\mathbf{A}^\top\mathbf{A})^{-1}\right]\mathbf{A}^\top\\ & = \mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top\\ & = \mathbf{p} \end{align} \]

得证。

举两个极端的例子：

• 如果 \(\mathbf{b}\in C(\mathbf{A})\)，则 \(\mathbf{Pb}=\mathbf{b}\)；
• 如果 \(\mathbf{b}\bot C(\mathbf{A})\)，则 \(\mathbf{Pb}=0\)。

一般情况下，\(\mathbf{b}\) 将会有一个垂直于 \(\mathbf{A}\) 的分量，有一个在 \(\mathbf{A}\) 列空间中的分量，投影的作用就是去掉垂直分量而保留列空间中的分量。

在第一个极端情况中，如果 \(\mathbf{b}\in C(\mathbf{A})\) 则有 \(\mathbf{b}=\mathbf{Ax}\)。带入投影矩阵\(\mathbf{p}=\mathbf{Pb}=\mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top\mathbf{Ax}=\mathbf{Ax}\)，得证。

在第二个极端情况中，如果 \(\mathbf{b}\bot C(\mathbf{A})\) 则有 \(\mathbf{b}\in N(\mathbf{A}^\top)\)，即 \(\mathbf{A}^\top\mathbf{b}=0\)。则 \(\mathbf{p}=\mathbf{Pb}=\mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top\mathbf{b}=0\)，得证。

向量 \(\mathbf{b}\) 投影后，有 \(\mathbf{b}=\mathbf{E}+\mathbf{p}, \mathbf{p}=\mathbf{Pb}, \mathbf{E}=(\mathbf{I}-\mathbf{p})\mathbf{b}\)，这里的 \(\mathbf{p}\) 是 \(\mathbf{b}\) 在 \(C(\mathbf{A})\) 中的分量，而 \(\mathbf{E}\) 是 \(\mathbf{b}\) 在 \(N(\mathbf{A}^\top)\) 中的分量。

最小二乘法

接下来看投影的经典应用案例：最小二乘法拟合直线(least squares fitting by a line)。

我们需要找到距离图中三个点 \((1, 1), (2, 2), (3, 2)\) 偏差最小的直线：\(b=C+Dt\)。

根据条件可以得到方程组

\[ \begin{cases} C+D & =1 \\ C+2D & =2 \\ C+3D & =2 \end{cases} \]

写作矩阵形式

\[ \begin{pmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \end{pmatrix} \begin{pmatrix} C \\ D \end{pmatrix} = \begin{pmatrix} 1 \\ 2 \\ 2 \end{pmatrix} \]

也就是我们的 \(\mathbf{Ax}=\mathbf{b}\)，很明显方程组无解。但是 \(\mathbf{A}^\top\mathbf{A}\hat{\mathbf{x}}=\mathbf{A}^\top\mathbf{b}\) 有解，是最小二乘法的核心方程。

我们需要在 \(\mathbf{b}\) 的三个分量上都增加某个误差 \(\mathbf{E}\)，使得三点能够共线，同时使得 \(e_1^2+e_2^2+e_3^2\) 最小，找到拥有最小平方和的解(即最小二乘)，即 \(\lvert\mathbf{Ax}-\mathbf{b}\rvert^2=\lvert\mathbf{E}\rvert^2\) 最小。此时向量 \(\mathbf{b}\) 变为向量 \(\mathbf{p}=\begin{pmatrix} p_1 \\ p_2 \\ p_3 \end{pmatrix}\)(在方程组有解的情况下，\(\mathbf{Ax}-\mathbf{b}=0\)，即 \(\mathbf{b}\) 在 \(\mathbf{A}\) 的列空间中，误差 \(\mathbf{E}\) 为零。)我们现在做的运算也称作线性回归，使用误差的平方和作为测量总误差的标准。

注：如果有另一个点，如 \((0, 100)\)，在本例中该点明显距离别的点很远，最小二乘将很容易被离群的点影响，通常使用最小二乘时会去掉明显离群的点。

现在我们尝试解出 \(\hat{\mathbf{x}}=\begin{pmatrix} \hat{C} \\ \hat{D} \end{pmatrix}\) 与 \(\mathbf{p}=\begin{pmatrix} p_1 \\ p_2 \\ p_3 \end{pmatrix}\)。

\[ \mathbf{A}^\top\mathbf{A}\hat{\mathbf{x}}=\mathbf{A}^\top\mathbf{b} \longrightarrow \begin{pmatrix} 3 & 6 \\ 6 & 14 \end{pmatrix} \begin{pmatrix}\hat{C} \\ \hat{D} \end{pmatrix}= \begin{pmatrix} 5 \\ 11 \end{pmatrix} \]

写作方程形式为 \(\begin{cases} 3\hat{C}+16\hat{D} & =5 \\ 6\hat{C}+14\hat{D} & =11 \end{cases}\)，也称作正规方程组(normal equations)。

回顾前面提到的“使得误差最小”的条件，\(e_1^2+e_2^2+e_3^2=(C+D-1)^2+(C+2D-2)^2+(C+3D-2)^2\)，对 \(C, D\) 分别求偏导，再令求得的偏导式为零，也有正规方程组，解得 \(\hat{C}=\frac{2}{3}, \hat{D}=\frac{1}{2}\)，则“最佳直线”为 \(y=\frac{2}{3}+\frac{1}{2}t\)，带回原方程组解得 \(p_1=\frac{7}{6}, p_2=\frac{5}{3}, p_3=\frac{13}{6}\)，即 \(e_1=-\frac{1}{6}, e_2=\frac{1}{3}, e_3=-\frac{1}{6}\)，于是我们得到 \(\mathbf{p}=\begin{pmatrix} \frac{7}{6} \\ \frac{5}{3} \\ \frac{13}{6} \end{pmatrix}, \mathbf{E}=\begin{pmatrix} -\frac{1}{6} \\ \frac{1}{3} \\ -\frac{1}{6} \end{pmatrix}\)，易看出 \(\mathbf{b}=\mathbf{p}+\mathbf{E}\)，且 \(\mathbf{p}\cdot\mathbf{E}=0\)。误差向量 \(\mathbf{E}\) 不仅垂直于投影向量 \(\mathbf{p}\)，它还垂直于列空间 \(\begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}, \begin{pmatrix} 1 \\ 2 \\ 3 \end{pmatrix}\)。

接下来我们观察 \(\mathbf{A}^\top\mathbf{A}\)，如果 \(\mathbf{A}\) 的各列线性无关，求证 \(\mathbf{A}^\top\mathbf{A}\) 是可逆矩阵。

先假设 \(\mathbf{A}^\top\mathbf{Ax}=0\)，两边同时乘以 \(\mathbf{x}^\top\)，有 \(\mathbf{x}^\top\mathbf{A}^\top\mathbf{Ax}=0\)，即 \((\mathbf{Ax})^\top(\mathbf{Ax})=0\)。一个矩阵乘其转置结果为零，则这个矩阵也必须为零(\((\mathbf{Ax})^\top(\mathbf{Ax})\) 相当于 \(\mathbf{Ax}\) 长度的平方)。则 \(\mathbf{Ax}=\mathbf{0}\)，结合题设中的“\(\mathbf{A}\) 的各列线性无关”，可知 \(\mathbf{x}=\mathbf{0}\)，也就是 \(\mathbf{A}^\top\mathbf{A}\) 的零空间中有且只有零向量，得证。

\(\mathbf{A}\) 的幂

求 \(\mathbf{u}_{k+1}=\mathbf{Au}_k\)

本例是一个一阶差分方程组(first-order system)，有 \(\mathbf{u}_1=\mathbf{Au}_0\)，\(\mathbf{u}_2=\mathbf{A}^2\mathbf{u}_0\)，\(\mathbf{u}_k=\mathbf{A}^k\mathbf{u}_0\)……要解此方程，需要将 \(\mathbf{u}_0\) 展开为矩阵 \(\mathbf{A}\) 特征向量的线性组合，即 \(\mathbf{u}_0=\Bigg[x_1x_2\dots x_n\Bigg]\begin{pmatrix}c_1\\c_2\\\vdots\\c_n\end{pmatrix}=\mathbf{Sc}\)，于是 \(\mathbf{Au}_0=\mathbf{S\Lambda S}^{-1}\mathbf{u}_0=\mathbf{S\Lambda S}^{-1}\mathbf{Sc}=\mathbf{S\Lambda c}\)。

接下来看斐波那契数列：\(F_{k+2}=F_{k+1}+F_{k}\)，但这不是 \(\mathbf{u}_{k+1}=\mathbf{Au}_{k}\) 的形式。令 \(\mathbf{u}_{k}=\begin{pmatrix}F_{k+1}\\F_{k}\end{pmatrix}\)，由 \(\begin{cases}F_{k+2}&=F_{k+1}+F_{k}\\F_{k+1}&=F_{k+1}\end{cases}\)，得 \(\begin{pmatrix}F_{k+2}\\F_{k+1}\end{pmatrix}=\begin{pmatrix}1&1\\1&0\end{pmatrix}\begin{pmatrix}F_{k+1}\\F_{k}\end{pmatrix}\)，于是 \(\mathbf{u}_{k+1}=\mathbf{Au}_{k}, \mathbf{A}=\begin{pmatrix}1&1\\1&0\end{pmatrix}\)。我们把二阶标量方程(second-order scalar problem)转化为了一阶向量方程组(first-order system)。

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

• 用求根公式解得特征值 \(\begin{cases}\lambda_1=\frac{1}{2}\left(1+\sqrt{5}\right)\approx{1.618}\\\lambda_2=\frac{1}{2}\left(1-\sqrt{5}\right)\approx{-0.618}\end{cases}\)，得到两个不同的特征值，也就是有两个线性无关的特征向量，则该矩阵可以被对角化。

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

\(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\) 。可以看出，这种问题与求解 \(\mathbf{Ax}=\mathbf{b}\) 不同，这是一个动态的问题，\(\mathbf{A}\) 的幂在不停的增长，而问题的关键就是这些特征值。

• 继续求解特征向量，\(A-\lambda I=\begin{pmatrix}1-\lambda&1\\1&1-\lambda\end{pmatrix}\)，因为有根式且矩阵只有二阶，我们直接观察 \(\begin{pmatrix}1-\lambda&1\\1&1-\lambda\end{pmatrix}\begin{pmatrix}?\\?\end{pmatrix}=0\)，由于 \(\lambda^2-\lambda-1=0\)，则其特征向量为 \(\begin{pmatrix}\lambda\\1\end{pmatrix}\)，即 \(x_1=\begin{pmatrix}\lambda_1\\1\end{pmatrix}, x_2=\begin{pmatrix}\lambda_2\\1\end{pmatrix}\)。

最后，计算初始项 \(\mathbf{u}_0=\begin{pmatrix}F_1\\F_0\end{pmatrix}=\begin{pmatrix}1\\0\end{pmatrix}\)，解得 \(c_1=\frac{\sqrt{5}}{5}, c_2=-\frac{\sqrt{5}}{5}\)。

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

• 再下来就是套用公式，即 \(\mathbf{A}\) 的 \(k\) 次方表达式 \(\mathbf{A}^k=\mathbf{S\Lambda}^k\mathbf{S}^{-1}\)，则有 \(\mathbf{u}_{99}=\mathbf{Au}_{98}=\dots=\mathbf{A}^{99}\mathbf{u}_{0}=\mathbf{S\Lambda}^{99}\mathbf{S}^{-1}\mathbf{Sc}=\mathbf{S\Lambda}^{99}\mathbf{c}\)，代入特征值和特征向量得 \(\mathbf{u}_{99}=\begin{pmatrix}F_{100}\\F_{99}\end{pmatrix}=\begin{pmatrix}\frac{1+\sqrt{5}}{2}&\frac{1-\sqrt{5}}{2}\\1&1\end{pmatrix}\begin{pmatrix}\left(\frac{1+\sqrt{5}}{2}\right)^{99}&0\\0&\left(\frac{1-\sqrt{5}}{2}\right)^{99}\end{pmatrix}\begin{pmatrix}\frac{\sqrt{5}}{5}\\-\frac{\sqrt{5}}{5}\end{pmatrix}=\begin{pmatrix}c_1\lambda_1^{100}+c_2\lambda_2^{100}\\c_1\lambda_1^{99}+c_2\lambda_2^{99}\end{pmatrix}\)，最终结果为 \(F_{100}=c_1\lambda_1^{100}+c_2\lambda_2^{100}\)。

• 原式的通解为 \(u_k=c_1\lambda^kx_1+c_2\lambda^kx_2\)。

微分方程和 \(\mathrm{e}^{\mathbf{A}t}\)

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

设方程组 \(\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}\)，其系数矩阵为 \(\mathbf{A}=\begin{pmatrix}-1&2\\1&-2\end{pmatrix}\)，且 \(\mathbf{u}(0)=\begin{pmatrix}u_1\\u_2\end{pmatrix}=\begin{pmatrix}1\\0\end{pmatrix}\)。

• 初始条件可以看做在开始时，一切都在 \(u_1\) 中，但随着时间的推移，有 \(\frac{\mathrm{d}u_2}{\mathrm{d}t}>0\)，因为 \(u_1\) 项初始为正，\(u_1\) 中的部分会流向 \(u_2\)。

• 第一步需要找到特征值与特征向量。\(\mathbf{A}\) 是一个奇异矩阵，所以一个特征值是 \(\lambda_1=0\)，另一个可以从迹 \(\mathrm{tr}(A)=-3\) 得到。

\(t\to\infty\) 时，\(\mathrm{e}^{-3t}\to 0\)，即特征值 \(\lambda_2=-3\) 将会逐渐消失；另一项 \(\mathrm{e}^{0t}=1\)，其并不随时间而改变。通常含有特征值 0 的矩阵会随着时间的推移达到稳态。

• 求特征向量。\(\lambda_1=0\) 时，求 \(\mathbf{A}\) 的零空间 \(\mathbf{x}_1=\begin{pmatrix}2\\1\end{pmatrix}\)；\(\lambda_2=-3\) 时，求 \(\mathbf{A}+3I\) 的零空间，\(\begin{pmatrix}2&2\\1&1\end{pmatrix}\) 的零空间为 \(\mathbf{x}_2=\begin{pmatrix}1\\-1\end{pmatrix}\)。

• 则方程组的通解为：\(\mathbf{u}(t)=c_1\mathrm{e}^{\lambda_1t}x_1+c_2\mathrm{e}^{\lambda_2t}x_2\)，通解的前后两部分都是该方程组的纯解，即方程组的通解就是两个与特征值、特征向量相关的纯解的线性组合。我们来验证一下，比如取 \(\mathbf{u}=\mathrm{e}^{\lambda_1t}x_1\) 带入 \(\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}=\mathbf{Au}\)，对 \(t\) 求导得到 \(\lambda_1\mathrm{e}^{\lambda_1t}x_1=\mathbf{A}\mathrm{e}^{\lambda_1t}x_1\)，化简得 \(\lambda_1x_1=\mathbf{A}x_1\)。

对比解 \(\mathbf{u}_{k+1}=\mathbf{Au}_k\) 时得到的 \(\mathbf{u}_k=c_1\lambda^kx_1+c_2\lambda^kx_2\)。

• 继续求 \(c_1,c_2\)，\(\mathbf{u}(t)=c_1\cdot 1\cdot\begin{pmatrix}2\\1\end{pmatrix}+c_2\cdot\mathrm{e}^{-3t}\cdot\begin{pmatrix}1\\-1\end{pmatrix}\)，已知 \(t=0\) 时，\(\begin{pmatrix}1\\0\end{pmatrix}=c_1\begin{pmatrix}2\\1\end{pmatrix}+c_2\begin{pmatrix}1\\-1\end{pmatrix}\)(\(\mathbf{Sc}=\mathbf{u}(0)\))，所以 \(c_1=\frac{1}{3}, c_2=\frac{1}{3}\)。

• 于是我们写出最终结果，\(\mathbf{u}(t)=\frac{1}{3}\begin{pmatrix}2\\1\end{pmatrix}+\frac{1}{3}\mathrm{e}^{-3t}\begin{pmatrix}1\\-1\end{pmatrix}\)。

稳定性：这个流动过程从 \(\mathbf{u}(0)=\begin{pmatrix}1\\0\end{pmatrix}\) 开始，初始值 1 的一部分流入初始值 0 中，最终达到稳态 \(\mathbf{u}(\infty)=\begin{pmatrix}\frac{2}{3}\\\frac{1}{3}\end{pmatrix}\)。所以，要使得 \(\mathbf{u}(t)\to 0\)，则 需要负的特征值。但如果特征值为复数呢？如 \(\lambda=-3+6i\)，我们来计算 \(\lvert \mathrm{e}^{(-3+6i)t}\rvert\)，其中 \(\lvert \mathrm{e}^{6it}\rvert\) 这部分等于 \(\lvert \cos 6t+i\sin 6t\rvert =1\)，模为 \(\cos^2\alpha+\sin^2\alpha=1\)，虚部在单位圆上转，所以只有实数部分才是重要的，即 需要实部为负数的特征值。实部会决定最终结果趋近于 0 或 \(\infty\)，虚部不过是一些小杂音。

收敛态：其中一个特征值实部为 0，而其他特征值的实部皆小于 0。

发散态：其中一个特征值实部大于 0。上面的例子中，如果将 \(\mathbf{A}\) 变为 \(-\mathbf{A}\)，特征值也会变号，结果发散。

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

总结：原方程组有两个相互耦合的未知函数，\(u_1, u_2\) 相互耦合，而特征值和特征向量的作则就是解耦，也就是对角化。回到原方程组 \(\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}=\mathbf{Au}\)，将 \(u\) 表示为特征向量的线性组合\(\mathbf{u}=\mathbf{Sv}\)，代入原方程有 \(\mathbf{S}\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=\mathbf{ASv}\)，两边同乘以 \(\mathbf{S}^{-1}\)，得 \(\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=\mathbf{S}^{-1}\mathbf{ASv}=\mathbf{\Lambda v}\)。以特征向量为基，将 \(\mathbf{u}\) 表示为 \(\mathbf{Sv}\)，得到关于 \(\mathbf{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}\)，这是一个各未知函数间没有联系的方程组，它们的解的一般形式为 \(\mathbf{v}(t)=\mathrm{e}^{\boldsymbol\Lambda t}\mathbf{v}(0)\)，则原方程组的解的一般形式为 \(\mathbf{u}(t)=\mathrm{e}^{At}\mathbf{u}(0)=\mathbf{S}\mathrm{e}^{\boldsymbol\Lambda t}\mathbf{S}^{-1}\mathbf{u}(0)\)。这里引入了指数部分为矩阵的形式。

指数矩阵 \(\mathrm{e}^{\mathbf{A}t}\)

在上面的结论中，我们见到了 \(\mathrm{e}^{\mathbf{A}t}\)。这种指数部分带有矩阵的情况称为指数矩阵。

理解指数矩阵的关键在于，将指数形式展开称为幂基数形式，就像 \(\mathrm{e}^x=1+\frac{x^2}{2}+\frac{x^3}{6}+\dots\) 一样，将 \(\mathrm{e}^{\mathbf{A}t}\) 展开成幂级数的形式为：

\[ \mathrm{e}^{\mathbf{A}t}=\mathbf{I}+\mathbf{A}t+\frac{(\mathbf{A}t)^2}{2}+\frac{(\mathbf{A}t)^3}{6}+\dots+\frac{(\mathbf{A}t)^n}{n!}+\dots \]

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

回到正题，我们需要证明 \(\mathbf{S}\mathrm{e}^{\boldsymbol\Lambda t}\mathbf{S}^{-1}=\mathrm{e}^{\mathbf{A}t}\)，继续使用泰勒级数：

\[ \begin{gather} \mathrm{e}^{\mathbf{A}t}=\mathbf{I}+\mathbf{A}t+\frac{(\mathbf{A}t)^2}{2}+\frac{(\mathbf{A}t)^3}{6}+\dots+\frac{(\mathbf{A}t)^n}{n!}+\dots \\ \mathrm{e}^{\mathbf{A}t}=\mathbf{SS}^{-1}+\mathbf{S\Lambda S}^{-1}t+\frac{\mathbf{S\Lambda}^2\mathbf{S}^{-1}}{2}t^2+\frac{\mathbf{S\Lambda}^3\mathbf{S}^{-1}}{6}t^3+\dots+\frac{\mathbf{S\Lambda}^n\mathbf{S}^{-1}}{n!}t^n+\dots \\ \mathrm{e}^{\mathbf{A}t}=\mathbf{S}\left(\mathbf{I}+\boldsymbol\Lambda t+\frac{\boldsymbol\Lambda^2t^2}{2}+\frac{\boldsymbol\Lambda^3t^3}{3}+\dots+\frac{\boldsymbol\Lambda^nt^n}{n}+\dots\right)\mathbf{S}^{-1} \\ \mathrm{e}^{\mathbf{A}t}=\mathbf{S}\mathrm{e}^{\boldsymbol\Lambda t}\mathbf{S}^{-1} \end{gather} \]

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

最后，我们来看看什么是 \(\mathrm{e}^{\boldsymbol\Lambda t}\)，我们之所以将 \(\mathrm{e}^{\mathbf{A}t}\) 变为对角矩阵，是因为对角矩阵简单、没有耦合，\(\mathrm{e}^{\boldsymbol\Lambda t}=\begin{pmatrix}\mathrm{e}^{\lambda_1t}&0&\dots&0\\0&\mathrm{e}^{\lambda_2t}&\dots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\dots&\mathrm{e}^{\lambda_nt}\end{pmatrix}\)。

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

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

继续推广，对于 \(5\) 阶微分方程 \(y'''''+by''''+cy'''+dy''+ey'+f=0\)，则可以写作 \(\begin{pmatrix}y'''''\\y''''\\y'''\\y''\\y'\end{pmatrix}=\begin{pmatrix}-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{pmatrix}\begin{pmatrix}y''''\\y'''\\y''\\y'\\y\end{pmatrix}\)，这样我们就把一个五阶微分方程化为 \(5\times 5\) 一阶方程组了，然后就是求特征值和特征向量了。

Markov 矩阵、Fourier 级数

Markov 矩阵

Markov 矩阵是指具有以下两个特性的矩阵：

1. 矩阵中的所有元素 大于等于 0(因为 Markov 矩阵与概率有关，而概率是非负的)；
2. 每一列的元素之和为 1。

对于 Markov 矩阵，我们关心幂运算过程中的稳态。与上一讲不同，指数矩阵关系特征值是否为 0，而幂运算要达到稳态需要特征值为 1。

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

1. Markov 矩阵必有特征值为 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+\dots+c_n\lambda_n^kx_n\)，从这个公式很容易看出幂运算的稳态。比如我们取 \(\lambda_1=1\)，其他特征值的绝对值均小于 1，经过 \(k\) 次迭代后，其他项都趋近于 0，\(k\to\infty\) 时，有稳态 \(u_k=c_1x_1\)，这也就是初始条件 \(u_0\) 的第一个分量。

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

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

如果 \(\det(A-\lambda I)=0\)，则有 \(\det(A-\lambda I)^\top=0\)，根据矩阵转置的性质有 \(\det(A^\top-\lambda I^\top)=0\)，即 \(\det(A^\top-\lambda I)=0\)。这正是 \(A^\top\) 特征值的计算式。

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

---

接下来介绍 Markov 矩阵的应用，我们用麻省和加州这两个州的人口迁移为例：

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

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

计算特征值：我们知道 Markov 矩阵的一个特征值为 \(\lambda_1=1\)，则另一个特征值可以直接从迹算出 \(\lambda_2=0.7\)。

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

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

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

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

Fourier 级数

在介绍 Fourier 级数之前，先来回顾一下投影。

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

写为矩阵形式有 \(\Bigg[q_1\ q_2\ \dots\ q_n\Bigg]\begin{pmatrix}x_1\\x_2\\\vdots\\x_n\end{pmatrix}=v\)，即 \(Qx=v\)。所以有 \(x=Q^{-1}v\)，而对标准正交基有 \(Q^\top=Q^{-1}\)，所以我们不需要计算逆矩阵可直接得出 \(x=Q^\top v\)。此时对于 \(x\) 的每一个分量有 \(x_i=q_i^\top v\)。

接下来介绍 Fourier 级数。先写出 Fourier 级数的展开式：

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

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

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

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

在本例中，由于 Fourier 级数使用正余弦函数，它们的周期都可以取 \(2\pi\)，所以有 \(f^\top g=\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\)，其余的三角函数族正交性结果可以参考Fourier 级数的“希尔伯特空间的解读”一节。

最后我们来看 \(\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\to a_1=\frac{1}{\pi}\int_0^{2\pi}f(x)\cos x\mathrm{d}x\)。

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

复数矩阵和快速 Fourier 变换

复数矩阵

复数向量 \(\mathbf{z}=\begin{pmatrix}z_1\\z_2\\\vdots\\z_n\end{pmatrix}\) 的每一个分量都是复数。此时 \(z\) 不再属于 \(\mathbb{R}^n\) 实向量空间，它现在属于 \(\mathbb{C}^n\) 复向量空间。

计算复向量的模

实向量的模算 \(\left|\mathbf{v}\right|=\sqrt{\mathbf{v}^\top\mathbf{v}}\)。对于复向量有 \(z^\top z=\begin{pmatrix}z_1&z_2&\dots&z_n\end{pmatrix}\begin{pmatrix}z_1\\z_2\\\vdots\\z_n\end{pmatrix}=z_1^2+z_2^2+\dots+z_n^2\)，这里 \(z_i\) 是复数，平方后虚部为负，求模时本应相加的运算变成了减法。(如向量 \(\begin{pmatrix}1 & i\end{pmatrix}\)，右乘其转置后结果为 0，但此向量的长度显然不是零。)

使用 \(\left|z\right|=\sqrt{\bar{z}^\top z}=\begin{pmatrix}\bar z_1&\bar z_2&\dots&\bar z_n\end{pmatrix}\begin{pmatrix}z_1\\z_2\\\vdots\\z_n\end{pmatrix}\)，即使用向量共轭的转置乘以原向量。(如向量 \(\begin{pmatrix}1 & i\end{pmatrix}\)，右乘其共轭转置后结果为 \(\begin{pmatrix}1&-i\end{pmatrix}\begin{pmatrix}1\\i\end{pmatrix}=2\)。)

我们把共轭转置乘以原向量记为 \(z^\mathrm{H}z\)，H 是 Hermite(厄米特)。

计算向量的内积

有了复向量模的计算公式，同理可得，对于复向量，内积不再是实向量的 \(\mathbf{y}^\top\mathbf{x}\) 形式，复向量内积应为 \(\mathbf{y}^\mathrm{H}\mathbf{x}\)。

对称性

对于实矩阵，\(A^\top=A\) 即可表达矩阵的对称性。对于复矩阵，则需要一次共轭 \(\bar{A}^\top=A\)，叫做厄米特矩阵，有性质 \(A^\mathrm{H}=A\)。

正交性

对于实标准正交向量：\(q_i^\top q_j=\begin{cases}0\quad i\ne j\\1\quad i=j\end{cases}\)。对于复向量我们需要共轭：\(\bar{q}_i^\top q_j=q_i^\mathrm{H}q_j=\begin{cases}0\quad i\ne j\\1\quad i=j\end{cases}\)。

对于实标准正交矩阵有 \(Q^\top Q=I\)。对于复矩阵则有 \(Q^\mathrm{H}Q=I\)。

正交(orthogonal)在复数情况下也有了新名字：酉(unitary)。酉矩阵与正交矩阵类似，满足 \(Q^\mathrm{H}Q=I\) 的性质。 Fourier 矩阵就是酉矩阵。

Fourier 矩阵

\(n\) 阶 Fourier 矩阵 \(F_n=\begin{pmatrix}1&1&1&\dots&1\\1&w&w^2&\dots&w^{n-1}\\1&w^2&w^4&\dots&w^{2(n-1)}\\\vdots&\vdots&\vdots&\ddots&\vdots\\1&w^{n-1}&w^{2(n-1)}&\dots&w^{(n-1)^2}\end{pmatrix}\)，对于每一个元素有 \((F_n)_{ij}=w^{(i-1)(j-1)}\quad i,j=1,2,\dots,n\)。原根 \(w\) 满足 \(w^n=1\)，\(w=\mathrm{e}^{\mathrm{i}2\pi/n}=\cos\frac{2\pi}{n}+i\sin\frac{2\pi}{n}\)。

• 我们现在来看四阶 Fourier 矩阵，有 \(w=i,\ w^2=-1,\ w^3=-i,\ w^4=1\)，\(F_4=\begin{pmatrix}1&1&1&1\\1&i&i^2&i^3\\1&i^2&i^4&i^6\\1&i^3&i^6&i^9\end{pmatrix}=\begin{pmatrix}1&1&1&1\\1&i&-1&-i\\1&-1&1&-1\\1&-i&-1&i\end{pmatrix}\)。

矩阵的四个列向量正交。注意到，\(F_4\) 的列向量并不是标准的，我们可以给矩阵乘上系数 \(\frac{1}{2}\)，得到标准正交矩阵 \(F_4=\frac{1}{2}\begin{pmatrix}1&1&1&1\\1&i&-1&-i\\1&-1&1&-1\\1&-i&-1&i\end{pmatrix}\)。此时有 \(F_4^\mathrm{H}F_4=I\)，于是该矩阵的逆矩阵也就是其共轭转置 \(F_4^\mathrm{H}\)。

快速 Fourier 变换(FFT)

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

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

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

• 第一个矩阵由单位矩阵 \(I\) 和对角矩阵 \(D=\begin{pmatrix}1&&&&\\&w&&&\\&&w^2&&\\&&&\ddots&\\&&&&w^{31}\end{pmatrix}\) 组成，我们称这个矩阵为修正矩阵，显然其计算量来自 \(D\) 矩阵，对角矩阵的计算量约为 32 即这个修正矩阵的计算量约为 32，单位矩阵的计算量忽略不计。

• 第二个矩阵是两个 \(F_{32}\) 与零矩阵组成的，计算量约为 \(2\times 32^2\)。

• 第三个矩阵通常记为 \(P\) 矩阵，这是一个置换矩阵，其作用是讲前一个矩阵中的奇数列提到偶数列之前，将前一个矩阵从 \(\Bigg[x_0\ x_1\ \dots\Bigg]\) 变为 \(\Bigg[x_0\ x_2\ \dots\ x_1\ x_3\ \dots\Bigg]\)，这个置换矩阵的计算量也可以忽略不计。(这里教授似乎在黑板上写错了矩阵，可以参考FFT、How the FFT is computed做进一步讨论。)

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

来看化简后计算量，\(2\left(2\left(2\left(2\left(2\left(2(1)^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}\)。

其他

矩阵的 F 范数

设矩阵 \(\mathbf{A}=(a_{i,j})_{m\times n}\)，则其 F 范数为：\(\lVert\mathbf{A}\rVert_F=\sqrt{\sum_{i,j}a_{i,j}^2}\)。它是向量的 L2 范数的推广。

迹的性质

\(\mathbf{A}\) 的 F 范数等于 \(\mathbf{AA}^\top\) 的迹的平方根：\(\lVert\mathbf{A}\rVert_F=\sqrt{\mathrm{tr}(\mathbf{AA}^\top)}\)。

叉积性质

\[ \mathbf{u}\times(\mathbf{v}\times\mathbf{w})=(\mathbf{u}\cdot\mathbf{w})\mathbf{v}-(\mathbf{u}\cdot\mathbf{v})\mathbf{w} \]

三维向量的混合积

\[ \lvert\mathbf{u}\mathbf{v}\mathbf{w}\rvert=(\mathbf{u}\times\mathbf{v})\cdot\mathbf{w}=\mathbf{u}\cdot(\mathbf{v}\times\mathbf{w}) =\begin{vmatrix} u_x & u_y & u_z \\ v_x & v_y & v_z \\ w_x & w_y & w_z \end{vmatrix} \]