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

Тип матричной факторизации

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

Определения

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

Пусть A — квадратная матрица. Факторизация LU относится к факторизации A с правильным порядком строк и/или столбцов или перестановками на два фактора — нижнюю треугольную матрицу 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 — снова нижняя и верхняя треугольные матрицы, а P — матрица перестановок , которая при умножении слева на A переупорядочивает строки A. Оказывается, что все квадратные матрицы могут быть факторизованы в этой форме ^[3] , и на практике факторизация численно устойчива. ^[4] Это делает LUP-разложение полезным на практике методом.

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

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

${\displaystyle PAQ=LU,}$

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

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

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

${\displaystyle A=LDU,}$

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

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

Выше мы требовали, чтобы A была квадратной матрицей, но все эти разложения можно обобщить и на прямоугольные матрицы. ^[6] В этом случае 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{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 была единичной треугольной матрицей (т. е. приравняла все элементы ее главной диагонали к единицам). Тогда система уравнений имеет следующее решение:

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

Подстановка этих значений в приведенное выше 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\end{bmatrix}}.}$

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

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

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

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

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

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

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

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

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

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

Это разложение называется разложением Холецкого . Если положительно определено, то разложение Холецкого существует и единственно. Более того, вычисление разложения Холецкого более эффективно и численно более стабильно , чем вычисление некоторых других LU-разложений. ${\displaystyle A}$

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

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

Алгоритмы

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

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

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

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

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

Обозначения

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

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

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

Процедура

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

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

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

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

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

Пример

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

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

${\displaystyle A}$ ${\displaystyle P}$

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

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

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

${\displaystyle A^{(1)}}$

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

${\displaystyle P}$

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

${\displaystyle row_{3}=row_{3}-(\ell _{3,2})\cdot row_{2}}$ ${\textstyle \ell _{3,2}={\frac {5}{6}}}$ ${\displaystyle A}$ ${\displaystyle A}$ ${\displaystyle U}$ ${\displaystyle P}$

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

${\displaystyle L}$

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

${\displaystyle PA=LU}$

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

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

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

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

$${\displaystyle 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}}$ ${\displaystyle A=LA^{(N-1)}=LU}$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Приложения

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

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

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

мы хотим решить уравнение для x , учитывая A и b . Предположим, мы уже получили LUP-разложение A такое, что , поэтому . ${\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-разложение. ^[17] ${\textstyle {\frac {2}{3}}n^{3}}$ ${\textstyle A}$ ${\textstyle n}$ ${\textstyle {\frac {4}{3}}n^{3}}$

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

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

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

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

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

Учитывая 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

/* ВХОД: 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            для ( я знак равно 0 ; я < N ; я ++ ) { maxA = 0,0 ; имакс = я ;               for ( k = i ; k < N ; k ++ ) if (( absA = fabs ( A [ k ][ i ])) > maxA ) { maxA = absA ; имакс = к ; }                       если ( maxA < Tol ) вернуть 0 ; //неудача, матрица вырождена       if ( imax != i ) { // поворот P j = P [ i ]; P [ я ] = P [ imax ]; п [ имакс ] знак равно j ;               // поворот строк A ptr = A [ i ]; А [ я ] = А [ imax ]; А [ имакс ] = ПТР ;          //подсчет пивотов, начиная с N (для определителя) P [ N ] ++ ; }   для ( 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 ; //разложение завершено }  /* ВВОД: A,P заполнены LUPDecompose; б - правый вектор; N - размерность * ВЫХОД: x - вектор решения A*x=b */ void LUPSolve ( double ** A , int * P , double * b , int N , double * x ) {            для ( int я знак равно 0 ; я < 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 я знак равно N - 1 ; я >= 0 ; я -- ) { for ( int k знак равно я + 1 ; k < N ; k ++ ) IA [ я ][ j ] -= A [ я ][ к ] * ИА [ к ][ j ];                            IA [ я ][ j ] /= А [ я ][ я ]; } } }    /* ВВОД: A,P заполнены LUPDecompose; Н - размерность. * ВЫХОД: Функция возвращает определитель исходной матрицы */ double LUPDeterminant ( double ** A , int * P , int N ) {        двойной дет = А [ 0 ][ 0 ];    for ( int я знак равно 1 ; я < N ; я ++ ) det *= A [ я ][ я ];            возврат ( P [ N ] - N ) % 2 == 0 ? дет : - дет ; }           

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

public class SystemOfLinearEquations { public double [] SolveUsingLU ( двойная [,] матрица , двойная [] 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 ); конец %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-разложение

Примечания

1. Шварценберг-Черни, А. (1995). «О матричной факторизации и эффективном решении методом наименьших квадратов». Серия дополнений по астрономии и астрофизике . 110 : 405. Бибкод : 1995A&AS..110..405S.
2. Дуайер, Пол С. (1951). Линейные вычисления . Нью-Йорк: Уайли.
3. Okunev & Johnson (1997), следствие 3.
4. Трефетен и Бау (1997), с. 166.
5. Трефетен и Бау (1997), с. 161.
6. Лэй, Дэвид К. (2016). Линейная алгебра и ее приложения. Стивен Р. Лэй, Джудит Макдональд (Пятое изд.). Харлоу. п. 142. ИСБН 978-1-292-09223-2. ОСЛК  920463015.{{cite book}}: CS1 maint: location missing publisher (link)
7. Риготти (2001), Второстепенный ведущий принцип
8. Horn & Johnson (1985), следствие 3.5.5
9. Хорн и Джонсон (1985), Теорема 3.5.2.
10. Нхиайи, Ли; Фан-Ямада, Туетдонг (2021). «Изучение возможного разложения LU». Североамериканский журнал GeoGebra . 9 (1).
11. Окунев и Джонсон (1997)
12. Домохозяин (1975)
13. Голуб и Ван Лоан (1996), с. 112, 119.
14. Кормен, Томас Х .; Лейзерсон, Чарльз Э .; Ривест, Рональд Л .; Штейн, Клиффорд (2001). Введение в алгоритмы . MIT Press и McGraw-Hill. ISBN 978-0-262-03293-3.
15. Шабат, Гил; Шмуэли, Янив; Айзенбуд, Ярив; Авербух, Амир (2016). «Рандомизированное LU-разложение». Прикладной и вычислительный гармонический анализ . 44 (2): 246–272. arXiv : 1310.7202 . дои :10.1016/j.acha.2016.04.006. S2CID  1900701.
16. Банч и Хопкрофт (1974)
17. Трефетен и Бау (1997), с. 152.
18. Голуб и Ван Лоан (1996), с. 121

Рекомендации

• Банч, Джеймс Р.; Хопкрофт, Джон (1974), «Треугольная факторизация и инверсия путем быстрого умножения матриц», Mathematics of Computation , 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), Матричный анализ , Издательство Кембриджского университета, 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.