对线性方程组的研究是促进线性代数这么学科诞生的重要原因。接下来，笔者首先给出线性方程组是否有解的理论分析，然后介绍如何通过消元法和迭代法求解它。
1 线性方程组
$\forall A \in \mathbb{R}^{m \times n}$，如何求解线性方程组 $A x = b$？
1.1 理论分析
将 $A$ 写作 $[a_1,\cdots,a_n]$，则 $Ax = \sum_{i} a_i x_i$ 是以 $a_1, ..., a_n$ 为基的线性组合。若 $b \in \operatorname{range}(A)$ 则存在 $x$ 使得 $\sum_{i} a_i x_i = b$。若 $a_1, ..., a_n$ 线性无关，则 $b$ 只有一种表示方法，否则有无数种。这就对应着原方程组有唯一解（$\operatorname{rank}(A)=n$）和有无穷解（$\operatorname{rank}(A)=\operatorname{rank}(A,b)$）。因此可总结如下：
-
若
$A$为非奇异阵（此时需要$m=n$），则$x = A^{-1} b$。根据克拉默法则可知$x_i = \frac{|A_i|}{|A|}$； -
若
$A$为奇异阵且$b \in \operatorname{range} (A)$，则$Ax=b$有无穷解； -
若
$A$为奇异阵且$b \neq \operatorname{range} (A)$，则$Ax=b$无解。
其中，一个矩阵是奇异的当且仅当它至少有一个特征值为零。
可是，克拉默法则是怎样得到的？下面给出了一个有趣的几何学解释。 首先思考如下等价变换：
$$ 1 \cdot x_2 = \operatorname{det} \Bigg( \left[ \begin{matrix} 1 & x_1 \\ 0 & x_2 \end{matrix} \right] \Bigg), \quad 1 \cdot x_1 = \operatorname{det} \Bigg( \left[ \begin{matrix} x_1 & 0 \\ x_2 & 1 \end{matrix} \right] \Bigg) $$
上述两个等式的左右分别是对同一个平行四边形的有向面积的不同计算方法。扩展到 $n$ 维线性空间，即 $x_i = \operatorname{det} (X)$，其中
$$ \begin{aligned} X_{:, i} &= [x_1, ..., x_n]^\textrm{T} \in \mathbb{R}^{n} \\ X_{:, j} &= \Bigg[ \underbrace{0, ..., 1, ..., 0}_{\textrm{仅第 }i\textrm{个位置为 }1} \Bigg]^\textrm{T}, \forall j \neq i \end{aligned} $$
对于向量 $x = [x_1, ..., x_n]^\textrm{T}$，假设经过某个变换 $A$，$x$ 变成了 $b$，上面的等价变换会怎样？首先，等式左边的 $x_i$ 被缩放了 $|A|$ 倍，而等式右边的 $\operatorname{det}(X)$ 变成了 $\operatorname{det}(X')$，其中 $X'_{:, j \neq i} = A_{:, j \neq i}$，这是因为基从旧的变成了 $\{ A_{:, j} \}_{j=1, ..., n}$；而 $X'_{:, i} = [b_1, ..., b_n]^\textrm{T}$。变形即可得到克拉默法则
$$ x_i = \frac{|X'|}{|A|} $$
1.2 高斯消元法
对记计算机而言，依据克拉默法则计算 $x$ 的复杂度为 $\mathcal{O} (n^4)$，这是不可接受的。计算机还会采用消元的方法求解线性方程组。最经典的方法是 高斯消元法。
如果 $A$ 是下三角矩阵，则 $x$ 非常容易计算。因为
$$ x_i = \frac{b_i - \sum_{j=1}^{i-1} a_{ij} x_j}{a_{ii}} $$
依次求解 $x_1,x_2,\cdots,x_n$ 即可，这就是 前向算法。这个过程要求 $\forall i: a_{ii} \neq 0$，复杂度为 $\mathcal{O} (n^2)$。
如果 $A$ 是上三角矩阵，则有 后向算法：
$$ x_i = \frac{b_i - \sum_{j=i+1}^n a_{ij} x_j}{ a_{ii}} $$
若 $A$ 既不是上三角矩阵，也不是下三角矩阵，怎么办？我们可以通如一系列的初等变换将其转换为上 / 下三角矩阵。
以上三角矩阵为例，我们有无数种选择可以将 A 通过初等变换转换为上三角矩阵。高斯消去法 作为一种方法，固定了一个程式：第 $i$ 轮高斯消去是将从第 $i+1$ 行到最后一行的第 $i$ 个元素依次通过和第 $i$ 行比较消为 $0$。 通过执行高斯消去法，我们可以将 $A$ 转换成下三角矩阵和上三角矩阵的乘积，即 $A = LU$，其中 $L$ 为下三角矩阵，$U$ 为上三角矩阵，这就是 LU 分解。那么对于 $LUx=b$，先后使用前向算法和后向算法即可得到 $x$，读者可自行验证。
1.3 迭代方法
高斯消去法的复杂度为 $\mathcal{O} (n^2)$。对于大规模稀疏矩阵而言，更好的求解策略是迭代。
因为 $x_i = \frac{b_i - \sum_{j \neq i} a_{ij} x_j}{a_{ii}}$，所以有 Jacobi 迭代方法 如下：
$$ x_i^{(k+1)} \leftarrow \frac{b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k)} - \sum_{j=i+1}^n a_{ij} x_j^{(k)}}{a_{ii}} $$
如果每次用最新的 $x_j |_{j=1,\cdots,i-1}$，就得到了 Gauss-Seidel 迭代方法：
$$ x_i^{(k+1)} \leftarrow \frac{b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k+1)} - \sum_{j=i+1}^n a_{ij} x_j^{(k)}}{a_{ii}} $$
接下来，我们采用 矩阵分解 的方式表达迭代方程。令 $D=\operatorname{diag}(A)$，$E$ 是 $A$ 的严格下三角部分的负，$F$ 是 $A$ 的严格下上角部分的负，则 $A=D-E-F$。则对于 Jacobi 迭代方法有
$$ x^{(k+1)} = D^{-1} (E + F) x^{(k)} + D^{-1} b $$
对于 Gauss-Seidel 迭代方法有
$$ x^{(k+1)} = (D-E)^{-1} F x^{(k)} + (D-E)^{-1} b $$
如果左右乘上松弛变量则有 $\omega (D-E-F) = \omega x$，从而得到松弛法迭代公式
$$ x^{(k+1)}=(D-\omega E)^{-1} \big(\omega F + (1-\omega)D \big) x^{(k)} + (D- \omega E)^{-1}\omega b $$
2 正交
2.1 正交集
正交集（orthogonal set）是两两正交的向量的集合。如果一个正交集不含零向量，则这个正交集线性无关。
证明思路：$\forall i: 0 = c_1 v_1 + \cdots + c_p v_p \leftrightarrow 0 = 0 ⋅ v_i =(c_1 v_1 + \cdots + c_p v_p) ⋅ v_i$。由 $v_i ⋅ v_j = 0|_{i \neq j}$ 可知 $c_i=0$。
因此正交集可以作为向量空间 $\mathbb{R}^p$ 的一组正交基。进一步地，可定义单位正交集（orthonormal set）。
2.2 正交投影
对任意向量 $y \in \mathbb{R}^p$，如何求解 $y$ 在一组正交基 ${v_1,\cdots,v_p}$ 各个基上的坐标（系数）？设 $y = \sum_{i=1}^p c_i v_i$，则
$$ c_i = \frac{y \cdot v_i}{v_i \cdot v_i} $$
这个结论亦可用于正交投影时各个分量上坐标的计算。
2.3 正交矩阵
正交矩阵（orthogonal matrix）是拥有单位正交的列向量和行向量的方阵（因此更合理的名称应当是单位正交矩阵，orthonormal matrix）。
正交矩阵使得单位正交基在正交变换之后仍然 保持单位长度、相互垂直（刚体运动）。正交矩阵 $U$ 满足 $U^{-1}=U^\textrm{T}$。这是一个非常优越的性质。因为很多运算涉及到求解矩阵的逆，计算机对一般的矩阵求逆的复杂度很高。如果能和正交矩阵扯上关系，就会方便处理很多。实际上 QR 分解和特征值分解就是这么干的。
$U^{-1}=U^\textrm{T}$ 的证明思路：$U=[u_1,\cdots,u_n]$，$U^\textrm{T}=[u_1^\textrm{T},\cdots,u_n^\textrm{T}]^\textrm{T}$，乘一下即得证。
若 $U$ 是列向量相互单位正交的矩阵且非方阵，虽然 $U$ 上没有定义矩阵的逆，但仍然有 $U^\textrm{T} U=I$。
酉矩阵（Unitary Matrix）是正交阵在复空间上的推广。
左乘正交阵是在实施旋转变换，这是因为 $|U^\textrm{T} U| = 1$。$Ux$ 将 $x$ 旋转过去，$U^{-1} Ux = U^\textrm{T} Ux = Ix = x$ 又将 $x$ 旋转回来了。
2.4 构造正交基
如何从任意一组基 $\{x_1,\cdots,x_n \}$ 中构造出一组正交基 $\{ q_1, ..., q_n \}$？
可通过如下方式执行：
$$ \begin{aligned} q_1 &= x_1 \\ q_2 &= x_2 - \frac{x_2 \cdot q_1}{q_1 \cdot q_1} q_1 \\ q_3 &= x_3 - \frac{x_3 \cdot q_1}{q_1 \cdot q_1} q_1 - \frac{x_2 \cdot q_2}{q_2 \cdot q_2} q_2 \\ & \cdots \\ q_n &= x_n - \sum_{i=1}^{n-1} \frac{x_i \cdot q_i}{q_i \cdot q_i} q_i \end{aligned} $$
如何理解这个过程？以 $x_2$ 为例：将 $x_2$ 正交分解，其中一个分量在 $q_1$ 上，系数为 $\frac{x_2 \cdot q_1}{q_1 \cdot q_1}$，另一个分量就是 $q_2$。后续的过程均是如此。每次得到的 $q_i$ 是 $x_i$ 在现有的正交基上分解完毕之后的补。若要求构造出一组单位正交基 $\{ q_1, ..., q_n \}$，只需在上述每个步骤后面附上单位化的操作：
$$ \forall j = 1, ..., n: \quad q_j = x_j - \sum_{i=1}^{j-1} (x_i \cdot q_i) q_i, \quad q_j = \frac{q_j}{|q_j|} $$
仔细观察上式，若将 $x_j \cdot q_i$ 记为 $r_{ij}$（$1 \leq i < q \leq n$），将 $\Vert q_j \Vert$ 记为 $r_{jj}$，则 $R=[r_{ij}]_{n \times n}$ 是一个上三角矩阵。将 $[q_1,\cdots,q_n]$ 记为 $Q$，则 $X=QR$。也就是说，我们将任意矩阵 $X$ 分解成了一个正交阵和上三角矩阵的乘积，这就是 $QR$ 分解。类似于高斯消去法，这个分解也可以拿来求解线性方程组：
- 对于
$QR x = b$，令$Rx = y$，则$y = Q^\textrm{T} b$； - 对于
$Rx = y$，根据前向算法得到$x$。
3 二次型
3.1 二次函数
对于一个二次函数（方程），增加非二次项，不会改变曲线的形状，只会在 $xy$ 轴上产生一些旋转、平移的效果。因此，研究二次函数（方程）的性质，只关心二次部分就够了。只包含二次部分的函数（方程）被称为二次齐次函数（方程）。
显然，我们可以通过对称矩阵来描述二次齐次函数（方程）。例如
$$ ax^2+2bxy+cy^2 \leftrightarrow \left[ \begin{matrix} x & y \end{matrix} \right] \left[ \begin{matrix} a & b \\ b & c \end{matrix} \right] \left[ \begin{matrix} x \\ y \end{matrix} \right] \leftrightarrow x^\textrm{T} A x $$
可以将对称阵 $A$ 正交对角化，得 $A=PDP^\textrm{T}$，只保留伸缩变换的部分。在图像上，这就是将二次型对应的函数（方程）的曲线扶正。以圆锥曲线之一椭圆为例，这个操作相当于把一个斜着的椭圆摆正并平移到坐标轴的中心。
3.2 正定
对于二次函数 $f(x)=x^\textrm{T} Ax$，若 $\forall x \in R∧x \neq 0$，$f(x)>0$，则 $f(x)$ 是正定二次型，$A$ 是正定矩阵，记为 $A \succ 0$。同理可定义半正定矩阵、负定矩阵、半负定矩阵。正定阵必然可逆，但是不一定可对角化。
判断一个对称阵 $A$ 是否是正定矩阵的重要方法是对其执行正交对角化 $A=PDP^\textrm{T}$，若 $D$ 中每一个特征值都大于零，则正定。同理可判定是否为半正定矩阵、负定矩阵、半负定矩阵。为什么可以这样判断？这是因为 $x^\textrm{T} Ax$ 和 $x^\textrm{T} Dx$ 本质是一样的，后者只是移除了其中旋转的部分。将 $x^\textrm{T} Dx$ 展开就得到了 $\sum_{i=1}^b d_{ii} x_i^2$，只有当所有的 $d_{ii}$ 都恒大于零时才可以保证 $\sum_{i=1}^b d_{ii} x_i^2$ 恒大于零。也就是说，需要 $A$ 的特征值都大于零。
正定 / 负定阵一定可逆。因为对应的特征值不可能为零。
3.3 多元二次函数的极值
$f(x)=x^\textrm{T} Ax$（s.t. $|(|x|)|=1$）在何时取得极值？
当 $A$ 为对称阵时，$f(x)$ 的最大值为 $A$ 的最大特征值，$f(x)$ 的最小值为 $A$ 的最小特征值。
这个结论可以通过拉格朗日乘子法得到。
证明思路：当 $A$ 为对称阵时，对 $A$ 执行正交对角化，则 $f(x)=(P^\textrm{T} x)^\textrm{T} D(P^\textrm{T} x)$，$P^\textrm{T}$ 作为正交阵只对 $x$ 产生旋转变换，不改变其二范数，故 $f(x)=f(y)=y^\textrm{T} Dy$，展开再使用拉格朗日乘子法即可。
4 奇异值分解
对称阵可以进行正交对角化。如何将这个优美的性质运用到任意矩阵 $A$ 上？
4.1 理论分析
$\forall A \in R^{m \times n}$，因为 $m \times n = (m \times m)⋅(m \times n)⋅(n \times n)$，不违背矩阵乘法规则，所以必然存在 1 个 $m \times n$ 的矩阵 $\Sigma$，一个 $m \times m$ 的矩阵 $U$，和一个 $n \times n$ 的矩阵 $V$ 使得 $A=U \Sigma V^\textrm{T}$。那么问题来了，$U$、$\Sigma$ 和 $V$ 分别怎么算？
我们来看看 $A^\textrm{T} A$ 和 $AA^\textrm{T}$。
因为 $A^\textrm{T} A = (U \Sigma V^\textrm{T})^\textrm{T} U \Sigma V^\textrm{T} = V (\Sigma^\textrm{T} \Sigma) V^\textrm{T}$ 且 $A^\textrm{T} A$ 是正交阵，所以 $V$ 是 $A^\textrm{T} A$ 的特征向量组成的正交阵！$\Sigma^\textrm{T} \Sigma$ 是 $A^\textrm{T} A$ 的特征值组成的对角阵。同理，因为 $A A^\textrm{T} = (U \Sigma V^\textrm{T}) (U \Sigma V^\textrm{T})^\textrm{T} = U (\Sigma \Sigma^\textrm{T}) U^\textrm{T}$ 且 $A A^\textrm{T}$ 是正交阵，所以 $U$ 是 $A A^\textrm{T}$ 的特征向量组成的正交阵，$\Sigma \Sigma^\textrm{T}$ 是 $A A^\textrm{T}$ 的特征值组成的对角阵。
以上这段分析就已经告诉我们 $U$ 和 $V$ 分别怎么算了。
那么 $\Sigma$ 呢？要知道 $\Sigma^\textrm{T} \Sigma$ 是一个 $n \times n$ 的对角阵而 $\Sigma \Sigma^\textrm{T}$ 是一个 $m \times m$ 的对角阵，而 $m$ 可以不等于 $n$。尽管如此，设 $k=\min \{m,n\}$，由 $(A^\textrm{T} A)^\textrm{T}=AA^\textrm{T}$ 可知二者前 $k$ 个对角元素必然相等，不妨设为 $\lambda_1,\cdots,\lambda_k$。
首先，必然有 $\lambda_1 \geq \cdots \geq \lambda_k \geq 0$，只有这样我们才可以对这些数开根号。
证明思路：因为 $\forall x \in \mathbb{R}^n$，$x^\textrm{T} (A^\textrm{T} A) x = (Ax)^\textrm{T} (Ax) =\Vert Ax \Vert_2 \geq 0$，所以 $A^\textrm{T} A$ 是半正定矩阵。同理可证 $A A^\textrm{T}$ 也是半正定矩阵。结合半正定矩阵的所有特征值必然均为非负数即可得。
设 $\operatorname{rank}(A)=r$，则我们称 $\sigma_i=\sqrt{\lambda_i} |_{i=1,\cdots,r}$ 为 $A$ 的奇异值（显然有 $r \leq k$）。现在我们来说 $\Sigma$ 怎么算：
$$ \Sigma = \left[ \begin{matrix} (\sqrt{\Sigma^\textrm{T} \Sigma})_{r \times r} & 0 \\ 0 & 0 \end{matrix} \right] = \left[ \begin{matrix} (\sqrt{\Sigma \Sigma^\textrm{T}})_{r \times r} & 0 \\ 0 & 0 \end{matrix} \right] = \left[ \begin{matrix} \operatorname{diag} (\sigma_1, ..., \sigma_r) & 0 \\ 0 & 0 \end{matrix} \right] $$
其中 $(\sqrt{\Sigma^\textrm{T} \Sigma})_{r \times r}$ 表示对 $\Sigma^\textrm{T} \Sigma$ 每个元素依次开根号并只保留前 $r \times r$ 的部分。
有了 $\{\sigma_i \}_{i=1, ..., r}$，我们还可以换一种方式计算 $U$：
令 $U = (u_1, \cdots, u_m)$，由 $A = U \Sigma V^\textrm{T}$ 可知
$$ \forall i=1, \cdots, r: \quad u_i = \frac{1}{\sigma_i} A v_i $$
至于 $u_i |_{i=r+1, \cdots ,m}$，可用 $\{u_1,\cdots,u_m\}$ 相互单位正交得到。
4.2 SVD 的执行步骤
执行奇异值分解的步骤如下：
Step 1：计算 $\operatorname{rank} (A) = r$；
Step 2：对 $A^\textrm{T} A$ 执行正交对角化：$A^\textrm{T} A = V (\Sigma^\textrm{T} \Sigma) V^\textrm{T}$，得到 $V$ 和 $\Sigma^\textrm{T} \Sigma$；
Step 3：对 $A A^\textrm{T}$ 执行正交对角化：$A A^\textrm{T} = U (\Sigma \Sigma^\textrm{T}) U^\textrm{T}$，得到 $U$ 和 $\Sigma \Sigma^\textrm{T}$；
Step 4：计算 $\Sigma$：
$$ \Sigma = \left[ \begin{matrix} (\sqrt{\Sigma^\textrm{T} \Sigma})_{r \times r} & 0 \\ 0 & 0 \end{matrix} \right] $$
或者
$$ \Sigma = \left[ \begin{matrix} (\sqrt{\Sigma \Sigma^\textrm{T}})_{r \times r} & 0 \\ 0 & 0 \end{matrix} \right] $$
4.3 SVD 的应用
类似普分解定理，$\forall A \in R^{m \times n}$ 有
$$ A = \sum_{i=1}^r \sigma_i u_i v_i^\textrm{T} $$
一共有 $r$ 个分量，按照 $\sigma_i$ 的大小降序排列。这里可以解释为什么 $\Sigma$ 只选取 $\sqrt{\Sigma^\textrm{T} \Sigma}$ 的前 $r \times r$ 的部分。这是因为秩代表了矩阵 $A$ 线性无关的列向量的个数，从而决定了分量的个数。如果我们只保留这 $r$ 个分量的前部分分量（例如前 20% 的分量），就能够实现数据压缩。
最后
著名的数学科普博主 3Blue1Brown 在他的油管和 B 站账号上均发布了名为 线性代数的本质 系列视频，其中提供了本文所述概念的动态演示，推荐读者观看。
转载申请
本作品采用 知识共享署名 4.0 国际许可协议 进行许可，转载时请注明原文链接。您必须给出适当的署名，并标明是否对本文作了修改。