Lecture 19 - Numerical Linear Algebra

References

• Lecture: https://www.stat.cmu.edu/~ryantibs/convexopt-F18/
• Reading: 最优化: 建模算法与理论

1. Complexity of Basic Operations

在进行优化过程中, 要进行大量的线性系统的求解运算. 首先介绍相关的数值计算中的运算的基本复杂度单位: Flops (Floating Point Operations). 粗略地说, 一次 flop 是指一次浮点数运算, 例如加法, 减法, 乘法, 除法等. 我们可以将各种计算操作的复杂度用 flop 来衡量. 下面是一些常见的线性代数操作及其复杂度分析.

对于向量乘法, $\mathbf{a},\mathbf{b}\in {\mathbb{R}}^n$, 及 $c\in \mathbb{R}$, 下面是一些基本的线性代数操作及其复杂度:

• $\mathbf{a}+\mathbf{b}$ 需要 $n$ 次加法, 因此复杂度为 $\mathcal{O}(n)$ flops.
• 标量乘法 $c\mathbf{a}$ 需要 $n$ 次乘法, 因此复杂度为 $\mathcal{O}(n)$ flops.
• 内积 ${\mathbf{a}}^{\top }\mathbf{b}={\sum }_{i=1}^n{a}_i{b}_i$ 需要 $n$ 次乘法和 $n-1$ 次加法, 因此复杂度为 $\mathcal{O}(2n)$ flops.

对于矩阵-向量运算 $\mathbf{A}\in {\mathbb{R}}^{m\times n}$ 和 $\mathbf{x}\in {\mathbb{R}}^n$, 下面是一些常见的矩阵运算及其复杂度:

• 矩阵-向量乘法 $\mathbf{Ax}$ 需要 $m$ 行和 $n$ 列的乘法, 相当于 $m$ 个内积, 每个内积是长度为 $n$ 的向量, 因此复杂度为 $\mathcal{O}(2mn)$ flops.
• 特别地, 若 $\mathbf{A}$ 是 $s$-稀疏的, 即整个矩阵中只有 $s$ 个非零元素, 则矩阵-向量乘法只需要处理这 $s$ 个非零元素, 大概的复杂度为 $\mathcal{O}(2s)$ flops.
• 若 $\mathbf{A}$ 是 $k$-banded, 即矩阵集中在对角线附近宽度为 $k$ 的区域内, 则矩阵-向量乘法的复杂度大约为 $\mathcal{O}(2nk)$ flops.
• 若 $\mathbf{A}$ 是一个秩为 $r$ 的 low-rank 矩阵, 其可以拆分为 $\mathbf{A}={\sum }_{i=1}^r{\mathbf{u}}_i{\mathbf{v}}_i^{\top }$, 则乘法 $\mathbf{Ax}$ 可以通过计算 ${\sum }_{i=1}^r{\mathbf{u}}_i({\mathbf{v}}_i^{\top }\mathbf{x})$ 来实现, 相当于 $r$ 次内积和 $r$ 次标量乘法, 因此复杂度为 $\mathcal{O}(2r(n+m))$ flops.
• 若 $\mathbf{A}$ 是一个 permutation矩阵, 即每行和每列只有一个非零元素, 相当于对向量进行重排, 则没有实际的浮点数运算, 只是相当于在内存中进行数据的重排, 故在 flops 的意义上复杂度为 $\mathcal{O}(0)$ flops.

对于矩阵乘法 $\mathbf{A}\in {\mathbb{R}}^{m\times n}$ 和 $\mathbf{B}\in {\mathbb{R}}^{n\times p}$, 一般的矩阵乘法需要 $m$ 行, $p$ 列向量分别进行长度为 $n$ 的内积, 因此复杂度为 $\mathcal{O}(2mnp)$ flops. 但是对于特殊结构的矩阵, 例如 $\mathbf{A}$ 是 $s$-稀疏的, 则矩阵乘法的复杂度大约为 $\mathcal{O}(2sp)$ flops,

对于矩阵-矩阵-向量运算 $\mathbf{A}\in {\mathbb{R}}^{m\times n}$, $\mathbf{B}\in {\mathbb{R}}^{n\times p}$ 和 $\mathbf{x}\in {\mathbb{R}}^p$, 计算 $\mathbf{ABx}$ 可以通过先计算 $\mathbf{y}=\mathbf{Bx}$, 然后计算 $\mathbf{Ay}$ 来实现. 其总的复杂度为 $\mathcal{O}(2np+2mn)$ flops. 但若采取错误的计算顺序, 则先计算 $\mathbf{AB}$ 的复杂度为 $\mathcal{O}(2mnp)$ flops, 然后再计算 $(\mathbf{AB})\mathbf{x}$ 的复杂度为 $\mathcal{O}(2mp)$ flops, 总的复杂度为 $\mathcal{O}(2mnp+2mp)$ flops, 这在 $n$ 很大的情况下会非常昂贵. 因此选择正确的计算顺序对于优化计算效率非常重要.

2. Solving Linear Systems

给定一个非奇异方阵 $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ 和一个向量 $\mathbf{b}\in {\mathbb{R}}^n$, 求解线性系统 $\mathbf{Ax}=\mathbf{b}$.

• 一般而言, 若无特殊结构, 则其通过 Gaussian elimination 等类似方法进行消元的一般复杂度为 $\mathcal{O}(n^3)$ flops.

• 若 $\mathbf{A}$ 是纯对角阵, 例如 $\mathbf{A}=\text{diag}(a_1,a_2,\dots ,a_n)$, 则求解线性系统 $\mathbf{Ax}=\mathbf{b}$ 只需要进行 $n$ 次标量除法: $x_i=b_i/a_i$ for $i=1,2,\dots ,n$, 因此复杂度为 $\mathcal{O}(n)$ flops.

• 若 $\mathbf{A}$ 是一个下三角矩阵 (上三角矩阵同理), 则其总的复杂度为 $\mathcal{O}(n^2)$ flops. 其 forward substitution (上三角为 backward substitution) 求解流程如下. 给定

  $$\begin{array}{r}A=\left[\begin{array}{ccccc}{a}_{11}& 0& 0& \cdots & 0\\ a_{21}& a_{22}& 0& \cdots & 0\\ a_{31}& a_{32}& a_{33}& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ a_{n1}& a_{n2}& a_{n3}& \cdots & a_{nn}\end{array}\right],\mathbf{b}=\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\\ ⋮\\ b_n\end{array}\right],\mathbf{x}=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ ⋮\\ x_n\end{array}\right]\end{array}$$

  则

  $$\begin{array}{rl}{x}_1& =b_1/a_{11}\\ x_2& =(b_2-a_{21}{x}_1)/a_{22}\\ x_3& =(b_3-a_{31}{x}_1-a_{32}{x}_2)/a_{33}\\ & ⋮\\ x_n& =(b_n-a_{n1}{x}_1-a_{n2}{x}_2-\cdots -a_{n,n-1}{x}_{n-1})/a_{nn}\end{array}$$

  其总共的计算量为 $1+2+\cdots +(n-1)=\frac{n(n-1)}{2}$ 次乘法和 $n$ 次除法, 因此复杂度为 $\mathcal{O}(n^2)$ flops.

• 若 $\mathbf{A}$ 是一个稀疏矩阵, 则其求解线性系统的复杂度取决于稀疏矩阵的结构. 有些矩阵在进行消元之后仍然很稀疏, 例如 $k$-banded 矩阵, 其求解线性系统的复杂度大约为 $\mathcal{O}(n{k}^2)$ flops. 但有些矩阵在进行消元之后会变得非常稠密, 该现象称为 fill-in, 其求解线性系统的复杂度可能会增加到 $\mathcal{O}(n^3)$ flops. 因此, 真正影响复杂度的不是原始矩阵的稀疏性,而是在进行分解之后的稀疏性. 这也是为什么在进行矩阵分解时, 需要选择合适的 pivoting strategy 来尽量减少 fill-in 的现象.

• 若 $\mathbf{A}$ 是正交矩阵, 即 ${\mathbf{A}}^{\top }\mathbf{A}=\mathbf{I}$, 则其求解线性系统的复杂度为 $\mathcal{O}(n^2)$ flops. 这是因为对于正交矩阵 $\mathbf{A}$, 求解线性系统 $\mathbf{Ax}=\mathbf{b}$ 可以通过计算 $\mathbf{x}={\mathbf{A}}^{\top }\mathbf{b}$ 来实现, 其复杂度为 $\mathcal{O}(n^2)$ flops.

2.1 Sherman-Morrison-Woodbury Formula

这里额外补充一个求解逆矩阵的技巧. 其核心在讨论的问题是, 如果我们已经知道了 ${\mathbf{A}}^{-1}$ 或已经可以高效求解$\mathbf{Ax}=\mathbf{b}$, 那么当 $\mathbf{A}$ 发生了一个 low-rank 的 perturbation, 例如 $\mathbf{A}+\mathbf{UCV}$ 后, 我们可以有更高效的方式来求解更新后的线性系统.

Theorem (Sherman-Morrison-Woodbury Formula). 给定 $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ 是一个非奇异矩阵, $\mathbf{U}\in {\mathbb{R}}^{n\times k}$, $\mathbf{C}\in {\mathbb{R}}^{k\times k}$ 和 $\mathbf{V}\in {\mathbb{R}}^{k\times n}$ 是任意的矩阵 ($\mathbf{C}$ 需要是非奇异的), 则

$${\mathbf{A}}_{n\times n}+{\mathbf{U}}_{n\times k}{\mathbf{C}}_{k\times k}{\mathbf{V}}_{k\times n}\text{可逆}⟺{\mathbf{C}}_{k\times k}^{-1}+{\mathbf{V}}_{k\times n}{\mathbf{A}}_{n\times n}^{-1}{\mathbf{U}}_{n\times k}\text{可逆}$$

且其逆矩阵可以通过以下公式计算:

$$(\mathbf{A}+\mathbf{UCV}{)}^{-1}={\mathbf{A}}^{-1}-{\mathbf{A}}^{-1}\mathbf{U}({\mathbf{C}}^{-1}+\mathbf{V}{\mathbf{A}}^{-1}\mathbf{U}{)}^{-1}\mathbf{V}{\mathbf{A}}^{-1}.$$

特别地,

• 若 $k=1$, $C=1$, 则 $\mathbf{A}+\mathbf{u}{\mathbf{v}}^{\top }$ 可逆当且仅当 $1+{\mathbf{v}}^{\top }{\mathbf{A}}^{-1}\mathbf{u}\ne 0$, 且其逆矩阵为

  $$(\mathbf{A}+\mathbf{u}{\mathbf{v}}^{\top }{)}^{-1}={\mathbf{A}}^{-1}-\frac{{\mathbf{A}}^{-1}\mathbf{u}{\mathbf{v}}^{\top }{\mathbf{A}}^{-1}}{1+{\mathbf{v}}^{\top }{\mathbf{A}}^{-1}\mathbf{u}}.$$

• 若 $\mathbf{C}=I$, 则 $\mathbf{A}+\mathbf{UV}$ 可逆当且仅当 $\mathbf{I}+\mathbf{V}{\mathbf{A}}^{-1}\mathbf{U}$ 可逆, 且其逆矩阵为

  $$(\mathbf{A}+\mathbf{UV}{)}^{-1}={\mathbf{A}}^{-1}-{\mathbf{A}}^{-1}\mathbf{U}(\mathbf{I}+\mathbf{V}{\mathbf{A}}^{-1}\mathbf{U}{)}^{-1}\mathbf{V}{\mathbf{A}}^{-1}.$$

SMW 分解的核心在于, 对于更新后的系统 $\mathbf{A}+\mathbf{UCV}$, 其总的维度仍然是 $n\times n$, 若直接计算其逆矩阵, 则浪费了大量的关于 low-rank perturbation 的结构信息. 通过 SMW 分解, 我们可以将其逆矩阵的计算转化为对 ${\mathbf{A}}^{-1}$ 的计算 (该计算我们已经知道了), 以及对一个 $k\times k$ 的矩阵 ${\mathbf{C}}^{-1}+\mathbf{V}{\mathbf{A}}^{-1}\mathbf{U}$ 的计算 (该矩阵的维度远小于 $n$), 从而大大降低了计算的复杂度 (其中 $k$ 是 perturbation 的秩, 通常远小于 $n$).

3. Matrix Factorizations for Solving Linear Systems

矩阵分解的思路如下. 在进行求解 $\mathbf{Ax}=\mathbf{b}$ 的过程中, 我们几乎从不会直接求解 $\mathbf{A}$ 的逆矩阵 ${\mathbf{A}}^{-1}$ 来计算. 更常见的做法是先将 $\mathbf{A}$ 分解成一些特殊矩阵的乘积:

$$\mathbf{A}={\mathbf{A}}_1{\mathbf{A}}_2\cdots {\mathbf{A}}_k⟺\mathbf{Ax}={\mathbf{A}}_1{\mathbf{A}}_2\cdots {\mathbf{A}}_k\mathbf{x}=\mathbf{b}$$

则此时我们可以通过依次求解以下线性系统来得到 $\mathbf{x}$:

$$\begin{array}{rl}{\mathbf{A}}_1({\mathbf{A}}_2\cdots {\mathbf{A}}_k\mathbf{x})& :={\mathbf{A}}_1{\mathbf{y}}_1=\mathbf{b}\text{(求解 y1)}\\ {\mathbf{A}}_2({\mathbf{A}}_3\cdots {\mathbf{A}}_k\mathbf{x})& :={\mathbf{A}}_2{\mathbf{y}}_2={\mathbf{y}}_1\text{(求解 y2)}\\ & ⋮\\ {\mathbf{A}}_k\mathbf{x}& ={\mathbf{y}}_{k-1}\text{(求解 x)}\end{array}$$

其中, 每一步的求解虽然形式上仍然要处理 ${\mathbf{y}}_i={\mathbf{A}}_i^{-1}{\mathbf{y}}_{i-1}$, 但由于 ${\mathbf{A}}_i$ 的特殊结构, 其求解的复杂度通常会大大降低. 常见的特殊结构就包括前面提到的对角矩阵, 三角矩阵, 稀疏矩阵, 正交矩阵等.

该效应尤其还在当我们要求解多个线性系统 $\mathbf{A}{\mathbf{x}}_i={\mathbf{b}}_i$ for $i=1,2,\dots ,m$ 的时候更加明显. 因为在这种情况下, 我们只需要进行一次矩阵分解 $\mathbf{A}={\mathbf{A}}_1{\mathbf{A}}_2\cdots {\mathbf{A}}_k$ (该复杂度通常为 $\mathcal{O}(n^3)$ flops), 然后对于每个线性系统, 只需要依次求解 $k$ 个特殊结构的线性系统 (每个复杂度通常为 $\mathcal{O}(n^2)$ flops), 因此总的复杂度为 $\mathcal{O}(n^3+m{n}^2)$ flops, 而如果直接求解每个线性系统, 则总的复杂度为 $\mathcal{O}(m{n}^3)$ flops, 当 $m$ 很大的时候, 前者的效率会远远高于后者.

3.1 QR Decomposition

对于任意的矩阵 $\mathbf{A}\in {\mathbb{R}}^{m\times n},m\ge n$, 其都存在如下标准 QR 分解:

$${\mathbf{A}}_{m\times n}={\mathbf{Q}}_{m\times m}{\mathbf{R}}_{m\times n}$$

其中 $\mathbf{Q}$ 是一个正交矩阵, $\mathbf{R}$ 是一个上三角矩阵. 然而这里, 当 $m>n$ 的时候, $\mathbf{R}$ 的下半部分是全零的, 因此会浪费大量的存储空间. 故在实际的数值计算中, 我们通常会使用一个更为紧凑的 QR 分解形式, 称为 reduced QR decomposition, 其本质在于去除 $\mathbf{R}$ 中的全零部分, 以及 $\mathbf{Q}$ 中对应的列. 其形式如下:

$${\mathbf{A}}_{m\times n}={\mathbf{Q}}_{m\times n}{\mathbf{R}}_{n\times n}$$

其中仍有 ${\mathbf{Q}}^{\top }\mathbf{Q}=\mathbf{I}$, 但 $\mathbf{Q}$ 不再是一个方阵; $\mathbf{R}$ 是一个上三角矩阵. 有数学定理可以保证, 只要 $\mathbf{A}$ 是列满秩的, 且额外要求 $\mathbf{R}$ 的对角线元素都为正, 则 reduced QR 分解是唯一的.

通过 Gram-Schmidt 等方法, 都可以实现 QR 分解. 其总体的复杂度约为 $\mathcal{O}(m{n}^2)$ flops.

Example (QR Decomposition for OLS). 给定标签 $\mathbf{y}\in {\mathbb{R}}^n$ 和特征矩阵 $\mathbf{X}\in {\mathbb{R}}^{n\times d}$, 考虑如下最小二乘问题:

$$\underset{\mathbf{\beta }\in {\mathbb{R}}^d}{\min }\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2$$

为求解该问题, 我们可以先对 $\mathbf{X}$ 进行 reduced QR 分解, 得到

$${\mathbf{X}}_{n\times d}={\mathbf{Q}}_{n\times d}{\mathbf{R}}_{d\times d},\text{where }{\mathbf{Q}}^{\top }\mathbf{Q}={\mathbf{I}}_d.$$

由于 $\mathbf{Q}\in {\mathbb{R}}^{n\times d}$ 还不是一个完整的正交矩阵, 故额外引入一个矩阵 $\tilde{\mathbf{Q}}\in {\mathbb{R}}^{n\times (n-d)}$, 得到 $\mathbf{P}:=[\mathbf{Q}\tilde{\mathbf{Q}}]\in {\mathbb{R}}^{n\times n}$ 是一个完整的正交矩阵, 其满足 ${\mathbf{P}}^{\top }\mathbf{P}={\mathbf{I}}_n$. 此时我们有如下观察:

$$\begin{array}{rl}\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2& =\Vert {\mathbf{P}}^{\top }(\mathbf{y}-\mathbf{X}\mathbf{\beta }){\Vert }_2^2\\ & ={\Vert \left[\begin{array}{c}{\mathbf{Q}}^{\top }\\ {\tilde{\mathbf{Q}}}^{\top }\end{array}\right]\left(\mathbf{y}-\mathbf{QR}\mathbf{\beta }\right)\Vert }_2^2\text{(Since X=QR and P=[Q Q]) }\\ & ={\Vert {\mathbf{Q}}^{\top }\mathbf{y}-{\mathbf{Q}}^{\top }\mathbf{QR}\mathbf{\beta }\Vert }_2^2+{\Vert {\tilde{\mathbf{Q}}}^{\top }\mathbf{y}-{\tilde{\mathbf{Q}}}^{\top }\mathbf{QR}\mathbf{\beta }\Vert }_2^2\text{(Since ∥[a;b]∥22=∥a∥22+∥b∥22) }\\ & ={\Vert {\mathbf{Q}}^{\top }\mathbf{y}-\mathbf{R}\mathbf{\beta }\Vert }_2^2+{\Vert {\tilde{\mathbf{Q}}}^{\top }\mathbf{y}\Vert }_2^2\text{(Since Q⊤Q=I and Q⊤Q=0) }\end{array}$$

故最小化 $\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2$ 等价于最小化 ${\Vert {\mathbf{Q}}^{\top }\mathbf{y}-\mathbf{R}\mathbf{\beta }\Vert }_2^2$ (因为 ${\Vert {\tilde{\mathbf{Q}}}^{\top }\mathbf{y}\Vert }_2^2$ 与 $\mathbf{\beta }$ 无关). 而由于 $\mathbf{R}$ 是一个上三角矩阵其一定是可逆的, 因此 ${\Vert {\mathbf{Q}}^{\top }\mathbf{y}-\mathbf{R}\mathbf{\beta }\Vert }_2^2$ 的最小值为 $0$, 其最优解为 $\hat{\mathbf{\beta }}={\mathbf{R}}^{-1}{\mathbf{Q}}^{\top }\mathbf{y}$. 不过这里指出, 由于 $\mathbf{R}$ 是一个上三角矩阵, 其求解 $\hat{\mathbf{\beta }}$ 的过程可以通过 backward substitution 来实现, 因此不需要显式地计算 ${\mathbf{R}}^{-1}$ 来得到 $\hat{\mathbf{\beta }}$, 从而大大降低了计算的复杂度.

其总体复杂度如下:

• 首先对 $\mathbf{X}$ 进行 reduced QR 分解的复杂度为 $\mathcal{O}(2n{d}^2-d^3/3)$ flops.
• 计算 ${\mathbf{Q}}^{\top }\mathbf{y}$ 的复杂度为 $\mathcal{O}(2nd)$ flops.
• 通过 backward substitution 求解 $\hat{\mathbf{\beta }}$ 的复杂度为 $\mathcal{O}(d^2)$ flops.

因此总的复杂度约为 $\mathcal{O}(2n{d}^2-d^3/3)$ flops. 该复杂度在 $n\gg d$ 的情况下, 主要由对 $\mathbf{X}$ 进行 reduced QR 分解的复杂度主导, 大约为 $\mathcal{O}(2n{d}^2)$ flops.

3.2 LU Decomposition & Cholesky Decomposition

3.2.1 LU Decomposition for General Non-Singular Matrices

• 对于一般的非奇异矩阵 $\mathbf{A}\in {\mathbb{R}}^{n\times n}$, 一个通用的 LU 分解为:

  $$\mathbf{PA}=\mathbf{LU}$$

  其中 $\mathbf{P}\in {\mathbb{R}}^{n\times n}$ 是一个 permutation 矩阵, 用来进行行交换. $\mathbf{L}\in {\mathbb{R}}^{n\times n}$ 是一个下三角矩阵, 其对角线元素为 $1$, $\mathbf{U}\in {\mathbb{R}}^{n\times n}$ 是一个上三角矩阵. $\mathbf{P}$ 的引入是为了保证在进行 Gaussian elimination 的过程中, 可以通过行交换来避免出现零 pivot 的情况, 从而保证分解的稳定性.

• 在进行 LU 分解之后, 求解线性系统 $\mathbf{Ax}=\mathbf{b}$ 的流程如下:

  $$\begin{array}{rl}\mathbf{Ax}=\mathbf{b}& ⟺\mathbf{PAx}=\mathbf{Pb}⟺\mathbf{LUx}=\mathbf{Pb}\end{array}$$

  故求解可分两步进行:

  $$\begin{array}{rl}\mathbf{Ly}& =\mathbf{Pb}\text{(forward substitution 求解 y)}\\ \mathbf{Ux}& =\mathbf{y}\text{(backward substitution 求解 x)}\end{array}$$

  对于分解后的系统, 求解的复杂度约为 $\mathcal{O}(3n^2)$ flops.

• LU 分解本身的一般复杂度为 $\mathcal{O}(\frac{2n^3}{3})$ flops.

3.2.2 Cholesky Decomposition for Symmetric Positive Definite Matrices

• 特别地, 对于 $\mathbf{A}$ 是一个对称正定矩阵的情况 ($\mathbf{A}={\mathbf{A}}^{\top }\succ 0$), 存在一个更为高效的分解方法, 称为 Cholesky 分解, 其形式如下:

  $$\mathbf{A}=\mathbf{L}{\mathbf{L}}^{\top }$$

  其中 $\mathbf{L}$ 是一个下三角矩阵, 其对角线元素为正. 根据待定系数法, 对于一个一般的 $n\times n$ 的矩阵 $\mathbf{A}$, 可以推得其递推公式为:

  $$L_{ij}=\{\begin{array}{l}\left(A_{ij}-{\sum }_{k=1}^{j-1}{L}_{ik}{L}_{jk}\right)/L_{jj},\text{for }i>j\\ \sqrt{A_{jj}-{\sum }_{k=1}^{j-1}{L}_{jk}^2},\text{for }i=j\end{array}$$

  其中其对角线元素需要进行开方运算.

• 或等价地, 通过引入一个对角矩阵 $\mathbf{D}$, 我们可以得到一个不用进行开方运算的 Cholesky 分解形式, 其形式如下:

  $$\mathbf{A}=\tilde{\mathbf{L}}\mathbf{D}{\tilde{\mathbf{L}}}^{\top }$$

  其中 $\tilde{\mathbf{L}}$ 是一个下三角矩阵, 其对角线元素为 $1$, $\mathbf{D}$ 是一个对角矩阵, 其对角线元素为正.

• 由于对称性的存在, 我们大概只需要处理一半的信息, 因此 Cholesky 分解的复杂度大约为 $\mathcal{O}(\frac{n^3}{3})$ flops, 大约是 LU 分解的一半.

• 在进行 Cholesky 分解之后, 求解线性系统 $\mathbf{Ax}=\mathbf{b}$ 的流程如下:

  $$\begin{array}{rl}\mathbf{Ax}=\mathbf{b}& ⟺\mathbf{L}{\mathbf{L}}^{\top }\mathbf{x}=\mathbf{b}\\ & ⟺\mathbf{Ly}=\mathbf{b}\text{(forward substitution 求解 y)}\\ & ⟺{\mathbf{L}}^{\top }\mathbf{x}=\mathbf{y}\text{(backward substitution 求解 x)}\end{array}$$

  对于分解之后的线性系统, 求解的复杂度约为 $\mathcal{O}(2n^2)$ flops.

Example (Cholesky Decomposition for OLS). 给定标签 $\mathbf{y}\in {\mathbb{R}}^n$ 和特征矩阵 $\mathbf{X}\in {\mathbb{R}}^{n\times d}$, 考虑如下最小二乘问题:

$$\underset{\mathbf{\beta }\in {\mathbb{R}}^d}{\min }\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2$$

假设 $\mathbf{X}$ 是列满秩的, 则其闭式解为:

$$\hat{\mathbf{\beta }}=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{y}$$

其计算复杂度为:

• 首先计算 ${\mathbf{X}}^{\top }\mathbf{X}$ 的复杂度为 $\mathcal{O}(n{d}^2)$ flops.
• 同时计算 ${\mathbf{X}}^{\top }\mathbf{y}$ 的复杂度为 $\mathcal{O}(2nd)$ flops.
• 对 ${\mathbf{X}}^{\top }\mathbf{X}$ 进行 Cholesky 分解的复杂度为 $\mathcal{O}(d^3/3)$ flops.
• 通过分解后的系统求解线性系统 $({\mathbf{X}}^{\top }\mathbf{X})\mathbf{\beta }={\mathbf{X}}^{\top }\mathbf{y}$, 其复杂度为 $\mathcal{O}(2d^2)$ flops.

因此总的复杂度为 $\mathcal{O}(n{d}^2+d^3/3)$ flops. 该复杂度在 $n\gg d$ 的情况下, 主要由计算 ${\mathbf{X}}^{\top }\mathbf{X}$ 的复杂度主导, 大约为 $\mathcal{O}(n{d}^2)$ flops.

3.2.3 Cholesky-like Decomposition for Symmetric Indefinite Matrices

• 对于对称但非定的矩阵 $\mathbf{A}={\mathbf{A}}^{\top }$, 其分解的形式如下: $$\mathbf{PA}{\mathbf{P}}^{\top }=\mathbf{LD}{\mathbf{L}}^{\top }$$ 其中 $\mathbf{P}$ 是一个 permutation 矩阵, $\mathbf{L}$ 是一个下三角矩阵, 其对角线元素为 $1$, $\mathbf{D}$ 是一个对角矩阵或分块对角矩阵, 其对角线元素可以是正的, 也可以是负的.

4. Matrix Factorizations for Numerical Stability

除了进行求解线性系统之外, 线性代数中另外一大类关心的问题是数值稳定性. 在正式展开之前, 回顾一下相关的代数知识.

4.1 Eigenvalue and Eigenvector

Definition (Eigenvalue and Eigenvector). 对于方阵 $\mathbf{A}\in {\mathbb{R}}^{n\times n}$, 若存在一个非零向量 $\mathbf{v}\in {\mathbb{R}}^n$ 和一个标量 $\lambda \in \mathbb{C}$, 使得

$$\mathbf{Av}=\lambda \mathbf{v}$$

则称 $\lambda$ 是 $\mathbf{A}$ 的一个特征值 (eigenvalue), $\mathbf{v}$ 是 $\mathbf{A}$ 的一个特征向量 (eigenvector).

• $\lambda$ 是 $\mathbf{A}$ 的一个特征值当且仅当 $\mathbf{A}-\lambda \mathbf{I}$ 是一个奇异矩阵, 即 $\det (\mathbf{A}-\lambda \mathbf{I})=0$.
• $p_A(\lambda ):=\det (\mathbf{A}-\lambda \mathbf{I})$ 被称为 $\mathbf{A}$ 的特征多项式 (characteristic polynomial), 其根就是 $\mathbf{A}$ 的特征值. 对于一个 $n$ 阶矩阵 $\mathbf{A}$ 其定有 $n$ 个特征值 (重根算作多个).
• 对于一个对角阵, 其特征值就是其对角线元素, 其特征向量就是标准基向量.
• 对于矩阵 $\mathbf{A}$, 其特征值, 行列式和 trace 之间的关系如下: $$\begin{array}{rl}\det (\mathbf{A})& =\prod _{i=1}^n{\lambda }_i\\ \mathrm{tr}(\mathbf{A})& =\sum _{i=1}^n{\lambda }_i\end{array}$$ 其中特征值按照重数进行计数.

Definition (Similarity). 对于两个矩阵 $\mathbf{A},\mathbf{B}\in {\mathbb{R}}^{n\times n}$, 若存在一个可逆矩阵 $\mathbf{X}\in {\mathbb{R}}^{n\times n}$, 使得

$$\mathbf{B}=\mathbf{XA}{\mathbf{X}}^{-1}$$

则称 $\mathbf{A}$ 和 $\mathbf{B}$ 是相似的 (similar).

• 矩阵的相似关系是一种等价关系. 其表示同一个线性变换在不同基下的矩阵表示.

• 相似的矩阵具有相同的特征值, 但不一定具有相同的特征向量.

  • $\det (\lambda \mathbf{I}-\mathbf{B})=\det (\lambda \mathbf{I}-\mathbf{XA}{\mathbf{X}}^{-1})=\det (\mathbf{X}(\lambda \mathbf{I}-\mathbf{A}){\mathbf{X}}^{-1})=\det (\lambda \mathbf{I}-\mathbf{A})$, 故 $\mathbf{A}$ 和 $\mathbf{B}$ 的特征多项式相同, 从而具有相同的特征值. 但由于 $\mathbf{B}$ 的特征向量 $\mathbf{v}$ 满足 $\mathbf{Bv}=\lambda \mathbf{v}$, 则 $\mathbf{A}$ 的特征向量为 ${\mathbf{X}}^{-1}\mathbf{v}$, 从而不一定与 $\mathbf{B}$ 的特征向量相同.

• 如果存在非奇异矩阵 $\mathbf{X}$ 使得 $\mathbf{XA}{\mathbf{X}}^{-1}=:\mathbf{D}=\text{diag}({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n)$ 是一个对角矩阵, 则称 $\mathbf{A}$ 是可对角化的 (diagonalizable), 其特征值为 ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n$, 其特征向量为 ${\mathbf{X}}^{-1}{e}_i$ for $i=1,2,\dots ,n$ (其中 $e_i$ 是标准基向量).

  • 一个矩阵 $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ 是可对角化的当且仅当 $\mathbf{A}$ 的全部特征向量 ${\mathbf{v}}_1,{\mathbf{v}}_2,\dots ,{\mathbf{v}}_n$ 是线性无关的.

4.2 Eigenvalue Decomposition & Singular Value Decomposition

事实上, 上述对角化的过程就是我们在数值线性代数中非常重要的一个工具, 称为 eigenvalue decomposition. 其形式如下.

Theorem (Eigenvalue Decomposition). 对于一个可对角化的矩阵 $\mathbf{A}\in {\mathbb{R}}^{n\times n}$, 其存在 $n$ 个线性无关的特征向量 ${\mathbf{v}}_1,{\mathbf{v}}_2,\dots ,{\mathbf{v}}_n$ 和对应的特征值 ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n$, 使得

$$\mathbf{A}=\mathbf{P}\mathbf{\Lambda }{\mathbf{P}}^{-1}=\left[\begin{array}{cccc}|& |& & |\\ {\mathbf{v}}_1& {\mathbf{v}}_2& \cdots & {\mathbf{v}}_n\\ |& |& & |\end{array}\right]\left[\begin{array}{cccc}{\lambda }_1& 0& \cdots & 0\\ 0& {\lambda }_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\lambda }_n\end{array}\right]{\left[\begin{array}{cccc}|& |& & |\\ {\mathbf{v}}_1& {\mathbf{v}}_2& \cdots & {\mathbf{v}}_n\\ |& |& & |\end{array}\right]}^{-1}$$

• 注意, 对于任意的可对角化矩阵, 其分解形式为 $$\mathbf{A}=\mathbf{P}\mathbf{\Lambda }{\mathbf{P}}^{-1}$$ 若对于实对称矩阵, 其一定是可以被对角化的, 并且更进一步是可被 orthogonally diagonalizable 的, 其分解形式为 $$\mathbf{A}=\mathbf{Q}\mathbf{\Lambda }{\mathbf{Q}}^{\top }$$
• 反过来, 如果对于一个方阵 $\mathbf{M}\in {\mathbb{R}}^{n\times n}$, 我们已经明确的知道了 $\mathbf{M}=\mathbf{Q}\mathbf{\Lambda }{\mathbf{Q}}^{\top }$ 的分解形式, 其中 $\mathbf{Q}\in {\mathbb{R}}^{n\times n}$ 满足 ${\mathbf{Q}}^{\top }\mathbf{Q}=\mathbf{I}$, $\mathbf{\Lambda }=\mathrm{diag}({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n)$ 是一个对角矩阵, 则我们可以直接得出 $\mathbf{M}$ 的特征值为 ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n$, 其特征向量为 $\mathbf{Q}$ 的列向量.

然而, 并不是任意矩阵都能进行 eigenvalue decomposition. 其要求矩阵必须是方阵, 可对角化的, 且其特征值必须是实数. 一个更为一般的分解方法是 singular value decomposition, 其形式如下.

Theorem (Singular Value Decomposition). 对于任意的矩阵 $\mathbf{A}\in {\mathbb{R}}^{m\times n}$, 其存在一个 $m\times m$ 的正交矩阵 $\mathbf{U}$, 一个 $n\times n$ 的正交矩阵 $\mathbf{V}$, 以及一个 $m\times n$ 的对角型矩阵 $\mathbf{\Sigma }$, 使得

$${\mathbf{A}}_{m\times n}={\mathbf{U}}_{m\times m}{\mathbf{\Sigma }}_{m\times n}{\mathbf{V}}_{n\times n}^{\top }$$

• 其中 $\mathbf{A}$ 若为方阵, 则 $\mathbf{\Sigma }$ 直接为对角矩阵 $\mathrm{diag}({\sigma }_1,{\sigma }_2,\dots ,{\sigma }_n)$. 若非方阵, 则相当于一个对角矩阵和一个纯零矩阵的拼接. 例如

  $$\mathbf{\Sigma }={\left[\begin{array}{cccc}{\sigma }_1& 0& \cdots & 0\\ 0& {\sigma }_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\sigma }_r\\ 0& 0& \cdots & 0\\ ⋮& ⋮& & ⋮\\ 0& 0& \cdots & 0\end{array}\right]}_{m\times n}\text{or}\mathbf{\Sigma }={\left[\begin{array}{ccccccc}{\sigma }_1& 0& \cdots & 0& 0& \cdots & 0\\ 0& {\sigma }_2& \cdots & 0& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& & ⋮\\ 0& 0& \cdots & {\sigma }_r& 0& \cdots & 0\end{array}\right]}_{m\times n}$$

• $\mathbf{\Sigma }$ 中的 ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_r>0$ 被称为 $\mathbf{A}$ 的奇异值 (singular value). 事实上, $\mathbf{A}$ 的奇异值 ${\sigma }_i$ 的平方即为 ${\mathbf{A}}^{\top }\mathbf{A}$ 的特征值, 其关系为

  $${\sigma }_i(\mathbf{A})=\sqrt{{\lambda }_i({\mathbf{A}}^{\top }\mathbf{A})}$$

  Proof of the singular value relation

  因为

  $$\mathbf{A}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{\top }$$

  则

  $${\mathbf{A}}^{\top }\mathbf{A}=\mathbf{V}{\mathbf{\Sigma }}^{\top }{\mathbf{U}}^{\top }\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{\top }=\mathbf{V}{\mathbf{\Sigma }}^{\top }\mathbf{\Sigma }{\mathbf{V}}^{\top }$$

  因此 ${\mathbf{A}}^{\top }\mathbf{A}$ 的特征值为 ${\sigma }_i^2$, 从而 ${\sigma }_i=\sqrt{{\lambda }_i({\mathbf{A}}^{\top }\mathbf{A})}$.

  $\square$

  • 非零奇异值的个数事实上就相当于衡量了输入空间在经过线性变换 $\mathbf{A}$ 之后的输出空间的维度. 例如, 如果某一个对应的维度的奇异值为 $0$, 则相当于在该维度的信息会被压缩成一个零向量, 从而导致输出空间的维度降低.
  • 总的而言, 奇异值是分析一个矩阵的奇异性和稳定性的一个非常重要的工具. 其相比于 rank, 能够提供更为细致的关于矩阵的结构信息. 例如, 对于 $\mathbf{A},\mathbf{B}\in {\mathbb{R}}^{n\times n}$, 而 $\mathrm{rank}(\mathbf{B})=r<n$ 是一个不满秩的矩阵而 $\mathrm{rank}(\mathbf{A})=n$ 是一个满秩的矩阵. 则若用 operator norm 来衡量 $\mathbf{A}$ 和 $\mathbf{B}$ 之间的距离, 并且记 ${\mathcal{R}}_k$ 为所有秩不超过 $k$ 的矩阵的集合, 则我们有如下观察: $$\mathrm{dist}(\mathbf{A},{\mathcal{R}}_r)=\underset{\mathbf{B}\in {\mathcal{R}}_r}{\min }\Vert \mathbf{A}-\mathbf{B}{\Vert }_2={\sigma }_{r+1}(\mathbf{A})$$
    • 该观察说明了, 对于一个满秩的矩阵 $\mathbf{A}$ 来说, 其距离所有不满秩的矩阵的距离就是其第 $r+1$ 大的奇异值. 因此, 若 ${\sigma }_{r+1}(\mathbf{A})$ 非常小, 则说明 $\mathbf{A}$ 非常接近于一个秩为 $r$ 的矩阵, 从而在数值计算中可能会出现一些不稳定的现象.

• 其对应的几何直观如下. SVD 衡量了任何一个线性变换 $\mathbf{A}$ 对于给定向量 $\mathbf{x}$ 的作用.

  $$\mathbf{Ax}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{\top }\mathbf{x}$$
  1. 旋转 / 建系 (${\mathbf{V}}^{\top }\mathbf{x}$): Recall, 给定一个正交矩阵 ${\mathbf{V}}^{\top }$, 其左乘作用在向量 $\mathbf{x}$ 上的结果就相当于求解向量 $\mathbf{x}$ 在 $\mathbf{V}$ 的空间中的坐标表示. 因此若记 $\mathbf{c}={\mathbf{V}}^{\top }\mathbf{x}$, 则 $\mathbf{c}$ 就是向量 $\mathbf{x}$ 在 $\mathbf{V}$ 的空间中的坐标表示.
  2. 伸缩 / 缩放 ($\mathbf{\Sigma }\mathbf{c}$): 由于 $\mathbf{\Sigma }$ 是一个对角型矩阵, 其左乘作用在向量 $\mathbf{c}$ 上的结果就相当于对 $\mathbf{c}$ 的每个坐标进行伸缩. 其中 ${\sigma }_i$ 就是第 $i$ 个坐标的伸缩因子.
  3. 旋转 / 摆放 ($\mathbf{U}\mathbf{\Sigma }\mathbf{c}$): 由于 $\mathbf{U}$ 是一个正交矩阵, 其左乘作用在向量 $\mathbf{\Sigma }\mathbf{c}$ 上的结果就相当于对刚刚旋转伸缩后的坐标再翻译回原空间中的坐标表示. 特别地, 若 $\mathbf{U}$ 和 $\mathbf{V}$ 的行列式不同号, 则会发生一个额外的翻转 (flip) 的现象.

• 奇异值分解还是和矩阵压缩密切相关的. 例如, 对于一个矩阵 $\mathbf{A}$, 其 SVD 分解为 $\mathbf{A}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{\top }$, 其可以写成:

  $$\mathbf{A}=\sum _{i=1}^r{\sigma }_i{\mathbf{u}}_i{\mathbf{v}}_i^{\top }$$
  • 因此, 若只保留前 $k$ 大的奇异值, 则可以得到一个秩为 $k$ 的矩阵 ${\mathbf{A}}_k={\sum }_{i=1}^k{\sigma }_i{\mathbf{u}}_i{\mathbf{v}}_i^{\top }$, 其是 $\mathbf{A}$ 的一个近似. 该近似在 Frobenius 范数意义下是最优的, 即 ${\mathbf{A}}_k=\arg {\min }_{\mathrm{rank}(\mathbf{M})\le k}\Vert \mathbf{A}-\mathbf{M}{\Vert }_F$. 即为 low-rank approximation 的一个重要工具.

Definition (Condition Number). 对于一个矩阵 $\mathbf{A}$, 其 condition number 定义为:

$$\kappa (\mathbf{A})=\frac{{\sigma }_{\max }(\mathbf{A})}{{\sigma }_{\min }(\mathbf{A})}$$

其中 ${\sigma }_{\max }(\mathbf{A})$ 和 ${\sigma }_{\min }(\mathbf{A})$ 分别是 $\mathbf{A}$ 的最大和最小(非零)奇异值.

• 若 $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ 可逆, 则可进一步得到 $\kappa (\mathbf{A})=\Vert \mathbf{A}{\Vert }_2\Vert {\mathbf{A}}^{-1}{\Vert }_2$.
• 若 $\mathbf{A}\in {\mathbb{R}}^{m\times n}$ 是一个矩形矩阵但满秩 (例如 OLS 中的标准情况), 则类似可以得到 $\kappa (\mathbf{A})=\Vert \mathbf{A}{\Vert }_2\Vert {\mathbf{A}}^{†}{\Vert }_2$, 其中 ${\mathbf{A}}^{†}$ 是 $\mathbf{A}$ 的 Moore-Penrose pseudo-inverse.

Example (Sensitivity of Linear Systems). 给定一个线性系统 $\mathbf{Ax}=\mathbf{b}$, 记其此时的解为 ${\mathbf{x}}_0$. 现在, 假设我们对 $\mathbf{A}$ 和 $\mathbf{b}$ 都引入了一个小的扰动, 得到新的线性系统 $(\mathbf{A}+\epsilon \Delta )\mathbf{x}=\mathbf{b}+\epsilon \delta$, 其中 $\epsilon$ 是一个非常小的标量, $\Delta$ 和 $\delta$ 分别是 $\mathbf{A}$ 和 $\mathbf{b}$ 的扰动矩阵. 记新的线性系统的解为 ${\mathbf{x}}_{\epsilon }$. 则有如下关系:

$$\frac{\Vert {\mathbf{x}}_{\epsilon }-{\mathbf{x}}_0{\Vert }_2}{\Vert {\mathbf{x}}_0{\Vert }_2}\le \kappa (\mathbf{A})|\epsilon |\left(\frac{\Vert \Delta {\Vert }_2}{\Vert \mathbf{A}{\Vert }_2}+\frac{\Vert \delta {\Vert }_2}{\Vert \mathbf{b}{\Vert }_2}\right)+o(\epsilon ).$$

Example (Sensitivity Analysis of OLS by Cholesky and QR Decompositions). 考虑同样的最小二乘问题. 这里想要说明, 虽然 Cholesky 分解本身的计算复杂度更快, 但其数值稳定性却不如 QR 分解. 其原因在于, Cholesky 分解需要对 ${\mathbf{X}}^{\top }\mathbf{X}$ 进行分解, 而

$$\kappa ({\mathbf{X}}^{\top }\mathbf{X})={\kappa }^2(\mathbf{X})$$

因此若原先的 $\mathbf{X}$ 已经是一个比较 ill-conditioned 的矩阵了, 则 ${\mathbf{X}}^{\top }\mathbf{X}$ 会进一步加剧这种 ill-conditioning 的现象, 从而导致 Cholesky 分解的数值稳定性较差. 反过来, QR 分解直接对 $\mathbf{X}$ 进行分解, 因此其数值稳定性更好. 可以推导证明, 对于 QR 分解, 其最终的误差敏感性大约为 $\kappa (\mathbf{X})+\Vert \mathbf{Y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2\cdot {\kappa }^2(\mathbf{X})$, 而后者在 OLS 能够被很好地拟合的情况下将取值较小, 此时的误差敏感性会明显优于 Cholesky 分解的 ${\kappa }^2(\mathbf{X})$.

5. Indirect Methods for Solving Linear Systems

前面介绍的诸如 QR, Cholesky, LU 等分解方法, 都是直接方法 (direct method) 的范畴. 其核心思想仍然是通过直接给出一个 closed-form 的解来求解线性系统. 但在一些特定的情况下, 例如当 $\mathbf{A}$ 是一个非常大的稀疏矩阵时, 矩阵分解的方法也会非常昂贵. 此时, 我们可以考虑一些迭代方法 (iterative method) 来求解线性系统. 其核心思想是通过不断地迭代来逐渐逼近线性系统的解.

5.1 Jacobi and Gauss-Seidel Methods

考虑正定矩阵 $\mathbf{A}\in {\mathbb{S}}_{++}^n$, Jacobi method 要求解

$$\mathbf{Ax}=\mathbf{b}.$$

其推导的迭代公式如下. 在第 $k$ 轮迭代中, ${\mathbf{x}}^{(k)}$ 的第 $i$ 个坐标的更新公式为:

$$x_i^{(k)}=\frac{1}{A_{ii}}\left(b_i-\sum _{j\ne i}{A}_{ij}{x}_j^{(k-1)}\right)$$

• 该方法的推导过程如下.

  • 考虑线性系统 $\mathbf{Ax}=\mathbf{b}$ 的第 $i$ 行:

    $$A_{i1}{x}_1+A_{i2}{x}_2+\cdots +A_{ii}{x}_i+\cdots +A_{in}{x}_n=b_i$$

    将 $x_i$ 提取则有

    $$A_{ii}{x}_i=b_i-\sum _{j\ne i}{A}_{ij}{x}_j$$

  • 因此可以用上式来更新 $x_i$. 由于在第 $k$ 轮迭代中, $x_j$ 的值还没有被更新, 因此只能使用上一轮迭代的值 $x_j^{(k-1)}$ 来进行更新.

• 然而, Jacobi method 并不保证收敛. 一个判断方法如下. 可以将 $\mathbf{A}$ 分解为 $\mathbf{A}=\mathbf{D}+\mathbf{R}$, 其中 $\mathbf{D}$ 是 $\mathbf{A}$ 的对角线部分, $\mathbf{R}$ 是 $\mathbf{A}$ 的非对角线部分. 定义 $T_J=-{\mathbf{D}}^{-1}\mathbf{R}$, 并计算其谱半径 $\rho (T_J)$ (谱半径为矩阵的特征值的绝对值的最大值). 则 Jacobi method 当且仅当 $\rho (T_J)<1$ 时收敛. 例如, 当 $\mathbf{A}$ 是一个严格对角占优矩阵时, 则 Jacobi method 一定收敛.

5.2 Gauss-Seidel method

• Gauss-Seidel method 的迭代公式如下. 在第 $k$ 轮迭代中, ${\mathbf{x}}^{(k)}$ 的第 $i$ 个坐标的更新公式为:

  $$x_i^{(k)}=\frac{1}{A_{ii}}\left(b_i-\sum _{j<i}{A}_{ij}{x}_j^{(k)}-\sum _{j>i}{A}_{ij}{x}_j^{(k-1)}\right)$$

  其中对于 $j<i$ 的部分, 使用了当前轮迭代已经更新的值 $x_j^{(k)}$, 而对于 $j>i$ 的部分, 则只能使用上一轮迭代的值 $x_j^{(k-1)}$ 来进行更新.

• Gauss-Seidel method 的收敛速度通常较快, 但由于具有强序列性, 因此不太适合并行计算.