1. Learning Objectives

By the end of this lesson, you will be able to:

• Implement and apply the Cholesky decomposition to generate correlated random returns for Monte Carlo simulation.

• Derive and apply the QR decomposition for solving linear systems and computing regression coefficients.

• Understand the LU decomposition and its application in solving systems of equations efficiently.

• Compute the pseudo-inverse using SVD and apply it to solve ill-conditioned linear systems.

• Apply the Woodbury Matrix Identity to efficiently invert large covariance matrices in factor models.

• Implement shrinkage estimators (Ledoit-Wolf) to produce well-conditioned covariance matrices.

• Understand the computational complexity of matrix operations and their implications for large-scale financial AI.

---

2. Cholesky Decomposition – Generating Correlated Random Variables

The Cholesky decomposition factors a symmetric positive definite matrix A into the product of a lower triangular matrix L and its transpose:
A = L L^T.

• L is lower triangular with positive diagonal entries.

• This is the matrix analogue of the square root of a positive number.

2.1 The Algorithm (Forward Substitution)
For A = {a_{ij}} with i, j = 1, ..., n:
L_{jj} = sqrt( a_{jj} - Σ_{k=1}^{j-1} L_{jk}^2 ).
L_{ij} = (1 / L_{jj}) * ( a_{ij} - Σ_{k=1}^{j-1} L_{ik} L_{jk} ) for i > j.
All other entries are zero.

2.2 Financial Application – Generating Correlated Asset Returns
We need to generate N correlated normal returns with covariance matrix Σ.

1. Compute Cholesky factor L such that Σ = L L^T.

2. Generate N independent standard normal variables z ∼ N(0, I).

3. Compute r = L z. Then Cov(r) = L Cov(z) L^T = L I L^T = Σ.

4. For log-normal prices: S_t = S_0 exp( (μ - σ²/2)Δt + σ sqrt(Δt) L z ).

2.3 AI Application – Monte Carlo Simulation for Option Pricing
To price a path-dependent option (e.g., Asian option):

1. Generate M paths of correlated asset returns using Cholesky.

2. Compute the payoff for each path.

3. Discount the average payoff: V = e^{-rT} * (1/M) Σ_{i=1}^{M} Payoff_i.
   Cholesky ensures the correlations between assets are preserved, which is essential for basket options and multi-asset derivatives.

2.4 Computational Complexity
Cholesky decomposition requires O(n³/3) floating-point operations. This is half the cost of LU decomposition. It is the preferred method for positive definite matrices.

---

3. LU Decomposition – Solving Systems of Linear Equations

LU decomposition factors a square matrix A into the product of a lower triangular matrix L and an upper triangular matrix U:
A = L U.

• L has ones on its diagonal.

• U is upper triangular.

3.1 The Algorithm (Doolittle’s Method)
For each column j = 1, ..., n:

1. Compute U_{ij} = A_{ij} - Σ_{k=1}^{j-1} L_{ik} U_{kj} for i ≥ j.

2. Compute L_{ij} = (1 / U_{jj}) * ( A_{ij} - Σ_{k=1}^{j-1} L_{ik} U_{kj} ) for i > j.

3.2 Solving Ax = b using LU

1. Factor A = L U.

2. Solve L y = b using forward substitution: y_1 = b_1 / L_{11}}, y_i = b_i - Σ_{j=1}^{i-1} L_{ij} y_j.

3. Solve U x = y using backward substitution: x_n = y_n / U_{nn}}, x_i = (y_i - Σ_{j=i+1}^{n} U_{ij} x_j) / U_{ii}}.

3.3 Financial Application – Portfolio Optimisation
The optimal portfolio weights w = Σ^{-1} μ require solving the linear system Σ w = μ.
Instead of inverting Σ (which is O(n³)), we can:

1. Perform LU decomposition of Σ: Σ = L U.

2. Solve L U w = μ using forward/backward substitution (O(n²) per solution).
   If we need to solve for multiple μ vectors (e.g., for different target returns), LU decomposition only needs to be done once.

3.4 Pivoting
If a zero pivot is encountered during LU decomposition, we perform row permutations (pivoting). This yields P A = L U, where P is a permutation matrix. In finance, covariance matrices are often ill-conditioned, and pivoting is essential for numerical stability.

---

4. QR Decomposition – Orthogonalisation and Regression

QR decomposition factors a matrix A ∈ R^{m x n} (with m ≥ n) into:
A = Q R
where:

• Q ∈ R^{m x m} is an orthogonal matrix (Q^T Q = I).

• R ∈ R^{m x n} is an upper triangular matrix.

4.1 Methods for QR Decomposition

• Gram-Schmidt Process: Orthogonalises columns iteratively. Numerically unstable.

• Householder Reflectors: Reflects vectors to zero out sub-diagonal entries. Stable and widely used.

• Givens Rotations: Rotates vectors to zero out entries. Used for sparse matrices.

4.2 Solving Least Squares Problems using QR
The least squares problem minimises ||Ax - b||_2. The solution is x = (A^T A)^{-1} A^T b.
Using QR: A = Q R. Then A^T A = R^T Q^T Q R = R^T R.
The normal equations become R^T R x = R^T Q^T b.
Since R is upper triangular, we solve R x = Q^T b using backward substitution.
Advantage: QR is numerically more stable than solving the normal equations directly (A^T A can be ill-conditioned).

4.3 Financial Application – Factor Model Estimation (Fama-French)
Estimate factor loadings (betas) by regressing asset returns on factor returns:
R_i = α_i + β_{i,1} F_1 + β_{i,2} F_2 + ... + β_{i,K} F_K + ε_i.
In matrix form: R = X β + ε.
The OLS estimator is β = (X^T X)^{-1} X^T R.
QR decomposition of X gives a stable estimate, even when factors are correlated (multicollinearity).

4.4 Financial Application – Kalman Filter
The Kalman filter is used for state estimation (e.g., tracking latent volatility). It requires solving least squares problems at each time step. QR decomposition is used in the Square Root Kalman Filter to maintain numerical stability over long horizons.

---

5. Singular Value Decomposition (SVD) – The Master Decomposition (Extended)

We covered SVD in Lesson 2.1. Here we focus on its practical applications in finance.

5.1 Pseudo-Inverse (Moore-Penrose)
For a matrix A = U Σ V^T, the pseudo-inverse is:
A^+ = V Σ^+ U^T
where Σ^+ is the diagonal matrix with 1/σ_i for non-zero singular values and 0 otherwise.
Application: Solve A x = b even when A is singular or rectangular. x = A^+ b gives the minimum norm solution.

5.2 Financial Application – Ill-Conditioned Covariance Matrices
In high-dimensional finance (N > T), the sample covariance matrix Σ_sample = (1/(T-1)) R^T R is singular (rank ≤ T < N). We cannot invert it.
Using SVD: R = U Σ V^T. Then Σ_sample = (1/(T-1)) V Σ^T Σ V^T.
The pseudo-inverse is: Σ_sample^+ = (T-1) V (Σ^T Σ)^+ V^T.
This allows us to compute portfolio weights even when N > T. However, the results are unstable and require regularisation.

5.3 Low-Rank Approximation (Eckart-Young)
The best rank-k approximation of A is A_k = U_k Σ_k V_k^T.
Financial Application: Factor models. If R is the returns matrix, its low-rank approximation R_k is the factor component. The residual R - R_k is idiosyncratic noise. This is the basis for Random Matrix Theory (RMT) denoising of covariance matrices.

5.4 Principal Component Analysis (PCA) via SVD
PCA is equivalent to SVD on the centered data matrix:
R_c = U Σ V^T.

• V are the principal components (eigenvectors of R_c^T R_c).

• U Σ are the principal component scores.
  Application: Feed the first K principal component scores into an LSTM instead of the raw returns. This reduces dimensionality and denoises the data.

---

6. The Woodbury Matrix Identity – Efficient Inversion

The Woodbury identity allows us to compute the inverse of a matrix that has been perturbed by a low-rank update:
(A + U C V)^{-1} = A^{-1} - A^{-1} U (C^{-1} + V A^{-1} U)^{-1} V A^{-1}.

6.1 Financial Application – Factor Models
Recall the factor model covariance matrix: Σ = B Λ B^T + Ψ, where Ψ is diagonal.
Using Woodbury:
Σ^{-1} = Ψ^{-1} - Ψ^{-1} B (Λ^{-1} + B^T Ψ^{-1} B)^{-1} B^T Ψ^{-1}.
If N = 1000 assets and K = 5 factors, the inversion cost reduces from O(N³) = 1e9 operations to O(K³) = 125 operations plus O(N K²) ≈ 25,000 operations. This is a speedup of ~40,000x.
AI Application: In reinforcement learning for portfolio optimisation, the agent needs to compute the inverse of the covariance matrix at every time step. The Woodbury identity makes this feasible.

6.2 Financial Application – Shrinkage Estimators
A shrinkage estimator combines the sample covariance S with a structured target T (e.g., the identity matrix or the one-factor model):
Σ_shrink = δ T + (1-δ) S.
The inverse can be computed efficiently using Woodbury if T is diagonal or low-rank.

---

7. Shrinkage Estimators – Conditioning the Covariance Matrix

The sample covariance matrix S is unbiased but has high variance, especially when N is large relative to T. Shrinkage reduces estimation error by biasing towards a structured target.

7.1 The Ledoit-Wolf Shrinkage Estimator
This is the most widely used shrinkage estimator in finance. It finds the optimal δ that minimises the Frobenius norm between the shrinkage estimator and the true covariance matrix.
Σ_shrink = δ * (Tr(S) / N) I + (1-δ) S.
The optimal δ is estimated from the data.
Implementation in Python: from sklearn.covariance import LedoitWolf; lw = LedoitWolf().fit(R); Σ_lw = lw.covariance_.

7.2 Benefits

• Positive Definite: Σ_shrink is always positive definite, even when N > T.

• Reduced Estimation Error: The condition number (ratio of largest to smallest eigenvalue) is reduced, making inversion stable.

• Improved Portfolio Performance: Shrinkage portfolios often outperform sample covariance portfolios out-of-sample.

7.3 AI Application – Neural Network Inputs
Instead of feeding raw returns to an AI model, feed the Ledoit-Wolf covariance matrix (or its Cholesky factor) as a feature. This provides a stable, well-conditioned representation of market risk.

---

8. Power Iteration and Eigenvalue Computation

In many financial applications (e.g., PCA for large datasets), we only need the largest eigenvalues and eigenvectors.

8.1 Power Iteration
For a matrix A, the power iteration converges to the eigenvector corresponding to the largest eigenvalue:

1. Initialise v_0 randomly.

2. For t = 1, 2, ...:

   • w_t = A v_{t-1}}.

   • v_t = w_t / ||w_t||_2.

   • λ_t = v_t^T A v_t.
     Convergence is linear with rate |λ_2 / λ_1|.

8.2 Financial Application – First Principal Component
The first principal component explains the largest fraction of variance. Power iteration on R^T R or R R^T gives the first eigenvector (the market factor) without computing the full SVD.
Implementation: Use scipy.sparse.linalg.eigs or sklearn.decomposition.PCA with svd_solver='randomized'.

8.3 Arnoldi Iteration (Lanczos for Symmetric Matrices)
For computing the top K eigenvalues and eigenvectors, Arnoldi iteration (or Lanczos for symmetric matrices) is used. This is the basis for scipy.sparse.linalg.eigsh.

---

9. Computational Complexity Summary

---

10. Summary for the AI Practitioner

1. Cholesky is used to generate correlated random returns for Monte Carlo simulation. Σ = L L^T, then r = L z.

2. LU decomposition solves linear systems Σ w = μ efficiently. Use forward/backward substitution instead of matrix inversion.

3. QR decomposition is numerically stable for regression problems (Fama-French factor estimation). Use R x = Q^T b.

4. SVD computes the pseudo-inverse for singular matrices (when N > T). Also used for low-rank approximation (factor models) and denoising.

5. Woodbury Identity inverts large covariance matrices in factor models. Σ^{-1} = Ψ^{-1} - Ψ^{-1} B (Λ^{-1} + B^T Ψ^{-1} B)^{-1} B^T Ψ^{-1}.

6. Ledoit-Wolf Shrinkage produces a stable, positive definite covariance matrix even when N > T. Essential for portfolio optimisation.

7. Power Iteration computes the first principal component without the full SVD. Used for the market factor.

8. Complexity matters. Choose algorithms based on the size of your financial dataset. For large portfolios, use Woodbury and low-rank methods.

---