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

LUP-разложение (LUP-декомпозиция) — представление данной матрицы $A$ в виде произведения $PA = LU$ где матрица $L$ является нижнетреугольной с единицами на главной диагонали, $U$ — верхнетреугольная общего вида, а $P$ — т. н. матрица перестановок получаемая из единичной матрицы путём перестановки строк или столбцов. Такое разложение можно осуществить для любой невырожденной матрицы. LUP-разложение используется для вычисления обратной матрицы по компактной схеме, вычисления решения системы линейных уравнений. По сравнению с алгоритмом LU-разложения алгоритм LUP-разложения может обрабатывать любые невырожденные матрицы и при этом обладает более высокой устойчивостью.

Алгоритм LUP-разложения

Пусть $A=(a_{ij})$, $L=(l_{ij})$, $U=(u_{ij})$. На практике как правило вместо матрицы перестановок P используют вектор перестановок получаемый из вектора $p_{i}=(1, 2, ..., n)$ путём перестановки элементов соответствующих номерам строк переставляемых в матрице P. Например, если

$P = \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \\ \end{pmatrix}$

то $p = (2, 1, 3)$ так как матрица P получена путём перестановки первой и второй строки. Вычисление LUP-разложения ведётся в несколько шагов. Пусть матрица C = A. На каждом i-м шаге сначала производится поиск опорного (ведущего) элемента — максимального по модулю элемента среди элементов i-го столбца находящихся не выше i-й строки — после чего строка с опорным элементом меняется местами с i-й строкой. Одновременно производится такой же обмен в матрице P. При этом если матрица невырождена то опорный элемент гарантированно будет отличен от нуля. После этого все элементы текущего i-го столбца находящиеся ниже i-й строки делятся на опорный. Далее из всех элементов $c_{jk}$ находящихся ниже i-й строки и i-го столбца (то есть таких что j>i и k>i) вычитается произведение $c_{ji}c_{ik}$. После этого счётчик i увеличивается на единицу и процесс повторяется пока i<n где n — размерность исходной матрицы. После того как все шаги будут выполнены матрица C будет представлять собой следующую сумму:

$C = L+U-E$

где E — единичная матрица.

В алгоритме используется три вложенных линейных цикла так что общую сложность алгоритма можно оценить как O(n³).

Реализация алгоритма на языке С++

Ниже представлен программный код приведённого выше алгоритма на языке С++. Здесь Matrix — некоторый контейнер поддерживающий операцию индексирования. Обратите внимание что отсчёт ведётся с нуля а не с единицы.

void LUP(const Matrix &A, Matrix &C, Matrix &P) {
    //n - размерность исходной матрицы
    const int n = A.Rows();
 
    C = A;
 
    //загружаем в матрицу P единичную матрицу
    P = IdentityMatrix();
 
    for( int i = 0; i < n; i++ ) {
        //поиск опорного элемента
        double pivotValue = 0;
        int pivot = -1;
        for( int row = i; row < n; row++ ) {
            if( fabs(C[ row ][ i ]) > pivotValue ) {
                pivotValue = fabs(C[ row ][ i ]);
                pivot = row;
            }
        }
        if( pivotValue == 0 ) {
            error("Матрица вырождена");
            return;
        }
 
        //меняем местами i-ю строку и строку с опорным элементом
        P.SwapRows(pivot, i);
        C.SwapRows(pivot, i);
        for( int j = i+1; j < n; j++ ) {
            C[ j ][ i ] /= C[ i ][ i ];
            for( int k = i+1; k < n; k++ ) 
                C[ j ][ k ] -= C[ j ][ i ] * C[ i ][ k ];
        }
    }
}
 
//теперь матрица C = L + U - E

Литература

• Кормен, Т., Лейзерсон, Ч., Ривест, Р., Штайн, К. Алгоритмы: построение и анализ = Introduction to Algorithms / Под ред. И. В. Красикова. — 2-е изд. — М.: Вильямс, 2005. — 1296 с. — ISBN 5-8459-0857-4

См. также

• LU-разложение
• Разложение матрицы