Numerical Linear Algebra with Iterative Methods
\(\det(A) \approx 0\) may say it is “almost” singular, but it is not scale-invariant
The condition number \(\kappa\), given matrix norm \(||\cdot||\) uses the matrix norm
\[ \text{cond}(A) \equiv \|A\| \|A^{-1}\|\geq 1 \]
Expensive to calculate, can show that given spectrum
\[ \text{cond}(A) = \left|\frac{\lambda_{max}}{\lambda_{min}}\right| \]
Intuition: if \(\text{cond}(A) = K\), then \(b \to b + \nabla b\) change in \(b\) amplifies to a \(x \to x + K \nabla b\) error when solving \(A x = b\).
See Matlab Docs on inv for why inv is a bad idea when \(\text{cond}(A)\) is huge
Stationary Iterative Methods reorganize the problem so it is a contraction mapping and then iterate
For matrices that are strictly diagonal dominant \[ |A_{ii}| \geq \sum_{j\neq i} |A_{ij}| \quad\text{for all } i = 1\ldots N \]
Note for our problem rows sum to 0 so if \(\rho > 0\) then \(\rho I - Q\) is strictly diagonally dominant
\[ A \equiv D + R \]
\[ D \equiv \begin{bmatrix} A_{11} & 0 & \ldots & 0\\ 0 & A_{22} & \ldots & 0\\ \vdots & \vdots & \vdots & \vdots\\ 0 & 0 & \ldots & A_{NN} \end{bmatrix}\,\, R \equiv \begin{bmatrix} 0 & A_{12} & \ldots & A_{1N} \\ A_{21} & 0 & \ldots & A_{2N} \\ \vdots & \vdots & \vdots & \vdots\\ A_{N1} & A_{N2} & \ldots & 0 \end{bmatrix} \]
Where \(D^{-1}\) is trivial since diagonal. To solve, take an iteration \(x^k\), starting from \(x^0\),
\[ x^{k+1} = D^{-1}(b - R x^k) \]
\[ \min_{\beta} \| X \beta -y \|^2 + \alpha \| \beta\|^2 \]
M = 1000
N = 10000
sigma = 0.1
beta = rand(M)
# simulate data
X = sprand(N, M, 0.1)
y = X * beta + sigma * randn(N)
beta_direct = X \ y
results = lsmr(X, y, log = true)
beta_lsmr = results[1]
@show norm(beta_direct - beta_lsmr)
println("$(results[end])")norm(beta_direct - beta_lsmr) = 1.0615184241157149e-5
Converged after 15 iterations.N = 4
a = 0.1
b = 0.05
Q = Tridiagonal(fill(b, N-1),
[-a; fill(-(a + b), N-2); -b],
fill(a, N-1))
# Find smallest magnitude eigenvalue (i.e. 0)
lambda, phi = eigs( Q', nev = 1, which = :SM, maxiter = 1000)
phi = real(phi) ./ sum(real(phi))
@show lambda
@show mean(phi);
Q'lambda = ComplexF64[-9.270990685613062e-18 + 0.0im]
mean(phi) = 0.254×4 Tridiagonal{Float64, Vector{Float64}}:
-0.1 0.05 ⋅ ⋅
0.1 -0.15 0.05 ⋅
⋅ 0.1 -0.15 0.05
⋅ ⋅ 0.1 -0.05LinearMap wrapper adds features required for algorithms (e.g. size(Q_adj_map and Q_adj_map * v overloads)\[ \begin{aligned} A x &= b\\ A P^{-1} P x &= b\\ A P^{-1} y &= b\\ P x &= y \end{aligned} \] That is, solve \((A P^{-1})y = b\) for \(y\), and then solve \(P x = y\) for \(x\).