LU-разложение

Материал из Википедии — свободной энциклопедии

LU-разложение — представление матрицы $A$ в виде $L U$ , где $L$ — нижняя треугольная матрица с единичной диагональю, а $U$ — верхняя треугольная. LU-разложение еще называют LU-факторизацией.

Алгоритм LU-разложения лежит в основе широко распространенного метода решения систем линейных алгебраических уравнений(СЛАУ).

Матрица $L$ является нижнетреугольной с единичной диагональю, поэтому её определитель равен 1. Матрица $U$ — верхнетреугольная матрица, значит ее определитель равен произведению элементов, расположенных на главной диагонали.

Будем использовать следующие обозначения для элементов матриц $L = (l i j)$ , $U = (u i j)$ , $i,j = 1\dots n$ ; причём диагональные элементы матрицы $L$ : $l i i = 1$ , $i = 1\dots n$ . Тогда, если известно LU-разложение матрицы, её определитель можно вычислить по формуле $det(A) = det(L U) = det(L)det(U) = det(U)$ = произведению эементов на диагонали матрицы U.

Найти матрицы $L$ и $U$ можно следующим образом(выполнять шаги следует строго по порядку, так как следующие элементы находятся с использованием предыдущих):

$u_{1j} = a_{1j},\ j = 1\dots n$
$l_{j1} = \frac {a_{j1}} {u_{11}},\ j = 2\dots n\ (u_{11} \ne 0)$

Для $i = 2\dots n$

$u_{ij} = a_{ij} - \sum_{k=1}^{i-1} l_{ik}u_{kj},\ j = i\dots n$
$l_{ji} = \frac {1} {u_{ii}}(a_{ji} - \sum_{k=1}^{i-1} l_{jk}u_{ki}),\ j = i+1\dots n$

В итоге мы получим матрицы — $L$ и $U$ . В программной реализации данного метода (компактная схема Гаусса) для представления матриц $L$ и $U$ можно обойтись всего одним массивом, в котором совмещаются матрицы $L$ и $U$ . Например вот так(для матрицы размером $3\times 3$ ):

$\begin{pmatrix} u_{11} & u_{12} & u_{13} \\ l_{21} & u_{22} & u_{23} \\ l_{31} & l_{32} & u_{33} \\ \end{pmatrix}$

Фрагмент программы на C# для нахождения матриц $L$ и $U$ .

//переменная n(размерность иcходной ''квадратной'' матрицы) должна получить значение до этого момента
            double[,] A = new double[n, n];
            double[,] L = new double[n, n];
            double[,] U = new double[n, n];
//до этого момента массив A должен быть полностью определён
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    U[0, i] = A[0, i];
                    L[i, 0] = A[i, 0] / U[0, 0];
                    double sum = 0;
                    for (int k = 0; k < i; k++)
                    {
                        sum += L[i, k] * U[k, j];
                    }
                    U[i, j] = A[i, j] - sum;
                    if (i > j)
                    {
                        L[j, i] = 0;
                    }
                    else
                    {
                        sum = 0;
                        for (int k = 0; k < i; k++)
                        {
                            sum += L[j, k] * U[k, i];
                        }
                        L[j, i] = (A[j, i] - sum) / U[i, i];
                    }
                }
            }
//после выполнения цикла в массиве L — элементы матрицы L, в массиве U — элементы матрицы U.
 
//Теперь можно вычислять определитель

Фрагмент программы на matlab для нахождения матриц $L$ и $U$ (оптимизированный вариант программы для C#). Существует встроенная функция matlab [L,U]=lu(A), она считает немного иначе. Данный пример можно использовать для перевода алгоритма на другие языки программирования.

function [L,U] = lu_my(A)
    n = size(A, 1); %размерность матрицы A
 
    %Инициализация матриц U и L
    U = zeros(n);
    L = eye(n);
    for i = 1:n
        U(1, i) = A(1, i);
        L(i, 1) = A(i, 1) / U(1, 1);
    end
 
    %Вычисление матриц U и L    
    for i = 2:n
        for j = i:n
            sumu = 0;
            suml = 0;
            for k = 1:i
                sumu = sumu + L(i, k) * U(k, j);
                if i < j, suml = suml + L(j, k) * U(k, i); end
            end
            U(i, j) = A(i, j) - sumu;
            if i < j, L(j, i) = (A(j, i) - suml) / U(i, i); end        
        end
    end
end