LU Facorization and FactorUpdate

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.

// 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}");

Ouput

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.

// 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}");

Ouput

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.