layout: true class: typo, typo-selection --- count: false class: nord-dark, middle, center LDLT Factorization ======================================= @luk036 👨‍💻 · 2026 📅 --- Cholesky factorization --------------------------- - Let $A$ be a symmetric positive matrix in $\mathbb{R}^{n\times n}$. - The Cholesky factorization of $A$ is the factorization $A = L L^T$, where $L$ is lower triangular with positive diagonal elements. - For example, $$\left[\begin{array}{cccc} 1 & 1 & 1 & 1 \\ 1 & 2 & 1 & 2 \\ 1 & 1 & 3 & 1 \\ 1 & 2 & 1 & 4 \end{array}\right] = \left[\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 0 & \sqrt{2} & 0 \\ 1 & 1 & 0 & \sqrt{2} \end{array}\right] \left[\begin{array}{cccc} 1 & 1 & 1 & 1 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & \sqrt{2} & 0 \\ 0 & 0 & 0 & \sqrt{2} \end{array}\right].$$ --- ### 🔷 Matrix Factorization Flow .mermaid[ <pre> graph TD A[Symmetric Matrix A] --> B{Choose method} B --> C[Cholesky A = LL^T] B --> D[LDLT A = LDL^T] C --> E[L: lower triangular] C --> F[+ sqrt operations] D --> G[L: unit diagonal] D --> H[D: diagonal matrix] D --> I[No sqrt needed] style A fill:#ff9800 style B fill:#f44336 style C fill:#2196f3 style D fill:#4caf50 style E fill:#2196f3 style F fill:#2196f3 style G fill:#4caf50 style H fill:#4caf50 style I fill:#4caf50 </pre> ] --- Basic algorithm ---------------- - Start with $a_{11} = l_{11}^2$, the basic algorithm of Cholesky factorization is: $$\begin{array}{l} 1~ \text{for}~i=1:n \\ 2~ \quad \text{for}~ j=1:i-1 \\ 3~ \qquad s = a_{ij} - \sum_{k=1}^{j-1} l_{ik} l_{jk} \\ 4~ \qquad l_{ij} = s / l_{jj} \\ 5~ \quad \text{end} \\ 6~ \quad l_{ii} = \bigl(a_{ii} - \sum_{k=1}^{i-1} l_{ik}^2\bigr)^{1/2} \\ 7~ \text{end} \end{array} $$ - Stop as soon as indefiniteness is detected (line 6). - Invoke $\frac{p^3}{2}$ FLOP's, where $p$ is the place the algorithm stops. - Invoke square root operations --- Witness of indefiniteness ------------------------- - In the case of failure, a vector $v$ can be constructed to certify that $v^T A v < 0$. - Let $L_{1:p}$ denote the partial sub-matrix $L(1:p, 1:p)$ where $p$ is the row of failure. - Then $v = [L_{1:p}^{-T} e_p, 0, \cdots, 0]^T$, where $e_p = [0, \cdots, 0, 1]^T \in \mathbb{R}^p$ - Start with $v = [0, \cdots, 0, l^{-1}_{11}]^T$, the basic algorithm is: $$\begin{array}{l} 1~ \text{for}~i = p - 1~\text{downto}~1 \\ 2~ \quad \text{for}~ k = i~\text{to}~p \\ 3~ \qquad v_i = v_i - l_{k,i} v_k \\ 4~ \quad \text{end} \\ 5~ \quad v_i = v_i / l_{ii} \\ 6~ \text{end} \end{array} $$ --- LDLT factorization ------------------ - The LDLT factorization of a symmetric positive definite matrix $A$ is the factorization $A = L D L^T$, where $L$ is lower triangular with unit diagonal elements and $D$ is a diagonal matrix. - For example, $$\left[\begin{array}{cccc} 1 & 1 & 1 & 1 \\ 1 & 2 & 1 & 2 \\ 1 & 1 & 3 & 1 \\ 1 & 2 & 1 & 4 \end{array}\right] = \left[\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 1 & 0 & 1 \end{array}\right] \left[\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 2 & 0 \\ 0 & 0 & 0 & 2 \end{array}\right] \left[\begin{array}{cccc} 1 & 1 & 1 & 1 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array}\right].$$ --- ### ⚡ Cholesky vs LDLT .mermaid[ <pre> graph LR subgraph Cholesky C1[p^3/2 FLOPs] --> C2[sqrt ops] C2 --> C3[A = LL^T] end subgraph LDLT L1[p^3/2 FLOPs] --> L2[No sqrt] L2 --> L3[A = LDL^T] end style C3 fill:#2196f3 style L3 fill:#4caf50 style C1 fill:#2196f3 style C2 fill:#2196f3 style L1 fill:#4caf50 style L2 fill:#4caf50 </pre> ] --- Naïve implementation -------------------- - Then, start with $a_{11} = d_{11}$, the basic algorithm of LDLT factorization is: $$\begin{array}{l} 1~ \text{for}~i=1:n \\ 2~ \quad \text{for}~ j=1:i-1 \\ 3~ \qquad s = a_{ij} - \sum_{k=1}^{j-1} d_{kk} l_{ik} l_{jk} \\ 4~ \qquad l_{ij} = s / d_{jj} \\ 5~ \quad \text{end} \\ 6~ \quad d_{ii} = a_{ii} - \sum_{k=1}^{j-1} d_{kk} l_{ik} l_{jk} \\ 7~ \text{end} \end{array} $$ - Invoke $p^3$ FLOP's, where $p$ is the place the algorithm stops. --- Storage representation ---------------------- First, we pack the solution and the intermediate storage on a single matrix $T$ such that: $$t_{ij} = \begin{cases} d_{ii} & \mathrm{if}\ i = j, \\ l_{ij} & \mathrm{if}\ i > j, \\ d_{ii} l_{ji} & \mathrm{if}\ j < i. \end{cases} $$ - For example, $$T = \left[\begin{array}{cccc} d_{11} & d_{11} l_{21} & d_{11} l_{31} & d_{11} l_{41} \\ l_{21} & d_{22} & d_{22} l_{32} & d_{22} l_{42} \\ l_{31} & l_{32} & d_{33} & d_{33} l_{43} \\ l_{41} & l_{42} & l_{43} & d_{44} \end{array}\right].$$ --- Improved implementation ----------------------- - Then, start with $a_{11} = t_{11}$, the improved implementation of LDLT factorization is: $$\begin{array}{l} 1~ \text{for}~i=1:n \\ 2~ \quad \text{for}~ j=1:i-1 \\ 3~ \qquad t_{ji} = a_{ij} - \sum_{k=1}^{j-1} t_{ik} t_{jk} \\ 4~ \qquad t_{ij} = t_{ji} / t_{jj} \\ 5~ \quad \text{end} \\ 6~ \quad t_{ii} = a_{ii} - \sum_{k=1}^{i-1} t_{ik} t_{ki} \\ 7~ \text{end} \end{array} $$ - Invoke $\frac{p^3}{2}$ FLOP's (same as Cholesky factorization's), where $p$ is the place the algorithm stops. --- Witness of indefiniteness ------------------------- - In the case of failure, a vector $v$ can be constructed to certify that $v^T A v < 0$. - Let $L_{1:p}$ denote the partial sub-matrix $L(1:p, 1:p)$ where $p$ is the row of failure. - Then $v = [L_{1:p}^{-T} e_p, 0, \cdots, 0]^T$, where $e_p = [0, \cdots, 0, 1]^T \in \mathbb{R}^p$ - Start with $v = e_p$, the basic algorithm is: $$\begin{array}{l} 1~ \text{for}~i = p - 1~\text{downto}~1 \\ 2~ \quad \text{for}~ k = i~\text{to}~p \\ 3~ \qquad v_i = v_i - t_{k,i} v_k \\ 4~ \quad \text{end} \\ 5~ \text{end} \end{array} $$ --- --- ### 🛡️ Indefiniteness Detection .mermaid[ <pre> graph TD A[Start factorization] --> B{diag > 0?} B -->|Yes| C[Continue next row] B -->|No| D[Stop at row p] D --> E[Construct witness v] E --> F[v^T A v < 0 proven] C --> G[Complete factorization] style A fill:#ff9800 style B fill:#f44336 style C fill:#2196f3 style D fill:#f44336 style E fill:#9c27b0 style F fill:#f44336 style G fill:#4caf50 </pre> ] --- count: false class: nord-dark, middle, center Q&A 🎤 =========