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

В численном анализе и линейной алгебре , разложение по нижнему и верхнему ( LU ) или факторизация множит матрицу как произведение нижней треугольной матрицы и верхней треугольной матрицы. Иногда продукт также включает в себя матрицу перестановок . LU-разложение можно рассматривать как матричную форму исключения Гаусса . Компьютеры обычно решают квадратные системы линейных уравнений с использованием LU-разложения, и это также ключевой шаг при обращении матрицы или вычислении определителя матрицы. LU-разложение было введено польским математиком Тадеушем Банахевичем в 1938 году.

Определения

Разложение LDU матрицы Уолша

Пусть A - квадратная матрица. LU разложение относится к факторизации А , с правильной строки и / или столбца упорядочений или перестановок, на два фактора - нижняя треугольная матрица L и верхней треугольной матрицы U :

${\ displaystyle A = LU.}$

В нижней треугольной матрице все элементы над диагональю равны нулю, в верхней треугольной матрице все элементы под диагональю равны нулю. Например, для матрицы A 3 × 3 ее LU-разложение выглядит следующим образом:

${\ displaystyle {\ begin {bmatrix} a_ {11} & a_ {12} & a_ {13} \\ a_ {21} & a_ {22} & a_ {23} \\ a_ {31} & a_ {32} & a_ {33} \ end {bmatrix}} = {\ begin {bmatrix} \ ell _ {11} & 0 & 0 \\\ ell _ {21} & \ ell _ {22} & 0 \\\ ell _ {31} & \ ell _ {32} & \ ell _ {33} \ end {bmatrix}} {\ begin {bmatrix} u_ {11} & u_ {12} & u_ {13} \\ 0 & u_ {22} & u_ {23} \\ 0 & 0 & u_ {33} \ end { bmatrix}}.}$

Без надлежащего упорядочивания или перестановок в матрице факторизация может не осуществиться. Например, легко проверить (развернув умножение матриц), что . Если , то , по крайней мере , один из , и должен быть равен нулю, что означает , что либо L или U является сингулярным . Это невозможно, если A невырожден (обратим). Это процедурная проблема. Его можно удалить, просто переупорядочив строки матрицы A так, чтобы первый элемент переставленной матрицы отличался от нуля. Та же проблема на последующих этапах факторизации может быть устранена таким же образом; см. основную процедуру ниже. ${\ textstyle a_ {11} = \ ell _ {11} u_ {11}}$${\ textstyle a_ {11} = 0}$${\ textstyle \ ell _ {11}}$${\ textstyle u_ {11}}$

Факторизация LU с частичным поворотом

Оказывается, что правильной перестановки в строках (или столбцах) достаточно для факторизации LU. Факторизация LU с частичным поворотом (LUP) часто относится к факторизации LU только с перестановками строк:

${\ displaystyle PA = LU,}$

где L и U снова нижняя и верхняя треугольные матрицы, а Р представляет собой матрицу перестановок , которые, когда левая умноженное на А , сортирует строки A . Оказывается, что все квадратные матрицы могут быть разложены на множители в этой форме, и на практике факторизация численно устойчива. Это делает разложение LUP полезным методом на практике.

Факторизация LU с полным поворотом

LU разложение с полным поворота включает в себя как строк и столбцов перестановки:

${\ displaystyle PAQ = LU,}$

где L , U и Р имеют значения , определенные выше, а Q представляет собой матрицу перестановок , что сортирует столбцы A .

Нижне-диагональное-верхнее разложение (LDU)

Нижний диагональное-верхний (LDU) разложение является разложением вида

${\ displaystyle A = LDU,}$

где D - диагональная матрица , а L и U - унитреугольные матрицы , что означает, что все элементы на диагоналях L и U равны единице.

Прямоугольные матрицы

Выше мы требовали, чтобы A была квадратной матрицей, но все эти разложения можно обобщить и на прямоугольные матрицы. В этом случае, L и D квадратные матрицы оба из которых имеют одинаковое число строк , как A , и U имеет точно такие же размеры, что и A . Верхний треугольник следует интерпретировать как имеющий только ноль входов ниже главной диагонали, которая начинается в верхнем левом углу. Аналогичным образом , более точный термин для U является то , что она является «ступенчатый вид по строкам» матрицы A .

Пример

Факторизуем следующую матрицу 2 на 2:

${\ displaystyle {\ begin {bmatrix} 4 & 3 \\ 6 & 3 \ end {bmatrix}} = {\ begin {bmatrix} \ ell _ {11} & 0 \\\ ell _ {21} & \ ell _ {22} \ end {bmatrix}} {\ begin {bmatrix} u_ {11} & u_ {12} \\ 0 & u_ {22} \ end {bmatrix}}.}$

Один из способов найти LU-разложение этой простой матрицы - просто решить линейные уравнения путем проверки. Расширение матричного умножения дает

${\ displaystyle {\ begin {align} \ ell _ {11} \ cdot u_ {11} +0 \ cdot 0 & = 4 \\\ ell _ {11} \ cdot u_ {12} +0 \ cdot u_ {22} & = 3 \\\ ell _ {21} \ cdot u_ {11} + \ ell _ {22} \ cdot 0 & = 6 \\\ ell _ {21} \ cdot u_ {12} + \ ell _ {22} \ cdot u_ {22} & = 3. \ end {выравнивается}}}$

Эта система уравнений недоопределена . В этом случае любые два ненулевых элемента матриц L и U являются параметрами решения и могут быть установлены произвольно на любое ненулевое значение. Следовательно, чтобы найти единственное LU-разложение, необходимо наложить некоторые ограничения на L- и U- матрицы. Например, мы можем удобно потребовать, чтобы нижняя треугольная матрица L была единичной треугольной матрицей (т. Е. Установить все элементы ее главной диагонали равными единицам). Тогда система уравнений имеет следующее решение:

${\ displaystyle {\ begin {align} \ ell _ {11} = \ ell _ {22} & = 1 \\\ ell _ {21} & = 1.5 \\ u_ {11} & = 4 \\ u_ {12 } & = 3 \\ u_ {22} & = - 1,5 \ end {align}}}$

Подставляя эти значения в разложение LU выше, получаем

${\ displaystyle {\ begin {bmatrix} 4 & 3 \\ 6 & 3 \ end {bmatrix}} = {\ begin {bmatrix} 1 & 0 \\ 1.5 & 1 \ end {bmatrix}} {\ begin {bmatrix} 4 & 3 \\ 0 & -1,5 \ конец {bmatrix}}.}$

Существование и уникальность

Квадратные матрицы

Любая квадратная матрица допускает факторизации LUP и PLU . Если это обратимый , то он допускает LU (или LDU ) разложение тогда и только тогда , когда все ее ведущие главные миноры отличны от нуля. Если - особая матрица ранга , то она допускает LU- факторизацию, если первые ведущие главные миноры отличны от нуля, хотя обратное неверно. ${\ textstyle A}$${\ textstyle A}$${\ textstyle A}$${\ textstyle k}$${\ textstyle k}$

Если квадратная обратимая матрица имеет LDU (факторизация со всеми диагональными элементами L и U, равными 1), то факторизация уникальна. В этом случае факторизация LU также уникальна, если мы требуем, чтобы диагональ (или ) состояла из единиц. ${\ textstyle L}$${\ textstyle U}$

Симметричные положительно определенные матрицы

Если является симметричным (или эрмитово , если сложна) положительно определенная матрица , мы можем устроить так , что U является сопряженной транспозицией из L . То есть мы можем записать A как

${\ Displaystyle A = LL ^ {*}. \,}$

Это разложение называется разложением Холецкого . Разложение Холецкого всегда существует и уникально - при условии, что матрица положительно определена. Кроме того, вычисление разложения Холецкого более эффективно и численно более стабильно, чем вычисление некоторых других LU-разложений.

Общие матрицы

Для (не обязательно обратимой) матрицы над любым полем известны точные необходимые и достаточные условия, при которых она имеет LU-факторизацию. Условия выражаются в виде рангов определенных подматриц. Алгоритм исключения Гаусса для получения разложения LU также был расширен на этот наиболее общий случай.

Алгоритмы

Закрытая формула

Когда LDU разложение существует и единственно, существует замкнутая (явная) формула для элементов L , D и U в терминах соотношения детерминант определенных подматриц исходной матрицы A . В частности, и для , этого отношение -й главной подматрицы к -му главной подматрице. Вычисление определителей требует больших вычислительных ресурсов , поэтому эта явная формула на практике не используется. ${\ textstyle D_ {1} = A_ {1,1}}$${\ textstyle i = 2, \ ldots, n}$${\ textstyle D_ {i}}$${\ textstyle i}$${\ textstyle (я-1)}$

Использование исключения Гаусса

Следующий алгоритм представляет собой модифицированную форму исключения Гаусса . Вычисление разложения LU с использованием этого алгоритма требует операций с плавающей запятой, игнорируя члены более низкого порядка. Частичное вращение добавляет только квадратичный член; это не относится к полному повороту. ${\ Displaystyle {\ tfrac {2} {3}} п ^ {3}}$

Для матрицы размера N × N определим . ${\ displaystyle A = (a_ {i, j}) _ {1 \ leq i, j \ leq N}}$${\ displaystyle A ^ {(0)}: = A}$

Мы исключаем элементы матрицы ниже главной диагонали в n -м столбце матрицы A ^{( n −1)} , добавляя к i-й строке этой матрицы n-ю строку, умноженную на

${\ displaystyle - \ ell _ {i, n}: = - {\ frac {a_ {i, n} ^ {(n-1)}} {a_ {n, n} ^ {(n-1)}} }, \ quad i = n + 1, \ dotsc, N.}$

Это можно сделать, умножив A ^{( n  - 1)} слева на нижнюю треугольную матрицу

${\ displaystyle L_ {n} = {\ begin {pmatrix} 1 &&&&& \\ & \ ddots &&&& \\ && 1 &&& \\ && - \ ell _ {n + 1, n} &&& \\ && \ vdots && \ ddots & \\ && - \ ell _ {N, n} &&& 1 \ end {pmatrix}},}$

то есть единичная матрица N × N, в которой n -й столбец заменен вектором${\ displaystyle {\ begin {pmatrix} 0 & \ dotsm & 0 & 1 & - \ ell _ {n + 1, n} & \ dotsm & - \ ell _ {N, n} \ end {pmatrix}} ^ {\textf {T} }.}$

Мы устанавливаем

${\ Displaystyle A ^ {(n)}: = L_ {n} A ^ {(n-1)}.}$

После N  - 1 шагов мы удалили все матричные элементы ниже главной диагонали, так что мы получили верхнюю треугольную матрицу A ^{( N  - 1)} . Находим разложение

${\ Displaystyle A = L_ {1} ^ {- 1} L_ {1} A ^ {(0)} = L_ {1} ^ {- 1} A ^ {(1)} = L_ {1} ^ {- 1} L_ {2} ^ {- 1} L_ {2} A ^ {(1)} = L_ {1} ^ {- 1} L_ {2} ^ {- 1} A ^ {(2)} = \ dotsm = L_ {1} ^ {- 1} \ dotsm L_ {N-1} ^ {- 1} A ^ {(N-1)}.}$

Обозначим верхнюю треугольную матрицу A ^{( N  - 1)} через U , и . Поскольку инверсия нижней треугольной матрицы L _n снова является нижней треугольной матрицей, а умножение двух нижних треугольных матриц снова является нижней треугольной матрицей, следует, что L является нижней треугольной матрицей. Более того, видно, что ${\ textstyle L = L_ {1} ^ {- 1} \ dotsm L_ {N-1} ^ {- 1}}$

${\ displaystyle L = {\ begin {pmatrix} 1 &&&&& \\\ ell _ {2,1} & \ ddots &&&& \\ & \ ddots & 1 &&& \\\ vdots & \ cdots & \ ell _ {n + 1, n} & 1 && \\ && \ vdots & \ ddots & \ ddots & \\\ ell _ {N, 1} & \ cdots & \ ell _ {N, n} & \ cdots & \ ell _ {N, N-1} & 1 \ end {pmatrix}}.}$

Мы получаем ${\ displaystyle A = LU.}$

Понятно, что для того, чтобы этот алгоритм работал, нужно иметь на каждом шаге (см. Определение ). Если в какой-то момент это предположение не срабатывает, нужно заменить n -ю строку другой строкой под ней, прежде чем продолжить. Вот почему LU-разложение в целом выглядит так . ${\ Displaystyle а_ {п, п} ^ {(п-1)} \ neq 0}$${\ Displaystyle \ ell _ {я, п}}$${\ displaystyle P ^ {- 1} A = LU}$

Обратите внимание, что разложение, полученное с помощью этой процедуры, является разложением Дулитла : главная диагональ L состоит только из 1 с. Если продолжить, удалив элементы выше главной диагонали, добавив несколько столбцов (вместо удаления элементов ниже диагонали путем добавления нескольких строк ), мы получим разложение Кроута , где главная диагональ U составляет 1 с. .

Другой (эквивалент) способ получения разложения Crout данной матрицы А , чтобы получить разложение Дулитло из транспонированной A . В самом деле, если это LU-разложение, полученное с помощью алгоритма, представленного в этом разделе, то, взяв и , мы получим разложение Краута. ${\ textstyle A ^ {\textf {T}} = L_ {0} U_ {0}}$${\ textstyle L = U_ {0} ^ {\textf {T}}}$${\ textstyle U = L_ {0} ^ {\textf {T}}}$${\ displaystyle A = LU}$

Через рекурсию

Cormen et al. описать рекурсивный алгоритм разложения LUP.

Для данной матрицы A пусть P ₁ - матрица перестановок такая, что

${\ displaystyle P_ {1} A = \ left ({\ begin {array} {c | ccc} a && w ^ {\textf {T}} & \\\ hline &&& \\ v && A '& \\ &&& \ end {array }}\Правильно)}$,

где , если в первом столбце A есть ненулевая запись ; или в противном случае возьмем P ₁ в качестве единичной матрицы. Теперь пусть , если ; или иначе. У нас есть ${\ textstyle а \ neq 0}$${\ textstyle c = 1 / a}$${\ textstyle а \ neq 0}$${\ textstyle c = 0}$

${\ displaystyle P_ {1} A = \ left ({\ begin {array} {c | ccc} 1 && 0 & \\\ hline &&& \\ cv && I_ {n-1} & \\ &&& \ end {array}} \ right) \ left ({\ begin {array} {c | c} a & w ^ {\textf {T}} \\\ hline & \\ 0 & A'-cvw ^ {\textf {T}} \\ & \ end {array} }\Правильно).}$

Теперь мы можем рекурсивно найти разложение LUP . Пусть . Следовательно ${\ textstyle P '\ left (A'-cvw ^ {\textf {T}} \ right) = L'U'}$${\ textstyle v '= P'v}$

${\ displaystyle \ left ({\ begin {array} {c | ccc} 1 && 0 & \\\ hline &&& \\ 0 && P '& \\ &&& \ end {array}} \ right) P_ {1} A = \ left ({ \ begin {array} {c | ccc} 1 && 0 & \\\ hline &&& \\ cv '&& L' & \\ &&& \ end {array}} \ right) \ left ({\ begin {array} {c | ccc} a && w ^ {\textf {T}} & \\\ hline &&& \\ 0 && U '& \\ &&& \ end {array}} \ right),}$

которая представляет собой разложение LUP из A .

Рандомизированный алгоритм

Можно найти приближение низкого ранга к разложению LU, используя рандомизированный алгоритм. Учитывая входную матрицу и желаемый низкий ранг , рандомизированный LU возвращает матрицы перестановок и нижние / верхние трапециевидные матрицы размера и, соответственно, так что с высокой вероятностью , где - константа, которая зависит от параметров алгоритма и является -й сингулярное значение входной матрицы . ${\ textstyle A}$${\ textstyle k}$${\ textstyle P, Q}$${\ textstyle L, U}$${\ textstyle m \ times k}$${\ textstyle k \ times n}$${\ textstyle \ left \ | PAQ-LU \ right \ | _ {2} \ leq C \ sigma _ {k + 1}}$${\ textstyle C}$${\ textstyle \ sigma _ {к + 1}}$${\ textstyle (к + 1)}$${\ textstyle A}$

Теоретическая сложность

Если две матрицы порядка n могут быть умножены за время M ( n ), где M ( n ) ≥ n ^a для некоторого n > 2, то разложение LU может быть вычислено за время O ( M ( n )). Это означает, например, что существует алгоритм O ( n ^2.376 ), основанный на алгоритме Копперсмита – Винограда .

Разложение по разреженной матрице

Были разработаны специальные алгоритмы факторизации больших разреженных матриц . Эти алгоритмы пытаются найти редкие факторы L и U . В идеале стоимость вычислений определяется количеством ненулевых элементов, а не размером матрицы.

Эти алгоритмы используют свободу обмена строками и столбцами для минимизации заполнения (записи, которые изменяются с начального нуля на ненулевое значение во время выполнения алгоритма).

Общая трактовка порядков, которые минимизируют заполнение, может быть решена с помощью теории графов .

Приложения

Решение линейных уравнений

Для системы линейных уравнений в матричной форме

${\ Displaystyle А \ mathbf {x} = \ mathbf {b},}$

мы хотим решить уравнение относительно x , учитывая A и b . Предположим, мы уже получили LUP-разложение A такое, что , so . ${\ textstyle PA = LU}$${\ textstyle LU \ mathbf {x} = P \ mathbf {b}}$

В этом случае решение осуществляется в два логических шага:

1. Сначала мы решаем уравнение относительно y .${\ textstyle L \ mathbf {y} = P \ mathbf {b}}$
2. Во-вторых, мы решаем уравнение относительно x .${\ textstyle U \ mathbf {x} = \ mathbf {y}}$

В обоих случаях мы имеем дело с треугольными матрицами ( L и U ), которые можно решить напрямую путем прямой и обратной подстановки без использования процесса исключения Гаусса (однако нам действительно нужен этот процесс или его эквивалент для вычисления самого LU- разложения).

Вышеупомянутую процедуру можно многократно применять для решения уравнения несколько раз для различных b . В этом случае быстрее (и удобнее) выполнить LU-разложение матрицы A один раз, а затем решить треугольные матрицы для различных b , вместо того, чтобы каждый раз использовать исключение Гаусса. Можно подумать, что матрицы L и U «закодировали» процесс исключения Гаусса.

Стоимость решения системы линейных уравнений приблизительно равна операциям с плавающей запятой, если матрица имеет размер . Это делает его вдвое быстрее, чем алгоритмы, основанные на QR-разложении , которое требует операций с плавающей запятой при использовании отражений Хаусхолдера . По этой причине обычно предпочтительнее разложение LU. ${\ textstyle {\ frac {2} {3}} п ^ {3}}$${\ textstyle A}$${\ textstyle n}$${\ textstyle {\ frac {4} {3}} п ^ {3}}$

Инвертирование матрицы

При решении системы уравнений, б , как правило , рассматривается как вектор с длиной , равной высоте матрицы A . В матричной инверсии , однако, вместо вектора Ь , мы имеем матрицу B , где B является в н матрице с размерностью р матрицей, так что мы пытаемся найти матрицу X (также н матрицы с размерностью р матрица):

${\ displaystyle AX = LUX = B.}$

Мы можем использовать тот же алгоритм , представленный ранее решить для каждого столбца матрицы X . Теперь предположим, что B - единичная матрица размера n . Он будет следовать , что результат X должен быть обратным A .

Вычисление определителя

Учитывая LUP-разложение квадратной матрицы A , определитель A может быть вычислен напрямую как ${\ textstyle A = P ^ {- 1} LU}$

${\ displaystyle \ det (A) = \ det \ left (P ^ {- 1} \ right) \ det (L) \ det (U) = (- 1) ^ {S} \ left (\ prod _ {i = 1} ^ {n} l_ {ii} \ right) \ left (\ prod _ {i = 1} ^ {n} u_ {ii} \ right).}$

Второе уравнение следует из того факта, что определитель треугольной матрицы является просто произведением ее диагональных элементов и что определитель матрицы перестановок равен (−1) ^S, где S - количество замен строк в разложении. .

В случае разложения LU с полным поворотом, также равняется правой части приведенного выше уравнения, если мы позволим S быть общим количеством обменов строк и столбцов. ${\ textstyle \ det (A)}$

Тот же метод легко применяется к разложению LU, устанавливая P равным единичной матрице.

Примеры кода

Пример кода C

/* INPUT: A - array of pointers to rows of a square matrix having dimension N
 *        Tol - small tolerance number to detect failure when the matrix is near degenerate
 * OUTPUT: Matrix A is changed, it contains a copy of both matrices L-E and U as A=(L-E)+U such that P*A=L*U.
 *        The permutation matrix is not stored as a matrix, but in an integer vector P of size N+1 
 *        containing column indexes where the permutation matrix has "1". The last element P[N]=S+N, 
 *        where S is the number of row exchanges needed for determinant computation, det(P)=(-1)^S    
 */
int LUPDecompose(double **A, int N, double Tol, int *P) {

    int i, j, k, imax; 
    double maxA, *ptr, absA;

    for (i = 0; i <= N; i++)
        P[i] = i; //Unit permutation matrix, P[N] initialized with N

    for (i = 0; i < N; i++) {
        maxA = 0.0;
        imax = i;

        for (k = i; k < N; k++)
            if ((absA = fabs(A[k][i])) > maxA) { 
                maxA = absA;
                imax = k;
            }

        if (maxA < Tol) return 0; //failure, matrix is degenerate

        if (imax != i) {
            //pivoting P
            j = P[i];
            P[i] = P[imax];
            P[imax] = j;

            //pivoting rows of A
            ptr = A[i];
            A[i] = A[imax];
            A[imax] = ptr;

            //counting pivots starting from N (for determinant)
            P[N]++;
        }

        for (j = i + 1; j < N; j++) {
            A[j][i] /= A[i][i];

            for (k = i + 1; k < N; k++)
                A[j][k] -= A[j][i] * A[i][k];
        }
    }

    return 1;  //decomposition done 
}

/* INPUT: A,P filled in LUPDecompose; b - rhs vector; N - dimension
 * OUTPUT: x - solution vector of A*x=b
 */
void LUPSolve(double **A, int *P, double *b, int N, double *x) {

    for (int i = 0; i < N; i++) {
        x[i] = b[P[i]];

        for (int k = 0; k < i; k++)
            x[i] -= A[i][k] * x[k];
    }

    for (int i = N - 1; i >= 0; i--) {
        for (int k = i + 1; k < N; k++)
            x[i] -= A[i][k] * x[k];

        x[i] /= A[i][i];
    }
}

/* INPUT: A,P filled in LUPDecompose; N - dimension
 * OUTPUT: IA is the inverse of the initial matrix
 */
void LUPInvert(double **A, int *P, int N, double **IA) {
  
    for (int j = 0; j < N; j++) {
        for (int i = 0; i < N; i++) {
            IA[i][j] = P[i] == j ? 1.0 : 0.0;

            for (int k = 0; k < i; k++)
                IA[i][j] -= A[i][k] * IA[k][j];
        }

        for (int i = N - 1; i >= 0; i--) {
            for (int k = i + 1; k < N; k++)
                IA[i][j] -= A[i][k] * IA[k][j];

            IA[i][j] /= A[i][i];
        }
    }
}

/* INPUT: A,P filled in LUPDecompose; N - dimension. 
 * OUTPUT: Function returns the determinant of the initial matrix
 */
double LUPDeterminant(double **A, int *P, int N) {

    double det = A[0][0];

    for (int i = 1; i < N; i++)
        det *= A[i][i];

    return (P[N] - N) % 2 == 0 ? det : -det;
}

Пример кода C #

public class SystemOfLinearEquations
{
    public double[] SolveUsingLU(double[,] matrix, double[] rightPart, int n)
    {
        // decomposition of matrix
        double[,] lu = new double[n, n];
        double sum = 0;
        for (int i = 0; i < n; i++)
        {
            for (int j = i; j < n; j++)
            {
                sum = 0;
                for (int k = 0; k < i; k++)
                    sum += lu[i, k] * lu[k, j];
                lu[i, j] = matrix[i, j] - sum;
            }
            for (int j = i + 1; j < n; j++)
            {
                sum = 0;
                for (int k = 0; k < i; k++)
                    sum += lu[j, k] * lu[k, i];
                lu[j, i] = (1 / lu[i, i]) * (matrix[j, i] - sum);
            }
        }

        // lu = L+U-I
        // find solution of Ly = b
        double[] y = new double[n];
        for (int i = 0; i < n; i++)
        {
            sum = 0;
            for (int k = 0; k < i; k++)
                sum += lu[i, k] * y[k];
            y[i] = rightPart[i] - sum;
        }
        // find solution of Ux = y
        double[] x = new double[n];
        for (int i = n - 1; i >= 0; i--)
        {
            sum = 0;
            for (int k = i + 1; k < n; k++)
                sum += lu[i, k] * x[k];
            x[i] = (1 / lu[i, i]) * (y[i] - sum);
        }
        return x;
    }
}

Пример кода MATLAB

function LU = LUDecompDoolittle(A)
    n = length(A);
    LU = A;
    % decomposition of matrix, Doolittle’s Method
    for i = 1:1:n
        for j = 1:(i - 1)
            LU(i,j) = (LU(i,j) - LU(i,1:(j - 1))*LU(1:(j - 1),j)) / LU(j,j);
        end
        j = i:n;
        LU(i,j) = LU(i,j) - LU(i,1:(i - 1))*LU(1:(i - 1),j);
    end
    %LU = L+U-I
end

function x = SolveLinearSystem(LU, B)
    n = length(LU);
    y = zeros(size(B));
    % find solution of Ly = B
    for i = 1:n
        y(i,:) = B(i,:) - L(i,1:i)*y(1:i,:);
    end
    % find solution of Ux = y
    x = zeros(size(B));
    for i = n:(-1):1
        x(i,:) = (y(i,:) - U(i,(i + 1):n)*x((i + 1):n,:))/U(i, i);
    end    
end

A = [ 4 3 3; 6 3 3; 3 4 3 ]
LU = LUDecompDoolittle(A)
B = [ 1 2 3; 4 5 6; 7 8 9; 10 11 12 ]'
x = SolveLinearSystem(LU, B)
A * x

Смотрите также

• Блочная декомпозиция LU
• Разложение Брюа
• Разложение Холецкого
• Неполная факторизация LU
• Уменьшение LU
• Разложение матрицы
• QR-разложение

Примечания

использованная литература

• Связка, Джеймс Р .; Hopcroft, Джон (1974), "Треугольная факторизационная и инверсия быстрого умножения матриц", Математика вычислений , 28 (125): 231-236, DOI : 10,2307 / 2005828 , ISSN  0025-5718 , JSTOR  2005828.
• Кормен, Томас Х .; Лейзерсон, Чарльз Э .; Ривест, Рональд Л .; Стейн, Клиффорд (2001), Введение в алгоритмы , MIT Press и McGraw-Hill, ISBN 978-0-262-03293-3.
• Голуб, Джин Х .; Ван Лоан, Чарльз Ф. (1996), Матричные вычисления (3-е изд.), Балтимор: Джонс Хопкинс, ISBN 978-0-8018-5414-9.
• Хорн, Роджер А .; Джонсон, Чарльз Р. (1985), Матричный анализ , Cambridge University Press, ISBN 978-0-521-38632-6. См. Раздел 3.5. N  - 1
• Хаусхолдер, Алстон С. (1975), Теория матриц в численном анализе , Нью-Йорк: Dover Publications , MR  0378371.
• Окунев, Павел; Джонсон, Чарльз Р. (1997), Необходимые и достаточные условия для существования LU-факторизации произвольной матрицы , arXiv : math.NA/0506382.
• Пул, Дэвид (2006), Линейная алгебра: современное введение (2-е изд.), Канада: Томсон Брукс / Коул, ISBN 978-0-534-99845-5.
• Нажмите, WH; Теукольский, С.А.; Феттерлинг, штат Вашингтон; Фланнери, Б.П. (2007), «Раздел 2.3» , Численные рецепты: Искусство научных вычислений (3-е изд.), Нью-Йорк: Cambridge University Press, ISBN. 978-0-521-88068-8.
• Trefethen, Lloyd N .; Бау, Дэвид (1997), Численная линейная алгебра , Филадельфия: Общество промышленной и прикладной математики, ISBN 978-0-89871-361-9.

внешние ссылки

использованная литература

• Разложение LU на MathWorld .
• Разложение LU на Math-Linux .
• LU-декомпозиция в Институте целостных численных методов
• Факторизация матрицы LU . Ссылка на MATLAB.

Компьютерный код

• LAPACK - это набор подпрограмм FORTRAN для решения задач плотной линейной алгебры.
• ALGLIB включает частичный перенос LAPACK на C ++, C #, Delphi и т. Д.
• Код C ++ , профессор Дж. Лумис, Дейтонский университет
• Код C , Библиотека исходных текстов математики
• LU в X10

Интернет-ресурсы

• WebApp описательно решает системы линейных уравнений с LU Decomposition
• Матричный калькулятор , bluebit.gr
• LU Decomposition Tool , uni-bonn.de
• LU Decomposition , Эд Пегг-младший , The Wolfram Demonstrations Project , 2007.