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 — снова нижняя и верхняя треугольные матрицы, а P — матрица перестановок при умножении слева на A переупорядочивает строки A. , которая Оказывается, все квадратные матрицы можно факторизовать в таком виде: и факторизация на практике численно устойчива. Это делает 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

/* ВХОД: A - массив указателей на строки квадратной матрицы размером N 
 * Tol - небольшое число допуска для обнаружения сбоя, когда матрица близка к вырождению 
 * ВЫХОД: Матрица A изменяется, она содержит копии обеих матриц LE и U как A=(LE)+U такой, что P*A=L*U. 
  * Матрица перестановок хранится не в виде матрицы, а в виде целочисленного вектора P размера N+1 
 *, содержащего индексы столбцов, где матрица перестановок имеет значение «1».   Последний элемент P[N]=S+N, 
 * где S — количество перестановок строк, необходимых для вычисления определителя, det(P)=(-1)^S 
 */ 
 int   LUPDecompose  (  double   **  A  ,   int   N  ,   двойной   Тол  ,   int   *  P  )   { 

     int   я  ,   j  ,   k  ,   imax  ;  
      двойной   maxA  ,   *  ptr  ,   absA  ; 

      для   (  я   знак равно   0  ;   я   <=   N  ;   я  ++  ) 
         п  [  я  ]   знак равно   я  ;    //Матрица перестановок единиц, P[N] инициализируется N 

     for   (  i   =   0  ;   i   <   N  ;   i  ++  )   { 
         maxA   =   0.0  ; 
          имакс   =   я  ; 

          for   (  k   =   i  ;   k   <   N  ;   k  ++  ) 
             if   ((  absA   =   fabs  (  A  [  k  ][  i  ]))   >   maxA  )   {  
                 maxA   =   absA  ; 
                  имакс   =   к  ; 
              } 

         if   (  maxA   <   Tol  )   возвращает   0  ;    // ошибка, матрица вырождена 

         if   (  imax   !=   i  )   { 
             // поворот P 
             j   =   P  [  i  ]; 
              п  [  я  ]   =   п  [  имакс  ]; 
              п  [  имакс  ]   знак равно   j  ; 

              // поворот строк A 
             ptr   =   A  [  i  ]; 
              А  [  я  ]   =   А  [  imax  ]; 
              А  [  имакс  ]   =   ПТР  ; 

              //считаем пивоты, начиная с N (для определителя) 
            П  [  Н  ]  ++  ; 
          } 

         for   (  j   знак равно   я   +   1  ;   j   <   N  ;   j  ++  )   { 
             A  [  j  ] [  я  ]   / =   A  [  я  ] [  я  ]; 

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

     вернуть   1  ;     //декомпозиция завершена 
 } 

 /* INPUT: A,P заполнены LUPDecompose;   б - правый вектор;   N - размерность 
 * ВЫХОД: x - вектор решения A*x=b 
 */ 
 void   LUPSolve  (  double   **  A  ,   int   *  P  ,   double   *  b  ,   int   N  ,   double   *  x  )   { 

     for   (  int   i   =   0  ;   i   <   N  ;   я  ++  )   { 
         Икс  [  я  ]   знак равно   б  [  п  [  я  ]]; 

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

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

          Икс  [  я  ]   / =   А  [  я  ] [  я  ]; 
      } 
 } 

 /* ВВОД: A,P заполнены LUPDecompose;   N — размерность 
 * ВЫХОД: IA — обратная исходной матрице 
 */ 
 void   LUPInvert  (  double   **  A  ,   int   *  P  ,   int   N  ,   double   **  IA  )   { 
  
     for   (  int   j   =   0  ;   j   <   N  ;   j  + +  )   { 
         for   (  int   я   знак равно   0  ;   я   <   N  ;   я  ++  )   { 
             IA  [  я  ][  j  ]   знак равно   п  [  я  ]   ==   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  ;   я   >=   0  ;   я  --  )   { 
             for   (  int   k   =   i   +   1  ;   k   <   N  ;   k  ++  ) 
                 IA  [  i  ][  j  ]   -=   A  [  i  ] [  к  ]   *   ИА  [  к  ][  j  ]; 

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

 /* INPUT: A,P заполнены LUPDecompose;   Н - размерность. 
  * ВЫХОД: Функция возвращает определитель исходной матрицы 
 */ 
 double   LUPDeterminant  (  double   **  A  ,   int   *  P  ,   int   N  )   { 

     double   det   =   A  [  0  ][  0  ]; 

      for   (  int   я   знак равно   1  ;   я   <   N  ;   я  ++  ) 
         det   *=   A  [  я  ][  я  ]; 

    возврат   (  P  [  N  ]   -   N  )   %   2   ==   0   ?    что   :   -  это  ; 
  } 

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

public   class   SystemOfLinearEquations 
 { 
     public   double  []   SolveUsingLU  (  double  [,]   матрица  ,   double  []   rightPart  ,   int   n  ) 
     { 
         // разложение матрицы 
         double  [,]   lu   =   new   double  [  n  ,   n  ]; 
          двойная   сумма   =   0  ; 
          for   (  int   я знак   равно   0  ;   я   <   n  ;   я  ++  ) 
         { 
             for   (  int   j   знак равно   я  ;   j   <   n  ;   j  ++  ) 
             { 
                 sum   =   0  ; 
                  for   (  int   k   =   0  ;   k   <   i  ;   k  ++  ) 
                     sum   +=   lu  [  i  ,   k  ]   *   lu  [  k  ,   j  ]; 
                  лу  [  я  ,   j  ]   =   матрица  [  я  ,   j  ]   -   сумма  ; 
              } 
             for   (  int   j   знак равно   я   +   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  ])   *   (  матрица  [  j  ,   i  ]   -   сумма  ); 
              } 
         } 

         // lu = L+UI 
         // находим решение Ly = b 
         double  []   y   =   new   double  [  n  ]; 
          для   (  int   я   знак равно   0  ;   я   <   п  ;   я  ++  ) 
         { 
             сумма   =   0  ; 
              for   (  int   k   =   0  ;   k   <   i  ;   k  ++  ) 
                 sum   +=   lu  [  i  ,   k  ]   *   y  [  k  ]; 
              y  [  я  ]   =   правая часть  [  я  ]   -   сумма  ; 
          } 
         // находим решение Ux = y 
         double  []   x   =   new   double  [  n  ]; 
          для   (  int   я   знак равно   п   -   1  ;   я   >=   0  ;   я  --  ) 
         { 
             сумма   =   0  ; 
              for   (  int   k   =   i   +   1  ;   k   <   n  ;   k  ++  ) 
                 sum   +=   lu  [  i  ,   k  ]   *   x  [  k  ]; 
              Икс  [  я  ]   =   (  1   /   lu  [  я  ,   я  ])   *   (  y  [  я  ]   -   сумма  ); 
          } 
         Вернуть   х  ; 
      } 
 } 

Пример кода MATLAB

функция   LU   =   LUDecompDoolittle  (  A  ) 
     n   =   длина  (  A  ); 
      ЛУ   =   А  ; 
      % разложения матрицы, метод Дулитла 
     для   i   =   1  :  1  :  n 
         для   j   =   1  :(  i   -   1  ) 
             LU  (  i  ,  j  )   =   (  LU  (  i  ,  j  )   -   LU  (  i  ,  1  :(  j   -   1)  ))  *  ЛУ  (  1  :(  j   -   1  ),  j  ))   /   LU  (  j  ,  j  ); 
          конец 
         j   знак равно   я  :  п  ; 
          ЛУ  (  я  ,  j  )   =   ЛУ  (  я  ,  j  )   -   ЛУ  (  я  ,  1  :  ( я   -   1  ))  *  ЛУ  (  1  :  ( я   -   1  ),  j  ); 
      end 
     %LU = L+UI 
 конечная 

 функция   x   =   SolveLinearSystem  (  LU, B  ) 
     n   =   длина  (  LU  ); 
      y   =   нули  (  размер  (  B  )); 
      % найти решение Ly = B 
     для   i   =   1  :  n 
         y  (  i  ,:)   =   B  (  i  ,:)   -   LU  (  i  ,  1  :  i  )  *  y  (  1  :  i  ,:); 
      end 
     % найти решение Ux = y 
     x   =   нули  (  размер  (  B  )); 
      для   i   =   n  :(  -  1  ):  1 
         x  (  i  ,:)   =   (  y  (  i  ,:)   -   LU  (  i  ,(  i   +   1  ):  n  )  *  x  ((  i   +   1  ):  n  ,: ))  /  ЛУ  (  я  ,   я  ); 
      конец     
 конец 

 А   знак равно   [   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.