LU Facorization and FactorUpdate — SepalSolver

1. LU Factorization: The Matrix Backbone

LU Factorization (Lower-Upper) is a method of decomposing a square matrix $A$ into the product of two triangular matrices.

A = LU

L (Lower Triangular): Has 1s on the diagonal and non-zero elements below the diagonal. It stores the "elimination steps."
U (Upper Triangular): Has non-zero elements on and above the diagonal. It is the "reduced" form of the matrix.

Why use it? Instead of re-calculating the inverse of a matrix (which is computationally expensive), you solve $Ly = b$ and $Ux = y$ using forward and backward substitution.

Example 1C#

// Create a sample 3x3 matrix
Matrix A = new double[,]
{
    { 4, 3, 2 },
    { 3, 2, 1 },
    { 2, 1, 3 }
};
// Perform LU decomposition
// P is the permutation matrix (for pivoting)
A.MakeLU();
Console.WriteLine($"Matrix A:{A}");
Console.WriteLine($"Matrix L:{A.L_lu}");
Console.WriteLine($"Matrix U:{A.U_lu}");
Console.WriteLine($"Matrix P:{A.pi}");

OutputFrom the book

Matrix A:
 4   3   2 
 3   2   1 
 2   1   3 

Matrix L:
   1.0000    0.0000    0.0000
   0.5000    1.0000    0.0000
   0.7500    0.5000    1.0000

Matrix U:
   4.0000    3.0000    2.0000
   0.0000   -0.5000    2.0000
   0.0000    0.0000   -1.5000

Matrix P:0,	2,	1

2. Rank-1 Updates: Adjusting the Matrix

In many real-time systems (like structural engineering or machine learning), the matrix $A$ changes slightly over time. A Rank-1 Update happens when we add the outer product of two vectors $u$ and $v$ to the original matrix:

\tilde{A} = A + uv^T

If we already have the LU factorization of $A$ , it is wasteful to re-compute the LU factorization of $\tilde{A}$ from scratch. Re-computing takes $O(n^3)$ operations, while updating takes only $O(n^2)$ .

3. The Sherman-Morrison Formula

While LU handles the decomposition, the Sherman-Morrison formula tells us how the inverse changes after a Rank-1 update:

(A + uv^T)^{-1} = A^{-1} - \frac{A^{-1}uv^TA^{-1}}{1 + v^TA^{-1}u}

4. usage in SepalSolverPython

This script demonstrates a Rank-1 update and compares the result to the standard matrix addition.

Example 2C#

// Create a sample 3x3 matrix
Matrix A = new double[,]
{
    { 4, 3, 2 },
    { 3, 2, 1 },
    { 2, 1, 3 }
};
// Vectors for Rank-1 Update
ColVec U = new double[] { 1, 0, 2 };
RowVec V = new double[] { 2, -1, 0 };
// Outer product (Rank-1 matrix)
Matrix UV = U*V;
// Updated Matrix
Matrix A_updated = A + UV;
// Make LU and then Update
A.MakeLU(); A.LU_Rank1Update(UV);
Console.WriteLine($"Matrix A:{A}");
Console.WriteLine($"Matrix L:{A.L_lu}");
Console.WriteLine($"Matrix U:{A.U_lu}");
Console.WriteLine($"Matrix P:{A.pi}");
// Verify with LU again
A_updated.MakeLU();
Console.WriteLine($"Updated Matrix A_tilde:{A_updated}");
Console.WriteLine($"Matrix L:{A_updated.L_lu}");
Console.WriteLine($"Matrix U:{A_updated.U_lu}");
Console.WriteLine($"Matrix P:{A_updated.pi}");

OutputFrom the book

Matrix A:
 6   2   2 
 3   2   1 
 6  -1   3 

Matrix L:
   1.0000    0.0000    0.0000
   1.0000    1.0000    0.0000
   0.5000   -0.3333    1.0000

Matrix U:
   6.0000    2.0000    2.0000
   0.0000   -3.0000    1.0000
   0.0000    0.0000    0.3333

Matrix P:0,	2,	1
Updated Matrix A_tilde:
 6   2   2 
 3   2   1 
 6  -1   3 

Matrix L:
   1.0000    0.0000    0.0000
   1.0000    1.0000    0.0000
   0.5000   -0.3333    1.0000

Matrix U:
   6.0000    2.0000    2.0000
   0.0000   -3.0000    1.0000
   0.0000    0.0000    0.3333

Matrix P:0,	2,	1

5. Applications in Industry

Optimization (Quasi-Newton Methods): Used to update Hessian approximations (BFGS algorithm).
Power Systems: When a transmission line trips, the admittance matrix undergoes a Rank-1 or Rank-2 update.
Signal Processing: Adaptive filtering where the correlation matrix is updated with each new data point.