数值分析期末复习

LU Decomposition

Basics

标准 LU 分解就是不带行交换的高斯消元的矩阵语言。

Full permutation： $P_{σ}$ 是一个行向量依次为 ${\vec{e_{σ (i)}}}^{T}$ 的矩阵。
Full row scalling： $S_{s} = diag (s_{1}, s_{2}, \dots, s_{n})$
Elimination： $E = I_{n} + c \vec{e_{i}} {\vec{e_{j}}}^{T}$ 将 $c$ 倍第 $j$ 列加到第 $i$ 列。

可逆矩阵 $A \in R^{n \times n}$ 存在标准 LU 分解的充要条件是，前 $n$ 个顺序主子式均不为零
不可逆矩阵

Example: Regression of multiple experiments

Any linear system:

vector
polynomial
general

相对于已经存在的采样（线性）系统而言，新增一条采样可以是冗余采样、矛盾采样、有效采样。

采样系统有三种可叠加信息：矛盾采样条数（是否一致）、冗余采样条数（是否冗余）、有效采样条数（是否完备）。
当且仅当一个采样系统一致且完备时，可以重构出唯一的原信息系统

$n$ 次采样（ $A$ 为方阵）:

满秩时，对应无冗余的一致完备采样系统：存在唯一解

多余 $n$ 次采样（ $A$ 为 Thin 矩阵）：

列满秩时：取其中 $n$ 条采样构成一个无矛盾无冗余的完备采样系统，考虑新增剩下的 $m - n$ 条采样，它们不可能是有效采样。
- 全是冗余采样时有唯一解。
- 存在矛盾采样时无解，我们退而求其次寻找最优近似解 $‖ A \vec{x} - \vec{b} ‖$ ，即超定最小二乘问题，这等价于求解恰定问题 $A^{T} A \vec{x} = A^{T} \vec{b}$

‖ A \vec{x} - \vec{b} ‖_{2}^{2} = ‖ A \vec{x} ‖_{2}^{2} - 2 (A^{T} \vec{b}) \cdot \vec{x} + ‖ \vec{b} ‖^{2}

\nabla_{\vec{x}} ‖ A \vec{x} - \vec{b} ‖_{2}^{2} = 2 A^{T} A \vec{x} - 2 A^{T} \vec{b}

\begin{matrix} (Normal Equation) & A^{T} A \vec{x} = A^{T} \vec{b} \end{matrix}

少于 $n$ 次采样（ $A$ 为 Fat 矩阵）：

对应不完备采样系统
- 存在矛盾采样时无解，同样可以求最优近似解
- 全是冗余采样时有无数解，此时可以限定求解最小范数解

方程角度：

一个方程对应一个限制，一共 $m$ 个限制：独立的方程对解空间施加一个维度的限制（解空间 $R^{n}$ 中的一个超平面），冗余的方程不影响解的数量和自由度，矛盾的方程导致无解。
超平面的平移不影响解的数量和自由度，所以可以不看 $\vec{b}$ 只看 $A$

线性变换角度：

超定方程无解的情况对应 $A$ 使一些维度被压平但是 $\vec{b}$ 在这些维度非零，欠定方程无解的情况类似，但是 $A$ 的列空间本来就是目标空间的一个子空间，矛盾的出现相当于让 $A$ 的列空间维度与目标空间维度之差进一步增加了。

Regularization of under-determined system

Tikhonov regularization、ridge regression、guassian prior：

min_{\vec{x}} ‖ A \vec{x} - \vec{b} ‖_{2}^{2} + α ‖ \vec{x} ‖_{2}^{2}

(A^{T} A + α I_{n}) \vec{x} = A^{T} \vec{b}

$α$ 越小，poor conditioned 的可能性越大； $α$ 越大，近似解的误差越大

Lasso、Laplace prior：使用 1-范数

Elastic Net：同时使用 1-范数和 2-范数

Example: Image Alignment

Gram 矩阵和 Cholesky Factorization

对称、半正定

对于正定对阵矩阵 $C$ ，其 LU 分解可以具有更加简洁的形式： $E_{k} \dots E_{1} C E_{1}^{T} \dots E_{k}^{T} = I_{n} \Rightarrow C = E_{1}^{- 1} \dots E_{k}^{- 1} (E_{k}^{T})^{- 1} \dots (E_{1}^{T})^{- 1} = L L^{T}$

$‖ L_{i :} ‖_{2}^{2} = C_{i i}$

half storage required

衡量向量和矩阵的大小、距离

向量 p-范数和矩阵 Frobenius 范数

‖ A ‖_{Fro} = \sqrt{\sum_{i = 1}^{m} \sum_{j}^{n} A_{i j}^{2}}

矩阵诱导范数

‖ A ‖_{p} = max_{‖ x ‖_{p} = 1} ‖ A x ‖_{p}

诱导 1-范数

\begin{aligned} ‖ A ‖_{1} & = max_{‖ x ‖_{1} = 1} ‖ A x ‖_{1} \\ = max_{‖ x ‖_{1} = 1} \sum_{i = 1}^{n} | \sum_{j = 1}^{m} A_{i j} x_{j} | \\ \leq max_{‖ x ‖_{1} = 1} \sum_{j = 1}^{m} (\sum_{i = 1}^{n} | A_{i j} |) | x_{j} | \\ \leq max_{j} \sum_{i = 1}^{n} | A_{i j} | \end{aligned}

诱导 infty-范数

\begin{aligned} ‖ A ‖_{\infty} & = max_{‖ x ‖_{\infty} = 1} ‖ A x ‖_{\infty} \\ = max_{‖ x ‖_{\infty} = 1} max_{i} | \sum_{j = 1}^{m} A_{i j} x_{j} | \\ \leq max_{i} max_{‖ x ‖_{\infty} = 1} \sum_{j = 1}^{m} | A_{i j} | | x_{j} | \\ \leq max_{i} \sum_{j = 1}^{m} | A_{i j} | \end{aligned}

诱导 2-范数、谱范数

\begin{aligned} ‖ A ‖_{2}^{2} & = max_{‖ \vec{x} ‖_{2} = 1} ‖ A \vec{x} ‖_{2}^{2} \\ = max_{‖ \vec{x} ‖_{2} = 1} {\vec{x}}^{T} A^{T} A \vec{x} \\ = max_{‖ \vec{x} ‖_{2} = 1} {\vec{x}}^{T} V Σ^{T} Σ V^{T} \\ = max_{‖ \vec{x} ‖_{2} = 1} {\vec{x}}^{T} Σ^{T} Σ \vec{x} \\ = max_{1 \leq i \leq n} σ_{i}^{2} x_{i}^{2} \\ \leq max_{1 \leq i \leq n} σ_{i}^{2} \end{aligned}

$A \vec{x} = \vec{b}$ 的条件数

限制为可逆方阵 $A$ ，考虑 $\vec{b}$ 的一个扰动 $δ \vec{b}$ 满足 $A (\vec{x} + δ \vec{x}) = \vec{b} + δ \vec{b}$ 也即 $δ \vec{x} = A^{- 1} δ \vec{b}$ ，根据诱导 2-范数的性质我们有 $‖ δ \vec{x} ‖_{2} \leq ‖ A^{- 1} ‖_{2} \cdot ‖ δ \vec{b} ‖_{2}$ ，结合 $‖ \vec{b} ‖_{2} \leq ‖ A ‖_{2} \cdot ‖ \vec{x} ‖_{2}$ 得到

\frac{‖ δ \vec{x} ‖_{2}}{‖ \vec{x} ‖_{2}} \leq (‖ A ‖_{2} \cdot ‖ A^{- 1} ‖_{2}) \frac{‖ δ \vec{b} ‖_{2}}{‖ \vec{b} ‖_{2}}

求解 Overdetermined

QR

求解 $A^{T} A \vec{x} = A^{T} \vec{b}$ 时 $L U$ 分解并不好用，条件数被平方，用 $Q R$ 分解

$A^{T} A \vec{x} = A^{T} \vec{b} \Leftrightarrow R^{T} Q^{T} Q R \vec{x} = R^{T} Q^{T} \vec{b} \Leftrightarrow R^{T} R \vec{x} = R^{T} Q^{T} \vec{b}$

当 $A$ 列满秩时 $R$ 一定是可逆矩阵，我们有 $\vec{x} = R^{- 1} Q^{T} \vec{b}$ ，非列满秩时 $R$ 一定不可逆

投影

考虑标准内积 $\cdot$ ， $\vec{x}$ 向 $\vec{b}$ 投影，投影向量 $\vec{p} = k {\frac{\vec{b}}{‖ b ‖}}_{2} = k \vec{e}$ 满足

\begin{aligned} k & = \arg min_{l} ‖ l \vec{e} - \vec{x} ‖_{2}^{2} \\ = \arg min_{l} l^{2} + ‖ \vec{x} ‖_{2}^{2} - 2 l \cdot \vec{e} \cdot \vec{x} \\ = \vec{e} \cdot \vec{x} \end{aligned}

\vec{p} = \frac{\vec{b} \cdot \vec{x}}{‖ b ‖_{2}^{2}} \vec{b} ⊥ (\vec{x} - \vec{p})

向量到向量的投影推广到向量（一维空间）到高维空间的投影：考虑空间 $S = span {\vec{e_{1}}, \vec{e_{2}}, \dots, \vec{e_{n}}}$ 一组标准正交基 $\vec{e_{1}}, \vec{e_{2}}, \dots, \vec{e_{n}}$ ， $\vec{x}$ 在分量 $\vec{e_{i}}$ 上的投影向量 ${proj}_{\vec{e_{i}}} \vec{x} = (\vec{e_{i}} \cdot \vec{x}) \vec{e_{i}}$

Gram-Schmidt 正交化

Gram-Schmidt 正交化是最自然的思路，每次取出一个向量并减去它到已经生成的子空间上的投影得到新的标准正交基。

每次计算投影向量都需要和 $a_{j}$ 做内积，计算 $a_{j}$ 过程中的误差会被逐步放大。比如 $a_{1}$ 本身参与了 $k - 1$ 次运算，受到 $a_{1}$ 误差影响的 $a_{2}$ 参与了 $k - 2$ 次运算，依次类推，我们最终得到 $p_{k}$ 大约叠加了 $Θ (k^{2})$ 次有 $a_{1}$ 误差参与的计算。 $a_{j}$ 的误差分量可能在与 $v_{i}$ 的运算过程中得到放大。

Modified Gram-Schmidt 正交化，从对每个 $v_{i}$ 都计算子空间的投影到构建子空间时就直接剔除 $v_{i}$ 的分量。得到的 $a_{k}$ 取决于 $v_{k}$ ， $v_{k}$ 受到 $a_{1}, a_{2}, \dots, a_{k - 1}$ 的影响，在 $v_{k}$ 被 $a_{i}$ 修正时， $v_{k}$ 已经失去了前 $i - 1$ 个基方向上的分量，在这些方向上不会放大 $a_{i}$ 的误差

Householder QR

$\vec{x}$ 关于 ${proj}_{\vec{b}} \vec{x}$ 轴对称向量为

\begin{aligned} 2 {proj}_{\vec{b}} \vec{x} - \vec{x} & = \frac{2 \vec{b} ({\vec{b}}^{T} \cdot \vec{x})}{‖ b ‖_{2}^{2}} - \vec{x} \\ = \frac{2 (\vec{b} \cdot {\vec{b}}^{T}) \vec{x}}{{\vec{b}}^{T} \cdot \vec{b}} - \vec{x} \\ = (\frac{2 (\vec{b} \cdot {\vec{b}}^{T})}{{\vec{b}}^{T} \cdot \vec{b}} - I) \vec{x} \\ = - H_{\vec{b}} \vec{x} \end{aligned}

$\vec{x}$ 关于镜面（法向量为 $\vec{b}$ ，单位法向量为 $\vec{v}$ ）的镜像向量为 $x - 2 {proj}_{\vec{b}} \vec{x} = H_{\vec{b}} \vec{x} = H_{\vec{v}} \vec{x} = (I - 2 \vec{v} \cdot {\vec{v}}^{T}) \vec{x}$

Eigen

Basics

$A$ 可对角化的充要条件是其特征空间生成 $R^{n}$ 。 $R$ 内 $A$ 可正交对角化的充要条件是 $A$ 为实对称矩阵。

对角化能极大简化求逆

PCA

最小化垂直重构误差等价于最大化投影方差。

判断数据之间的相关性（PCA 主成分分析）：

min_{\vec{v}} \sum_{i} ‖ \vec{x_{i}} - {proj}_{\vec{v}} \vec{x_{i}} ‖_{2}^{2} s . t . ‖ {\vec{v}}_{2} ‖ = 1

等价于

\max_{\vec{v}}\lVert X^T \vec{v} \rVert _{2}^2 \quad s.t. \lVert \vec{v} \rVert {#2_} {2} = 1

等价于求解 $X X^{T}$ 的主特征向量。求出的 $\vec{v}$

主特征向量

采用迭代法 $A^{k} \vec{v} \approx λ_{1}^{k} c_{1} {\vec{x}}_{1}$ 而不是解 $det (A - λ I) = 0$ 然后解 $(A - λ I) \vec{v} = \vec{0}$

Power iteration 使用迭代 $\vec{ω_{k}} = A \vec{v_{k - 1}}, \vec{v_{k}} = \vec{ω_{k}} / ‖ \vec{ω_{k}} ‖$ 求最大特征值。
Inverse iteration 使用迭代 $\vec{ω_{k}} = A^{- 1} \vec{v_{k - 1}}, \vec{v_{k}} = \vec{ω_{k}} / ‖ \vec{ω_{k}} ‖$ 求最小特征值，求逆用 $A = L U$ 代替。
Shifted $(A - σ I) \vec{v} = (λ - σ) \vec{v}$

配合 Inverse iteration 求最接近于 $σ$ 的特征值
当 $\frac{λ_{2}}{λ_{2}}$ 接近 $1$ 时用 $\frac{λ_{2} - σ}{λ_{1} - σ}$ 来替代，加速收敛
Rayleigh Quotient Iteration 把 ${\arg \min}_{λ} ‖ A \vec{v_{0}} - λ \vec{v_{0}} ‖_{2}^{2} = \frac{{\vec{v_{0}}}^{T} \vec{A} v_{0}}{‖ \vec{v_{0}} ‖_{2}^{2}}$ 作为 $σ$ 的逐轮最优估计值来加速迭代。

SVD

{argcrit}_{\vec{v}} ‖ A \vec{v} ‖_{2} / ‖ \vec{v} ‖_{2}

等价于

{argcrit}_{\vec{v}} ‖ A \vec{v} ‖_{2}^{2} s . t . ‖ \vec{v} ‖_{2} = 1

等价于求 $A^{T} A$ 的特征向量。

$A^{T} A$ 是对阵半正定矩阵，所以所有特征值都 $\geq 0$ 并所有特征向量正交
$A^{T} A$ 的一组特征向量 $v_{i}$ 和 $A A^{T}$ 的一组特定特征向量 $u_{i}$ 可以通过 $\vec{u} = A \vec{v}$ 建立一一映射
$‖ \vec{u_{i}} ‖_{2} = \sqrt{λ_{i}}$ ，利用这一点归一化 $u_{i}$
${\vec{u_{i}}}^{T} A \vec{v_{j}} = \sqrt{λ_{i}} δ_{i j}$
$A \in R^{m \times n}$ ， $k = card ({λ_{i}})$ ，我们可以自然地构造 $U Σ V^{T} : (m \times k) \cdot (k \times k) \cdot (k \times n)$ ，如果 $k$ 比较小，我们可以通过补全正交的零特征值特征向量得到 $U Σ V^{T} : (m \times m) \cdot (m \times n) \cdot (n \times n)$ ，补全后 $U$ 的列向量称为左奇异向量， $Σ$ 的对角线元素称为奇异值

Example: Overdetermined Equation and Pseudo Inverse

求解超定方程时施加 2-范数限制

min ‖ \vec{x} ‖_{2}^{2} s . t . A^{T} A \vec{x} = A^{T} \vec{b}

限制条件等价于 $V Σ^{T} Σ V^{T} \vec{x} = V Σ^{T} U^{T} \vec{b} \Leftrightarrow Σ^{T} Σ (V^{T} \vec{x}) = Σ^{T} (U^{T} \vec{b})$ ，令 $\vec{y} = V^{T} \vec{x}, \vec{d} = U^{T} \vec{b}$ 转化为 $σ_{i}^{2} y_{i} = σ d_{i}$ ，取 $Σ_{i j}^{+} = \frac{1}{σ_{i}} δ_{i j} if σ_{i} \neq 0 else C_{i} δ_{i j} where C_{i} is arbitrary$ 那么有 $\vec{x} = V \vec{y} = V Σ^{+} U^{T} \vec{b}$ ，其中 $A^{+} = V Σ^{+} U^{T}$ 称为 $A$ 的伪逆。

伪逆对于超定方程给出最小化 2-范数的解，对欠定方程给出在 2-范数意义下的最佳近似解。

Example: Low Rank Approximation

求 $A \vec{x}$ 近似解时，可以舍弃较小 $σ_{i}$ 的谱；求 $A^{+} \vec{x}$ 的近似解时可以舍弃较大 $σ_{i}$ 的谱

$Eckart-Young-Mirsky$ 定理，用之秩至多为 $k$ 的矩阵 $\tilde{A}$ 来逼近 $A$ ，在 Frobenius 范数和谱范数的意义下，按照奇异值大小取前 $k$ 大的谱作为 $\tilde{A}$ 是所有谱组合中的的最优解

Proof？

Example: Least Square with Tikonov Regularization

Tikonov 正则化的最小二乘法 $(A^{T} A + α I) \vec{x} = A^{T} \vec{b}$ 的解 $\vec{x} = V D U^{T} \vec{b}, D_{i i} = σ_{i} / (σ_{i}^{2} + α^{2})$

Detail？

Example: Rigid Alignment

min_{R^{T} R = I_{3}, \vec{t} \in R^{3}} \sum_{i} ‖ R \vec{x_{1 i}} + \vec{t} - \vec{x_{2 i}} ‖_{2}^{2}

通过最小二乘处理 $\vec{t}$

min_{R^{T} R = I_{3}} ‖ R X_{1} - X_{2}^{t} ‖_{F r o}^{2}

Orthogonal Procrustes Theorem：

${argmin}_{R^{T} R = I} ‖ A R - B ‖_{F r o}^{2} = U V^{T} where A^{T} B = U Σ V^{T}$
${argmin}_{R^{T} R = I} ‖ R A - B ‖_{F r o}^{2} = U V^{T} where B A^{T} = U Σ V^{T}$
保持 $R = U V^{T}$ 的形式分别需要 $Tr (R V Σ U^{T}) = Tr (R B^{T} A)$ 、 $Tr (R V Σ U^{T}) = Tr (R A B^{T})$

Detail？

Example: Eigenfaces(Minimum Reconstruction Error)

$C \in R^{n \times d} \land C^{T} C = I_{d \times d}$ 最小化 $‖ X - C C^{T} X ‖_{F r o}$ 的取值是 $U$ 的前 $d$ 列， $X = U Σ V^{T}$

Detail？

求根

二分法

根的存在性和算法收敛性：
- 充分条件： $f \in C^{0} [a, b]$ 且 $f (a) f (b) < 0$ ，则必定存在一个根 $x^{⋆} \in [a, b]$ 且二分法必定收敛。
一阶收敛速度

牛顿迭代法

$x_{n + 1} = x_{n} - \frac{f (x_{n})}{f^{'} (x_{n})}$
收敛性：
- 局部收敛的充分条件： $f \in C^{2} [a, b]$ ， $x^{⋆}$ 为单根（ $f (x^{⋆}) = 0 \land f^{'} (x^{⋆}) \neq 0$ ），必定存在一个 $x^{⋆}$ 的 $δ$ 邻域，只要初始值在邻域内，牛顿迭代法必定收敛。
- 全局收敛的充分条件： $f$ 在 $[a, b]$ 上单调、具有相同凹凸性，只要存在根（在此条件下等价于端点函数值异号），当 $f (x_{0}) f^{″} (x_{0}) > 0$ 时牛顿迭代必定收敛。
- 不收敛的根本原因包括 $f$ 局部性质差或者 $x_{0}$ 的选择差。
收敛速度为 $2$ 阶
- 泰勒展开即证

切线法

$x_{n + 1} = x_{n} - f (x_{n}) \cdot \frac{x_{n} - x_{n - 1}}{f (x_{n}) - f (x_{n - 1})} = \frac{f (x_{n}) x_{n - 1} - f (x_{n - 1}) x_{n}}{f (x_{n}) - f (x_{n - 1})}$
收敛性
- 局部收敛的充分条件
实际计算时在接近根时 $f (x_{n}) - f (x_{n - 1})$ 会发生灾难性抵消，数值稳定性不好
收敛速度为 $\frac{\sqrt{5} + 1}{2}$ 阶
- 记 $I (a_{1}, a_{2}, \dots, a_{n}) = [min {a_{1}, a_{2}, \dots, a_{n}}, max {a_{1}, a_{2}, \dots, a_{n}}]$
- 用一阶商差和二阶商差（首尾相接的割线斜率的变化率）改写并利用商差中值定理

\begin{aligned} e_{n + 1} & = e_{n} - \frac{f (x_{n})}{f [x_{n}, x_{n - 1}]} \\ = e_{n} - \frac{e_{n} f [x_{n}, x^{⋆}]}{f [x_{n}, x_{n - 1}]} \\ = \frac{e_{n} (f [x_{n - 1}, x_{n}] - f [x_{n}, x^{⋆}])}{f [x_{n}, x_{n - 1}]} \\ = \frac{e_{n} e_{n - 1} f [x_{n - 1}, x_{n}, x^{⋆}]}{f [x_{n}, x_{n - 1}]} \\ = e_{n} e_{n - 1} \cdot \frac{f^{″} (η_{n})}{f^{'} (ξ_{n})} where η_{n} \in I (x_{n - 1}, x_{n}, x^{⋆}), ξ_{n} \in I (x_{n - 1}, x_{n}) \end{aligned}

LU Decomposition

Basics

Example: Regression of multiple experiments

Regularization of under-determined system

Example: Image Alignment

Gram 矩阵和 Cholesky Factorization

衡量向量和矩阵的大小、距离

向量 p-范数和矩阵 Frobenius 范数

矩阵诱导范数

Ax→=b→ 的条件数

求解 Overdetermined

QR

投影

Gram-Schmidt 正交化

Householder QR

Eigen

Basics

PCA

主特征向量

SVD

Example: Overdetermined Equation and Pseudo Inverse

Example: Low Rank Approximation

Example: Least Square with Tikonov Regularization

Example: Rigid Alignment

Example: Eigenfaces(Minimum Reconstruction Error)

求根

不动点迭代求不动点

$A \vec{x} = \vec{b}$ 的条件数