CS205A Lecture 4 Designing linear systems (incl. least-squares); special structure (Cholesky, sparsity)

距离上次更新已经 2200 天了，文章内容可能已经过时。

课程主页：https://graphics.stanford.edu/courses/cs205a-13-fall/schedule.html

这次回顾第四讲，这一讲结束了LU分解的内容，介绍敏感性和病态性。

Chapter 3 设计和分析线性系统

回归

假设我们的函数有如下形式

$$f(\overrightarrow{x})=\sum _{i=1}^m{a}_i{f}_i(\overrightarrow{x})$$

如果我们有$m$个观测值，那么可以得到如下线性系统：

$$\left(\begin{array}{cccc}{f}_1({\overrightarrow{x}}^{(1)})& f_2({\overrightarrow{x}}^{(1)})& \dots & f_m({\overrightarrow{x}}^{(1)})\\ f_1({\overrightarrow{x}}^{(2)})& f_2({\overrightarrow{x}}^{(2)})& \dots & f_m({\overrightarrow{x}}^{(2)})\\ ⋮& ⋮& ⋮& ⋮\\ f_1({\overrightarrow{x}}^{(m)})& f_2({\overrightarrow{x}}^{(m)})& \dots & f_m({\overrightarrow{x}}^{(m)})\end{array}\right)\left(\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_m\end{array}\right)=\left(\begin{array}{c}{y}^{(1)}\\ y^{(2)}\\ ⋮\\ y^{(m)}\end{array}\right)$$

我们的任务是求解$(a_1,\dots ,a_m{)}^T$，在继续之前，来看两个具体例子：

多项式回归

如果

$$f(x)=\sum _{i=0}^n{a}_i{x}^i$$

那么此时称为多项式回归，对应的线性系统为：

$$\left(\begin{array}{ccccc}1& x^{(1)}& (x^{(1)}{)}^2& \dots & (x^{(1)}{)}^n\\ 1& x^{(2)}& (x^{(2)}{)}^2& \dots & (x^{(2)}{)}^n\\ ⋮& ⋮& ⋮& \dots & ⋮\\ 1& x^{(n)}& (x^{(n)}{)}^2& \dots & (x^{(n)}{)}^n\end{array}\right)\left(\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_n\end{array}\right)=\left(\begin{array}{c}{y}^{(1)}\\ y^{(2)}\\ ⋮\\ y^{(n)}\end{array}\right)$$

Mini-Fourier

另一种常见的形式为

$$f(x)=a\cos (x+\varphi )$$

最小二乘法

考虑上述问题，我们自然希望

$$A\overrightarrow{x}\approx \overrightarrow{b}$$

这等价于

$$\underset{\overrightarrow{x}}{min}||A\overrightarrow{x}-\overrightarrow{b}||$$

最小化上述问题等价于最小化

$$g(\overrightarrow{x})=||A\overrightarrow{x}-\overrightarrow{b}|{|}^2$$

化简可得

$$\begin{array}{rl}g(\overrightarrow{x})& =||A\overrightarrow{x}-\overrightarrow{b}|{|}^2\\ & =(A\overrightarrow{x}-\overrightarrow{b}{)}^T(A\overrightarrow{x}-\overrightarrow{b})\\ & ={\overrightarrow{x}}^T{A}^{T}A\overrightarrow{x}-{\overrightarrow{b}}^{T}A\overrightarrow{x}-{\overrightarrow{x}}^T{A}^T\overrightarrow{b}+{\overrightarrow{b}}^T\overrightarrow{b}\\ & ={\overrightarrow{x}}^T{A}^{T}A\overrightarrow{x}-2{\overrightarrow{x}}^T{A}^T\overrightarrow{b}+{\overrightarrow{b}}^T\overrightarrow{b}\end{array}$$

对上式关于$\overrightarrow{x}$求梯度可得

$$\begin{array}{rl}{\mathrm{\nabla }}_{\overrightarrow{x}}g(\overrightarrow{x})& ={\mathrm{\nabla }}_{\overrightarrow{x}}({\overrightarrow{x}}^T{A}^{T}A\overrightarrow{x}-2{\overrightarrow{x}}^T{A}^T\overrightarrow{b}+{\overrightarrow{b}}^T\overrightarrow{b})\\ & =2A^{T}A\overrightarrow{x}-2A^T\overrightarrow{b}\end{array}$$

令上式为$0$可得

$$A^{T}A\overrightarrow{x}=A^T\overrightarrow{b}$$

该方程被称为正规方程，下面对$A^{T}A$这个矩阵做更多的讨论。

(半)正定矩阵和Cholesky分解

讨论之前先介绍一些基本概念：

对称性

$$B\mathrm{是}\mathrm{对}\mathrm{称}\mathrm{矩}\mathrm{阵}\mathrm{如}\mathrm{果}{B}^T=B$$

显然$A^{T}A$是对称矩阵。

(半)正定性

$$\begin{gathered} B\mathrm{是}\mathrm{半}\mathrm{正}\mathrm{定}\mathrm{矩}\mathrm{阵}\mathrm{如}\mathrm{果}\mathrm{任}\mathrm{取}\overrightarrow{x}，\mathrm{我}\mathrm{们}\mathrm{有}{\overrightarrow{x}}^{T}B\overrightarrow{x}\ge 0； \\ B\mathrm{是}\mathrm{正}\mathrm{定}\mathrm{矩}\mathrm{阵}，\mathrm{如}\mathrm{果}\mathrm{当}\overrightarrow{x}\ne \overrightarrow{0}，\mathrm{我}\mathrm{们}\mathrm{有}{\overrightarrow{x}}^{T}B\overrightarrow{x}>0 \end{gathered}$$

注意到

$${\overrightarrow{x}}^T{A}^{T}A\overrightarrow{x}=(A\overrightarrow{x}{)}^T(A\overrightarrow{x})=||A\overrightarrow{x}|{|}_2^2\ge 0$$

所以$A^{T}A$是半正定矩阵，事实上，如果$A$可逆，那么$A^{T}A$是正定矩阵。

有了这些准备工作，考虑正定对称系统$C\overrightarrow{x}=\overrightarrow{d}$的求解，假设$C\in {\mathbb{R}}^{n\times n}$且$C$是对称矩阵，将$C$写为分块形式：

$$C=\left(\begin{array}{cc}{c}_{11}& {\overrightarrow{v}}^T\\ \overrightarrow{v}& \tilde{C}\end{array}\right)$$

其中$\overrightarrow{v}\in {\mathbb{R}}^{n-1},\tilde{C}\in {\mathbb{R}}^{(n-1)\times (n-1)}$。由$C$的定义，不难得到：

$$\begin{array}{rl}{\overrightarrow{e}}_1^{T}C\overrightarrow{e}& =\left(\begin{array}{cccc}1& 0& \dots & 0\end{array}\right)\left(\begin{array}{cc}{c}_{11}& {\overrightarrow{v}}^T\\ \overrightarrow{v}& \tilde{C}\end{array}\right)\left(\begin{array}{c}1\\ 0\\ \dots \\ 0\end{array}\right)\\ & =c_{11}\\ & >0\end{array}$$

其中最后一个不等号是因为正定性。回顾高斯消元法，这说明我们不需要对第一列选主元，那么使用高斯消元法，我们可以找到forward substitution矩阵$E$，其中

$$E=\left(\begin{array}{cc}1/\sqrt{c_{11}}& {\overrightarrow{0}}^T\\ \overrightarrow{r}& I_{(n-1)\times (n-1)}\end{array}\right)$$

使得

$$EC=\left(\begin{array}{cc}\sqrt{c_{11}}& {\overrightarrow{v}}^T/\sqrt{c_{11}}\\ \overrightarrow{0}& D\end{array}\right)$$

其中$D\in {\mathbb{R}}^{(n-1)\times (n-1)}$。接着，由$C$的对称性，考虑使用列变换：

$$\begin{array}{rl}EC{E}^T& =\left(\begin{array}{cc}\sqrt{c_{11}}& {\overrightarrow{v}}^T/\sqrt{c_{11}}\\ \overrightarrow{0}& D\end{array}\right)\left(\begin{array}{cc}1/\sqrt{c_{11}}& \overrightarrow{r}\\ {\overrightarrow{0}}^T& I_{(n-1)\times (n-1)}\end{array}\right)\\ & =\left(\begin{array}{cc}1& \overrightarrow{0}\\ \overrightarrow{0}& D\end{array}\right)\end{array}$$

不难发现$D$也是正定对称矩阵，对$D$重复这个过程，最终不难得到

$$\begin{gathered} (E_k...E_1)C(E_k...E_1{)}^T=I \\ C=(E_k...E_1{)}^{-1}((E_k...E_1{)}^{-1}{)}^T \end{gathered}$$

记

$$L=(E_k...E_1{)}^{-1}$$

那么

$$C=L{L}^T$$

注意到上述$L$并不唯一——对于正交矩阵$Q$，我们有

$$C=LQ{Q}^T{L}^T=(LQ)(LQ{)}^T$$

所以上述分解有无数种形式，特别的，如果我们取$L$为下三角矩阵，那么

$$C=L{L}^T$$

被称为Cholesky分解，下面推导Cholesky分解的具体形式。

由定义，我们有

$$\left(\begin{array}{cccc}{L}_{11}& 0& \dots & 0\\ L_{21}& L_{22}& \dots & 0\\ ⋮& ⋮& \ddots & 0\\ L_{n1}& L_{n2}& \dots & L_{nn}\end{array}\right)\left(\begin{array}{cccc}{L}_{11}& L_{21}& \dots & L_{n1}\\ 0& L_{22}& \dots & L_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & L_{nn}\end{array}\right)=C$$

考虑$C_{ij}$，我们有

$$C_{ij}=\sum _{k=1}^j{L}_{ik}{L}_{kj}=\sum _{k=1}^j{L}_{ik}{L}_{jk},j\le i$$

所以

$$L_{ij}=(C_{ij}-\sum _{k=1}^{j-1}{L}_{ik}{L}_{jk})/L_{jj},j<i$$

以及

$$L_{ii}=\sqrt{C_{ii}-\sum _{k=1}^{i-1}{L}_{ik}^2},j=i$$

所以，对$i,j$做循环即可计算得到$L_{ij}$，伪代码如下：

// Takes as input an n-by -n matrix A[i,j]
// Edits A in place to obtain the Cholesky factorization in its lower triangle
for i from 1 to n {
    // Back - substitute to find l_i
    for j from 1 to i -1 { // element j of l_i
        sum = 0;
        for k from 1 to j -1
    		sum += A[i,k]*A[j,k];
    	A[i,j] = (A[i,j]-sum )/A[j,j];
    }
    // Apply the formula for l_kk
    normSquared = 0
    for k from 1 to i -1
    	normSquared += A[i,k ]^2;
    A[i,i] = sqrt (A[i,i] - normSquared );
}