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

В численном анализе и алгебре линейной нижне-верхнее ( LU ) разложение или факторизация факторизуют матрицу как произведение нижней треугольной матрицы и верхней треугольной матрицы (см. матричное разложение ). Иногда продукт включает в себя матрицу перестановок также . LU-разложение можно рассматривать как матричную форму исключения Гаусса . Компьютеры обычно решают квадратные системы линейных уравнений с использованием LU-разложения, и это также ключевой шаг при инвертировании матрицы или вычислении определителя матрицы. Разложение LU было предложено польским астрономом Тадеушем Банакевичем в 1938 году. Цитируем: «Похоже, что Гаусс и Дулитл применили метод [исключения] только к симметричным уравнениям. Более поздние авторы, например, Эйткен, Банакевич, Дуайер и Краут… подчеркивали использование этого метода или его вариаций в связи с несимметричными проблемами… Банакевич… видел суть… в том, что основная проблема на самом деле одна. матричной факторизации, или «разложения», как он это называл». Его также иногда называют LR- разложением (факторы на левую и правую треугольные матрицы).

Определения

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

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

$A=LU.$

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

${\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}}.$

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

факторизация с частичным

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

$PA=LU,$

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

факторизация с полным

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

$PAQ=LU,$

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

по нижней диагонали-верхней (

Разложение «нижняя диагональ-верхняя» (LDU) — это разложение вида

$A=LDU,$

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

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

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

Пример

Мы факторизуем следующую матрицу 2х2:

${\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-разложение этой простой матрицы — просто решить линейные уравнения путем проверки. Расширение матричного умножения дает

${\begin{aligned}\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{aligned}}$

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

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

Подстановка этих значений в приведенное выше LU-разложение дает

${\begin{bmatrix}4&3\\6&3\end{bmatrix}}={\begin{bmatrix}1&0\\1.5&1\end{bmatrix}}{\begin{bmatrix}4&3\\0&-1.5\end{bmatrix}}.$

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

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

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

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

Вообще любая квадратная матрица $A_{n\times n}$ может иметь одно из следующих значений:

1. уникальная LU-факторизация (как упоминалось выше);
2. бесконечно много LU-факторизаций, если два или более любых первых ( n -1) столбцов линейно зависимы или любой из первых ( n -1) столбцов равен 0;
3. нет LU-факторизации, если первые ( n −1) столбцы ненулевые и линейно независимые и хотя бы один ведущий главный минор равен нулю.

В случае 3 можно аппроксимировать LU-факторизацию, изменив диагональную запись. $a_{jj}$ к $a_{jj}\pm \varepsilon$ чтобы избежать нулевого ведущего главного минора.

положительно

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

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

Это разложение называется разложением Холецкого . Если $A$ положительно определена, то разложение Холецкого существует и единственно. Более того, вычисление разложения Холецкого более эффективно и численно более стабильно, чем вычисление некоторых других 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 (i-1)}$-я главная подматрица. Вычисление определителей требует больших вычислительных затрат , поэтому эта явная формула на практике не используется.

Использование

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

Общее объяснение

Обозначения

Учитывая размера N × N матрицу $A=(a_{i,j})_{1\leq i,j\leq N}$, определять $A^{(0)}$ как оригинальная, немодифицированная версия матрицы $A$. Надстрочный индекс в скобках (например, $(0)$) матрицы $A$ это версия матрицы. Матрица $A^{(n)}$ это $A$ матрица, в которой элементы ниже главной диагонали уже исключены до 0 методом исключения Гаусса для первого $n$ столбцы.

Ниже приведена матрица, которая поможет нам запомнить обозначения (где каждый $*$ представляет любое действительное число в матрице):

$A^{(n-1)}={\begin{pmatrix}*&&&\cdots &&&*\\0&\ddots &&&&\\&\ddots &*&&&\\\vdots &&0&a_{n,n}^{(n-1)}&&&\vdots \\&&\vdots &a_{i,n}^{(n-1)}&*\\&&&\vdots &\vdots &\ddots \\0&\cdots &0&a_{i,n}^{(n-1)}&*&\cdots &*\end{pmatrix}}$

Процедура

В ходе этого процесса мы постепенно изменяем матрицу $A$ используя операции со строками, пока она не станет матрицей $U$ в котором все элементы ниже главной диагонали равны нулю. При этом мы одновременно создадим две отдельные матрицы $P$ и $L$, такой, что $PA=LU$.

Мы определяем окончательную матрицу перестановок $P$ как единичную матрицу, в которой все одни и те же строки поменяны местами в том же порядке, что и $A$ матрица, пока она преобразуется в матрицу $U$. Для нашей матрицы $A^{(n-1)}$, мы можем начать с замены строк, чтобы обеспечить желаемые условия для n-го столбца. Например, мы можем поменять местами строки, чтобы выполнить частичный поворот, или сделать это, чтобы установить элемент поворота. $a_{n,n}$ на главной диагонали к ненулевому числу, чтобы мы могли завершить исключение Гаусса.

Для нашей матрицы $A^{(n-1)}$, мы хотим установить каждый элемент ниже $a_{n,n}^{(n-1)}$ до нуля (где $a_{n,n}^{(n-1)}$ – элемент n-го столбца главной диагонали). Ниже мы обозначим каждый элемент $a_{n,n}^{(n-1)}$ как $a_{i,n}^{(n-1)}$ (где $i=n+1,\dotsc ,N$). Чтобы установить $a_{i,n}^{(n-1)}$ равным нулю, мы устанавливаем $row_{i}=row_{i}-(\ell _{i,n})\cdot row_{n}$ для каждой строки $i$. Для этой операции ${\textstyle \ell _{i,n}:={a_{i,n}^{(n-1)}}/{a_{n,n}^{(n-1)}}}$. Как только мы выполнили операции со строками для первого $N-1$ столбцов, мы получили верхнетреугольную матрицу $A^{(N-1)}$ который обозначается $U$.

Мы также можем создать нижнюю треугольную матрицу, обозначенную как ${\textstyle L}$, путем непосредственного ввода ранее рассчитанных значений $\ell _{i,n}$ по формуле ниже.

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

Пример

Если нам дана матрица

$$A={\begin{pmatrix}0&5&{\frac {22}{3}}\\4&2&1\\2&7&9\\\end{pmatrix}},$$

мы выберем реализацию частичного поворота и, таким образом, поменяем местами первую и вторую строки, чтобы наша матрица $A$ и первая итерация нашего $P$ матрица соответственно становится

$$A^{(0)}={\begin{pmatrix}4&2&1\\0&5&{\frac {22}{3}}\\2&7&9\\\end{pmatrix}},\quad P^{(0)}={\begin{pmatrix}0&1&0\\1&0&0\\0&0&1\\\end{pmatrix}}.$$

После того, как мы поменяли строки местами, мы можем исключить элементы ниже главной диагонали в первом столбце, выполнив

$${\begin{alignedat}{0}row_{2}=row_{2}-(\ell _{2,1})\cdot row_{1}\\row_{3}=row_{3}-(\ell _{3,1})\cdot row_{1}\end{alignedat}}$$

такой, что,

$${\begin{alignedat}{0}\ell _{2,1}={\frac {0}{4}}=0\\\ell _{3,1}={\frac {2}{4}}=0.5\end{alignedat}}$$

После вычитания этих строк мы получили $A^{(1)}$ матрица

$$A^{(1)}={\begin{pmatrix}4&2&1\\0&5&{\frac {22}{3}}\\0&6&8.5\\\end{pmatrix}}.$$

Поскольку мы реализуем частичный поворот, мы меняем местами вторую и третью строки нашей производной матрицы и текущую версию нашей $P$ матрица соответственно, чтобы получить

$$A^{(1)}={\begin{pmatrix}4&2&1\\0&6&8.5\\0&5&{\frac {22}{3}}\\\end{pmatrix}},\quad P^{(1)}={\begin{pmatrix}0&1&0\\0&0&1\\1&0&0\\\end{pmatrix}}.$$

Теперь мы исключим элементы ниже главной диагонали во втором столбце, выполнив $row_{3}=row_{3}-(\ell _{3,2})\cdot row_{2}$ такой, что ${\textstyle \ell _{3,2}={\frac {5}{6}}}$. Поскольку в нашей текущей итерации ниже главной диагонали не существует ненулевых элементов. $A$ после вычитания этой строки это вычитание строки дает наш окончательный результат $A$ матрица (обозначается как $U$) и финал $P$ матрица:

$$A^{(2)}=A^{(N-1)}=U={\begin{pmatrix}4&2&1\\0&6&8.5\\0&0&0.25\\\end{pmatrix}},\quad P={\begin{pmatrix}0&1&0\\0&0&1\\1&0&0\\\end{pmatrix}}.$$

Также переключив соответствующие строки, мы получим наш окончательный результат. $L$ матрица:

$$L={\begin{pmatrix}1&0&0\\\ell _{3,1}&1&0\\\ell _{2,1}&\ell _{3,2}&1\\\end{pmatrix}}={\begin{pmatrix}1&0&0\\0.5&1&0\\0&{\frac {5}{6}}&1\\\end{pmatrix}}$$

Теперь эти матрицы имеют такое отношение, что $PA=LU$.

Отношения, когда строки не меняются местами

Если мы вообще не меняли местами строки во время этого процесса, мы можем выполнять операции со строками одновременно для каждого столбца. $n$ установив $A^{(n)}:=L_{n}^{-1}A^{(n-1)},$ где $L_{n}^{-1}$ представляет собой единичную матрицу размера N × N , в которой ее n -й столбец заменен транспонированным вектором ${\begin{pmatrix}0&\dotsm &0&1&-\ell _{n+1,n}&\dotsm &-\ell _{N,n}\end{pmatrix}}^{\textsf {T}}.$ Другими словами, нижняя треугольная матрица

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

Выполнение всех операций над строками для первого $N-1$ столбцы с использованием $A^{(n)}:=L_{n}^{-1}A^{(n-1)}$ формула эквивалентна нахождению разложения

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

Обозначим ${\textstyle L=L_{1}\dotsm L_{N-1}}$ так что $A=LA^{(N-1)}=LU$.

Теперь вычислим последовательность $L_{1}\dotsm L_{N-1}$. Мы знаем, что $L_{i}$ имеет следующую формулу.

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

Если есть две нижние треугольные матрицы с единицами на главной диагонали, и ни одна из них не имеет ненулевого элемента ниже главной диагонали в том же столбце, что и другая, то мы можем включить все ненулевые элементы в одно и то же место в произведении. из двух матриц. Например:

$\left({\begin{array}{ccccc}1&0&0&0&0\\77&1&0&0&0\\12&0&1&0&0\\63&0&0&1&0\\7&0&0&0&1\end{array}}\right)\left({\begin{array}{ccccc}1&0&0&0&0\\0&1&0&0&0\\0&22&1&0&0\\0&33&0&1&0\\0&44&0&0&1\end{array}}\right)=\left({\begin{array}{ccccc}1&0&0&0&0\\77&1&0&0&0\\12&22&1&0&0\\63&33&0&1&0\\7&44&0&0&1\end{array}}\right)$

Наконец, умножьте $L_{i}$ вместе и сгенерировать объединенную матрицу, обозначенную как ${\textstyle L}$ (как упоминалось ранее). Использование матрицы ${\textstyle L}$, мы получаем $A=LU.$

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

LU Разложение Краута

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

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

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

Кормен и др. описать рекурсивный алгоритм LUP-разложения.

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

$P_{1}A=\left({\begin{array}{c|ccc}a&&w^{\textsf {T}}&\\\hline &&&\\v&&A'&\\&&&\end{array}}\right)$,

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

$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^{\textsf {T}}\\\hline &\\0&A'-cvw^{\textsf {T}}\\&\end{array}}\right).$

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

$\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^{\textsf {T}}&\\\hline &&&\\0&&U'&\\&&&\end{array}}\right),$

что является LUP-разложением A .

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

Можно найти аппроксимацию низкого ранга для LU-разложения, используя рандомизированный алгоритм. Учитывая входную матрицу ${\textstyle A}$ и желаемый низкий ранг ${\textstyle k}$, рандомизированный LU возвращает матрицы перестановок ${\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 _{k+1}}$ это ${\textstyle (k+1)}$-е сингулярное значение входной матрицы ${\textstyle A}$.

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

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

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

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

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

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

Приложения

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

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

$A\mathbf {x} =\mathbf {b} ,$

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

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

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

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

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

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

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

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

$AX=LUX=B.$

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

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

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

$\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 — количество перестановок строк в разложении.

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

Тот же метод легко применим к 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,:) - LU(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,:) - LU(i,(i + 1):n)*x((i + 1):n,:))/LU(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-факторизация
• Сокращение ЛУ
• Разложение матрицы
• QR-разложение

Ссылки

• Банч, Джеймс Р.; Хопкрофт, Джон (1974), «Треугольная факторизация и инверсия путем быстрого умножения матриц», Mathematics of Computation , 28 (125): 231–236, doi : 10.2307/2005828 , hdl : 1813/6003 , 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), Матричный анализ , Издательство Кембриджского университета, ISBN  978-0-521-38632-6 . См. раздел 3.5. Н - 1
• Хаусхолдер, Алстон С. (1975), Теория матриц в численном анализе , Нью-Йорк: Dover Publications , MR   0378371 .
• Окунев Павел; Джонсон, Чарльз Р. (1997), Необходимые и достаточные условия существования LU-факторизации произвольной матрицы , arXiv : math.NA/0506382 .
• Пул, Дэвид (2006), Линейная алгебра: современное введение (2-е изд.), Канада: Томсон Брукс/Коул, ISBN  978-0-534-99845-5 .
• Пресс, WH; Теукольский, С.А.; Феттерлинг, WT; Фланнери, BP (2007), «Раздел 2.3» , Численные рецепты: искусство научных вычислений (3-е изд.), Нью-Йорк: Cambridge University Press, ISBN  978-0-521-88068-8 .
• Трефетен, Ллойд Н .; Бау, Дэвид (1997), Численная линейная алгебра , Филадельфия: Общество промышленной и прикладной математики, ISBN  978-0-89871-361-9 .
• Риготти, Лука (2001), ECON 2001 - Введение в математические методы, лекция 8

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

Ссылки

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

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

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

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

• WebApp описательное решение систем линейных уравнений с LU-разложением
• Матричный калькулятор с шагами, включая LU-разложение ,
• Инструмент разложения LU , uni-bonn.de
• Разложение LU Эда Пегга-младшего , Демонстрационный проект Вольфрама , 2007.