Cholesky factorization, sparse matrices
Structured Linear Systems
Gaussian elimination and/or LU can solve all the example problems above. But these systems can have special properties that make them easier or stabler to solve.
Today's example: Positive definite, sparsity.
min square problem: \[ \min_{x} ||Ax-b||^2_2 \rightarrow A^TAx = A^Tb \] \[ M:=A^TA, \] 1. Symmetric \[ M^T = (A^TA)^T = A^TA = M \] 2. Positive (Semi-)Definite
B is positive semidefinite if for all \(x \in \mathbb{R}^n, x^TMx \geq 0\).
B is positive definite if \(x^TMx > 0\) whenever \(x \neq 0\).
\[ \text{Take } v \in \mathbb{R}^n, v^TMv = v^TA^TAv = ||Av||^2_2 \geq 0 \]
A ridiculously important Matrix
\[ A^TA \] "Gram matrix"
Pivoting for SPD \(C\)
Solve \(Cx = d\) for \(C\) SPD.
\[ C = \begin{bmatrix} c_{11} & v^T \\ v & \tilde{C} \end{bmatrix}, v \in \mathbb{R}^{n-1} \]
Forward substitution
\[ E=\begin{bmatrix} \frac{1}{\sqrt{c_{11}}} & 0^T \\ r & I_{(n-1)\times (n-1)} \end{bmatrix} \] It can be a problem if \(c_{11} < 0\).
Proposition: \(c_{11} > 0\). Proof: \[ e_1^TCe_1 = c_{11} > 0, \] \(e_1\) is the first standard basis vector,\(e_1 = [1, 0, \cdots, 0]^T\).
\[ EC = \begin{bmatrix} \frac{1}{\sqrt{c_{11}}} & 0^T \\ r & I_{(n-1)\times (n-1)} \end{bmatrix} \begin{bmatrix} c_{11} & v^T \\ v & \tilde{C} \end{bmatrix} = \begin{bmatrix} \sqrt{c_{11}} & \frac{v^T}{\sqrt{c_{11}}} \\ c_{11}r + v & rv^T + \tilde{C} \end{bmatrix} \]
By definition, since \(EC\) is the result of a forward substitution, it is lower triangular. So \(c_{11}r + v = 0\).
To maintain symmetry, we can post-multiply by \(E^T\).
\[ E^TCE = \begin{bmatrix} \sqrt{c_{11}} & \frac{v^T}{\sqrt{c_{11}}} \\ 0 & rv^T + \tilde{C} \end{bmatrix} \begin{bmatrix} \frac{1}{\sqrt{c_{11}}} & r^T \\ 0 & I_{(n-1)\times (n-1)} \end{bmatrix} = \begin{bmatrix} 1 & \sqrt{c_{11}}r^T + \frac{v^T}{\sqrt{c_{11}}} \\ 0 & rv^T + \tilde{C} \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & rv^T + \tilde{C} \end{bmatrix} \]
Since \(E^TCE\) is symmetric, \(\sqrt{c_{11}}r^T + \frac{v^T}{\sqrt{c_{11}}}\) is zero.
What enabled this?
- Positive definiteness $$ existence of \(\sqrt{c_{11}}\)
- Symmetry \(\rightarrow\) apply \(E\) to both sides.
正定矩阵可以不是对称的。
Cholesky Factorization
\[ E_1^T E_2^T \cdots E_n^T C E_n \cdots E_2 E_1 = I\\ C = LL^T, L = E_1^{-1}E_2^{-1} \cdots E_n^{-1} \] If we replace \(c_{11}\) with \(1\), The 1 in the first row and column of \(ECE^T\) will be another value. So finally we get a diagonal matrix instead of a identity matrix. \[ C = LDL^T\] Reason to do this: 1. Get rid of square roots 2. \(c_{11}\) might be close to zero.
separate row \(k\) and column \(k\).
\[ L= \begin{bmatrix} L_{11} & \vec{0} & 0 \\ \vec{l_k^T} & l_{kk} & \vec{0} \\ L_{31} & \vec{l_{k}'} & L_{33} \end{bmatrix} \]
\[ \begin{aligned} LL^T &= \begin{bmatrix} L_{11} & \vec{0} & 0 \\ \vec{l_k^T} & l_{kk} & \vec{0} \\ L_{31} & \vec{l_{k}'} & L_{33} \end{bmatrix} \begin{bmatrix} L_{11}^T & \vec{l_k} & L_{31}^T \\ \vec{0} & l_{kk} & \vec{l_k'}^T \\ 0 & \vec{0} & L_{33}^T \end{bmatrix}\\ &= \begin{bmatrix} .. & .. & .. \\ \vec{l_k}^T L_{11}^T & ||\vec{l_k}||^2 + l_{kk}^2 & ..\\ .. & .. & .. \end{bmatrix}\\ &= \begin{bmatrix} c_{11} & .. & .. \\ c_{k}^T & c_{kk} & ..\\ c_{31} & c_{k}' & .. \end{bmatrix}\\ \end{aligned} \]
\[ L_{11}l_k = c_k\text{:LT Solve}, O((k-1)^2)\\ c_{kk} = ||\vec{l_k}||^2 + l_{kk}^2 \\ \Rightarrow l_{kk}^2 = c_{kk} - ||\vec{l_k}||^2 \\ \Rightarrow l_{kk} = \sqrt{c_{kk} - ||\vec{l_k}||^2}\\ \]
\(l_{kk}\)取正负都可以，但是为了方便一般取正的。
算法是逐行进行的。首先\(L_{11}\)可以直接求出，然后从第二行开始，根据公式依次：
- 求解下三角线性方程组，\(L_{11}l_k = c_k, O((k-1)^2)\)得到\(l_k\)。
- 计算\(l_{kk} = \sqrt{c_{kk} - ||\vec{l_k}||^2}\)
因此整个算法的复杂度为\(O(n^3)\)。n行乘以每行的复杂度\(O(n^2)\)。
Interpretation of Cholesky
What is \(x^TCx\)?
\[ x^TCx = x^TLL^Tx = ||L^Tx||^2_2 \geq 0 \]
Storing Sparse Matrices
Want \(O(n)\) storage if we have \(O(n)\) nonzeros.
Examples:
- List of triplets \((r, c, val)\)
- For each row, matrix[r] holds a dictionary \(c \rightarrow A[r][c]\)
直接使用高斯消元法的话，会导致稀疏矩阵变成稠密矩阵。
Avoiding fill
- Common strategy: Permute rows/columns
- Mostly heuristics constructions
- Minimizing fill in Cholesky is NP-complete
- Alternative strategy: Avoid Gaussian altogether
Band Matrices
Cyclic Matrices
Cholesky factorization, sparse matrices