LU Algorithm
The LU algorithm, also known as LU decomposition or LU factorization, is a widely-used numerical method in linear algebra and applied mathematics for solving linear systems of equations. This algorithm involves the decomposition of a square matrix into two lower and upper triangular matrices, denoted as L and U respectively. These matrices are then used to simplify the complex system of linear equations, making it easier to compute the solution. LU decomposition has numerous applications in various fields, such as engineering, physics, and computer science, for solving problems related to partial differential equations, optimization, and large-scale simulations.
The key idea behind the LU algorithm is to perform Gaussian elimination on the original matrix, which helps to eliminate elements below the main diagonal and create a triangular structure. During this process, the elimination steps are captured by the L matrix, which contains multipliers used for elimination, and the U matrix, which holds the resulting upper triangular matrix. Once the decomposition is complete, the original matrix A can be expressed as the product of L and U, i.e., A = LU. To solve the linear system Ax = b, two simpler triangular systems are formed: Ly = b and Ux = y. These systems can be solved efficiently using forward and backward substitution, ultimately yielding the solution vector x. This method is advantageous for large-scale problems, as it reduces the computational complexity and enables efficient use of memory resources.
using System;
namespace Algorithms.Numeric.Decomposition
{
/// <summary>
/// LU-decomposition factors the "source" matrix as the product of lower triangular matrix
/// and upper triangular matrix.
/// </summary>
public static class LU
{
/// <summary>
/// Performs LU-decomposition on "source" matrix.
/// Lower and upper matrices have same shapes as source matrix.
/// Note: Decomposition can be applied only to square matrices.
/// </summary>
/// <param name="source">Square matrix to decompose.</param>
/// <returns>Tuple of lower and upper matrix.</returns>
/// <exception cref="ArgumentException">Source matrix is not square shaped.</exception>
public static (double[,] L, double[,] U) Decompose(double[,] source)
{
if (source.GetLength(0) != source.GetLength(1))
{
throw new ArgumentException("Source matrix is not square shaped.");
}
var pivot = source.GetLength(0);
var lower = new double[pivot, pivot];
var upper = new double[pivot, pivot];
for (var i = 0; i < pivot; i++)
{
for (var k = i; k < pivot; k++)
{
double sum = 0;
for (var j = 0; j < i; j++)
{
sum += lower[i, j] * upper[j, k];
}
upper[i, k] = source[i, k] - sum;
}
for (var k = i; k < pivot; k++)
{
if (i == k)
{
lower[i, i] = 1;
}
else
{
double sum = 0;
for (var j = 0; j < i; j++)
{
sum += lower[k, j] * upper[j, i];
}
lower[k, i] = (source[k, i] - sum) / upper[i, i];
}
}
}
return (L: lower, U: upper);
}
/// <summary>
/// Eliminates linear equations system represented as A*x=b, using LU-decomposition,
/// where A - matrix of equation coefficients, b - vector of absolute terms of equations.
/// </summary>
/// <param name="matrix">Matrix of equation coefficients.</param>
/// <param name="coefficients">Vector of absolute terms of equations.</param>
/// <returns>Vector-solution for linear equations system.</returns>
/// <exception cref="ArgumentException">Matrix of equation coefficients is not square shaped.</exception>
public static double[] Eliminate(double[,] matrix, double[] coefficients)
{
if (matrix.GetLength(0) != matrix.GetLength(1))
{
throw new ArgumentException("Matrix of equation coefficients is not square shaped.");
}
var pivot = matrix.GetLength(0);
var upperTransform = new double[pivot, 1]; // U * upperTransform = coefficients
var solution = new double[pivot]; // L * solution = upperTransform
(double[,] l, double[,] u) = LU.Decompose(matrix);
for (var i = 0; i < pivot; i++)
{
double pivotPointSum = 0;
for (var j = 0; j < i; j++)
{
pivotPointSum += upperTransform[j, 0] * l[i, j];
}
upperTransform[i, 0] = (coefficients[i] - pivotPointSum) / l[i, i];
}
for (var i = pivot - 1; i >= 0; i--)
{
double pivotPointSum = 0;
for (var j = i; j < pivot; j++)
{
pivotPointSum += solution[j] * u[i, j];
}
solution[i] = (upperTransform[i, 0] - pivotPointSum) / u[i, i];
}
return solution;
}
}
}
