2022-02-01

Computational

稀疏矩阵分解

稀疏因子

$$
\delta(A)=\frac{\sum_{i,j}[A_{i,j} \ne 0]}{n\cdot m}
$$

当 $\delta \le 0.05$ 时，认为矩阵是稀疏的

稀疏矩阵存储

行优先存储（CSR）

列优先存储（CSC）

填入元(fill-in)

高斯消元中的新增位置

$$
\begin{bmatrix}
\color{red}X & \color{red}X & \color{red}X \\
X & & X \\
X & X &
\end{bmatrix}
\Rightarrow
\begin{bmatrix}
\color{red}X & \color{red}X & \color{red}X \\
\color{green}X & \color{blue}X & \color{green}X \\
\color{green}X & \color{green}X & \color{blue}X
\end{bmatrix}
$$

Cholesky分解中的update矩阵

$$
B \Rightarrow B-\frac{\boldsymbol{v}\boldsymbol{v}^T}{x}
$$

求解方阵方程

$$
\begin{cases}
A\boldsymbol{x}=\boldsymbol{b} \\
A=LU
\end{cases}
\Rightarrow
LU\boldsymbol{x}=\boldsymbol{b}
\Rightarrow
\begin{cases}
L\boldsymbol{y}=\boldsymbol{b} \\
U\boldsymbol{x}=\boldsymbol{y}
\end{cases}
$$

只需要解决两个问题

$A=LU$
$\boldsymbol{y}=L^{-1}\boldsymbol{b}$

重排序

$$
A\boldsymbol{x}=\boldsymbol {b} \Leftrightarrow PAP^TP\boldsymbol{x}=P\boldsymbol{b} \\
PAP^T=LU
$$

消去图中按照度数从小到大排序

Markowitz积，按照度数平方排序，即fill-in个数

嵌套割集(Nested dissection)

每次取当前图的平均割集，并放到最后

$$
A - S - B
\Rightarrow
\begin{bmatrix}
A & & X \\
&B& X\\
X &X& S \\
\end{bmatrix}
$$

近似最小度(Appr. min. degree, AMD)

带宽缩小排序(Reverse Cuthill-McKee, Sloan, …)

非零元向对角线居中

Lsolve算法

用于解决 $\boldsymbol{y}=L^{-1}\boldsymbol{b}$

$$
y(j)=\frac{b(j)-\sum_{i=1}^{j-1}L(j,i) \cdot y(i)}{L(j,j)}
$$

y = b
forall(i = 1:n, y(i) /= 0)
    y(i) = y(i) / L(i, i)    ! 如果没有单位对角化L
    forall(j = i+1:n, L(j, i) /= 0)
        y(j) = y(j) - L(j, i) * y(i)
    end
end

重排序？若 $L(j, i) \ne 0$，则先遍历 $j$ 后遍历 $i$

若 $L(j,i) \ne 0$，连边 $j \to i$，由于 $j \ge i$，所以得到一张DAG，找拓扑序即可

假设 $L$ 已经进行主对角线单位化，则 $y(j)=b(j)-\sum_{i=1}^{j-1}L(j,i) \cdot y(i)$

考虑fill-in的关系：

$$
\begin{cases}
b(j) \ne 0 \Rightarrow y(j) \ne 0 \\
y(i) \ne 0 \land L(j,i) \ne 0 \Rightarrow y(j) \ne 0
\end{cases}
$$

只需要计算：$\mathcal{B}={i \mid b(i) \ne 0},\mathcal{Y}={i \mid y(i) \ne 0}$

其中 $\mathcal{Y}=\mathrm{Reach}_{G(L)}(\mathcal{B})$

这一部分的时间复杂度是 $O(|\mathcal{B}|+f)$

Left-looking算法

用于解决 $A=LU$

考虑从左到右依次确定 $L,U$ 的前 $k$ 列的值

$$
\begin{bmatrix}
A_{11} & a_{12} & A_{13} \\
a_{21} & a_{22} & a_{23} \\
A_{31} & a_{32} & A_{33}
\end{bmatrix}
=
\begin{bmatrix}
L_{11} & & \\
l_{21} & 1 & \\
L_{31} & l_{32} & L_{33}
\end{bmatrix}
\begin{bmatrix}
U_{11} & u_{12} & U_{13} \\
& u_{22} & U_{23} \\
& & U_{33}
\end{bmatrix}
$$

归纳性的，假设已经计算完前 $k-1$ 列的LU分解，考虑计算第 $k$ 列

$$
\begin{bmatrix}
a_{12} \\
a_{22} \\
a_{32}
\end{bmatrix}
=
\begin{bmatrix}
L_{11} & & \\
l_{21} & 1 & \\
L_{31}& l_{32} & L_{33} \\
\end{bmatrix}
\begin{bmatrix}
u_{12} \\
u_{22} \\
0 \\
\end{bmatrix}
$$

这里 $l_{32},L_{33},u_{12},u_{22}$ 都是未知的，同时要求求出 $l_{32},u_{12},u_{22}$

考虑到 $a_{32}=L_{31} \cdot u_{12}+l_{32} \cdot u_{22}$，所以只需要用 $l_{32}=\frac{a_{32}-L_{31}\cdot u_{12}}{u_{22}}$ 计算 $l_{32}$ 即可

即解方程：

$$
\begin{bmatrix}
a_{12} \\
a_{22} \\
a_{32}
\end{bmatrix}
=
\begin{bmatrix}
L_{11} & & \\
l_{21} & 1 & \\
L_{31}& 0 & I \\
\end{bmatrix}
\begin{bmatrix}
x_1 \\
x_2 \\
x_3 \\
\end{bmatrix}
$$

其中：

$$
\begin{cases}
u_{11}=x_1 \\
u_{22}=x_2 \\
l_{32}=\frac{x_3}{x_2}
\end{cases}
$$

初始情况：$\boldsymbol{a} _ 1=\boldsymbol{l} _ 1 \cdot u _ {11}$，其中 $L _ {11}=1,u _ {11}=A _ {11},L _ {k1}=\frac{A _ {k1}}{A _ {11}}$

Up-looking算法

用于解决实对称矩阵Cholesky分解

由对称性可知，只需要沿主对角线依次求 $L$ 的 $k$ 阶顺序主子阵 $L_k$

$$
L _ {k-1}\hat{\boldsymbol{l}} _ k=\boldsymbol{a} _ k(1:k-1) \\
\left\Vert\begin{bmatrix}
\hat{\boldsymbol{l}} _ k \\
L_{kk}
\end{bmatrix}\right\Vert^2
=a_{kk}
\Rightarrow
L _ {kk}=\sqrt{a _ {kk}-\hat{\boldsymbol{l}} _ k^T \cdot \hat{\boldsymbol{l}} _ k} \\
\hat{\mathcal{L}} _ k=\mathrm{Reach} _ {G_{k-1}}(\mathcal{A} _ k)
$$

实对称矩阵Cholesky分解

$$
A=A^T \Rightarrow A=LL^T
$$

考虑：

$$
A=\begin{bmatrix}
x & \boldsymbol{v}^T \\
\boldsymbol{v} & B
\end{bmatrix}=
\begin{bmatrix}
1 & \\
\frac{\boldsymbol{v}}{x} & I
\end{bmatrix}
\begin{bmatrix}
x & \boldsymbol{v}^T \\
& B-\frac{\boldsymbol{v}\boldsymbol{v}^T}{x}
\end{bmatrix}
=
\begin{bmatrix}
1 & \\
\frac{\boldsymbol{v}}{x} & I
\end{bmatrix}
\begin{bmatrix}
x & \\
& B-\frac{\boldsymbol{v}\boldsymbol{v}^T}{x}
\end{bmatrix}
\begin{bmatrix}
1 & \frac{\boldsymbol{v}^T}{x} \\
& I
\end{bmatrix}
$$