---
title: "Numerical Linear Algebra"
source: https://www.jemoka.com/posts/kbhnumerica_linear_algebra/
---

FLOP A FLOP is a &ldquo;floating-point operations&rdquo;. One addition, subtraction, multiplication, or division of two floating-point numbers.
Basic Linear Algebra Subroutines BLAS Level 1 Vector-Vector operations
inner product \(x^{T}y\) which is \(2n-1\) flops sum \(x+y\) and scale \(\alpha x\) which is \(n\) flops BLAS Level 2 Matrix-vector product
\(y = Ax\) with \(A \in \mathbb{R}^{m \times n}\)
standard multiply \(m\qty(2n - 1)\) \(2N\) if \(A\) is sparse with \(N\) non-zero elements \(2p\qty(n+m)\) if \(A\) is given as \(A = UV^{T}\), where \(U \in \mathbb{R}^{m \times p}\) and \(V \in \mathbb{R}^{n \times p}\), which is called a &ldquo;banded&rdquo; matrix (sparse everywhere except for a strip nexht to diagonal) BLAS Level 3 Matrix-matrix operations
for \(A \in \mathbb{R}^{m\times n}\) and \(B \in \mathbb{R}^{n \times p}\)
\(mp\qty(2n-1)\) flops (on the order of \(2 mnp\)) less if \(A\) and or \(B\) is false \(\frac{1}{2}m\qty(m+1)\qty(2n-1) \approx m^{2}n\) if \(m = p\) and \(AB\) is symmetric BLAS Speeds CPU single thread does at around 1-10 GFlops/s CPU multi-tread 10-100 GFlops/s GPU typically 100 GFlops/s - 1 TFloaps/s GPUs uses &ldquo;blocking&rdquo; (divide and concur)
Solving Linear Equations solution of \(Ax = b\) is \(x = A^{-1}b\)
no one actually materializes \(A^{-1}\) for general methods, \(O\qty(n^{3})\) for \(A\) with half-bandwidth \(k\), we have \(O\qty(k^{2}n)\) For special structures, we can have:
special matricies lower-triangular &mdash; forward substitution upper-triangular &mdash; do it backwards orthogonal matrices &mdash; \(2n^{2}\) flops to compute with householder transformation factorization Factor \(A\) as a product of simple matricies \(A = A_1 &hellip; A_{k}\). And then we can compute:
\begin{equation} x = A^{-1}b = A_{k}^{-1} \dots A_{1}^{-1} b \end{equation}
Cost of solver is dominated by this factoring instead of actually the solving. We can solve a modest number of equations with the same coefficinets at the same cost as one (because you are just factoring most of the time, you can just cache the factored results).
LU Factorization Every nonsingular matrix \(A\) can be factored as:
\begin{equation} A = PLU \end{equation}
where \(P\) is a permutation, \(L\) is lower triangular, \(U\) is upper triangular. Solving procedure:
LU factorization permutation (shuffle \(X\) around backwards) forward substitution backward substitution Sparse LU Factorization For \(A\) sparse and invertible, we can factor it into four matricies \(A = P_1 L U P_2\).
Doing these row and column permutations allows \(L\) and \(U\) to be sparse. They are then chosen heuristically to yield sparse \(L\), \(U\).
Cholesky Factorization Suppose \(A\) is positive definite.
Factor \(A\) into \(A = LL^{T}\). \(L\) is lower-triangular with positive diagonal entries.
Sparse Cholesky Factorization for sparse positive definite \(A\), factor as \(A = PLL^{T}P^{T}\) adding permutation matrix \(P\) makes sparser \(L\) This is same as:
permuting rows and columns of \(A\) to get \(\tilde{A} = P^{T}AP\) and then Cholesky Interestingly we can figure out \(P\) just with the sparsity pattern of \(A\), with no value.
LDL-transpose Factorization Every nonsingular symmetric matrix \(A\) can be factored as:
\begin{equation} A = PLDL^{T} P^{T} \end{equation}
with \(P\) permutaiton matrix, \(L\) lower triangular, \(D\) block diagonal with \(1 \times 1\) or \(2 \times 2\) diagonal blocks.
Structured Sub-Blocks For:
\begin{equation} \mqty(A_{11} &amp; A_{12} \\ A_{21} &amp;A_{22}) \mqty(x_1 \\ x_2) = \mqty(b_1 \\ b_2) \end{equation}
We can first solve:
Form \(A_{11}^{-1}A_{12}\) and \(A_{11}^{-1}b_1\) (factor \(A_{11}\), solve with each column of \(A_{12}\) and \(b\)) Form Shure&rsquo;s complement: \(S = A_{22} - A_{21}A^{-1}_{11} A_{12}\) and \(\tilde{b} = b_2 - A_{21} A_{11}^{-1}b_1\) Determine \(x_2\) by solving \(S x_2 = \tilde{b}\) Determine \(x_1\) by solving \(A_{11}x_1 = b_1 - A_{12} x_2\) Matrix Inversion Lemma \begin{equation} \qty(A+BC)^{-1} = A^{-1} - A^{-1}B\qty(I + CA^{-1}B)^{-1}CA^{-1} \end{equation}
