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. То есть мы можем записать A как

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

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

${\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 Crout

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

Другой (эквивалентный) способ получения разложения Краута заданной матрицы 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 равным единичной матрице.

Примеры кода

Пример кода на языке С

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

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

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

функция  LU = LUDecompDoolittle ( A ) n = длина ( A ); LU = A ; % разложение матрицы, метод Дулиттла для i = 1 : 1 : n для j = 1 :( i - 1 ) LU ( i , j ) = ( LU ( i , j ) - LU ( i , 1 :( j - 1 )) * LU ( 1 :( j - 1 ), j )) / LU ( j , j ); конец j = i : n ; LU ( i , j ) = LU ( i , j ) - LU ( i , 1 :( i - 1 )) * LU ( 1 :( i - 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 ,:); конец % найти решение Ux = y x = нули ( размер ( B )); для i = n :( - 1 ): 1 x ( i ,:) = ( y ( i ,:) - LU ( i ,( i + 1 ): n ) * x (( i + 1 ): n ,:)) / LU ( i , i ); конец конец                                       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 = РешитьЛинейнуюСистему ( LU , B ) A * x                                  

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

• Блок LU-декомпозиции
• разложение Брюа
• Разложение Холецкого
• Разложение матрицы Кроута
• Неполная LU-факторизация
• Сокращение LU
• Матричное разложение
• QR-разложение

Примечания

1. Шварценберг-Черни, А. (1995). «О матричной факторизации и эффективном решении методом наименьших квадратов». Серия приложений к астрономии и астрофизике . 110 : 405. Bibcode : 1995A&AS..110..405S.
2. Дуайер (1951).
3. Okunev & Johnson (1997), Следствие 3.
4. Трефетен и Бау (1997), с. 166.
5. Трефетен и Бау (1997), с. 161.
6. Lay, Lay & McDonald (2021), стр. 133, 2.5: Факторизация матриц.
7. Риготти (2001), ведущий минор.
8. Horn & Johnson (1985), Следствие 3.5.5
9. Хорн и Джонсон (1985), Теорема 3.5.2.
10. Nhiayi, Ly; Phan-Yamada, Tuyetdong (2021). «Исследование возможного разложения LU». Североамериканский журнал GeoGebra . 9 (1).
11. Окунев и Джонсон (1997).
12. Домохозяин (1975).
13. Голуб и Ван Лоан (1996), стр. 112, 119.
14. Кормен и др. (2009), стр. 819, 28.1: Решение систем линейных уравнений.
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), «Треугольная факторизация и инверсия с помощью быстрого умножения матриц», Математика вычислений , 28 (125): 231–236, doi : 10.2307/2005828 , hdl : 1813/6003 , ISSN  0025-5718, JSTOR  2005828.
• Кормен, Томас Х .; Лейзерсон, Чарльз Э .; Ривест, Рональд Л .; Штейн, Клиффорд (2009), Введение в алгоритмы (3-е изд.), MIT Press и McGraw-Hill, ISBN 978-0-262-03293-3.
• Дуайер, Пол С. (1951), Линейные вычисления , Нью-Йорк: Wiley.
• Голуб, Джин Х.; Ван Лоан, Чарльз Ф. (1996), Матричные вычисления (3-е изд.), Балтимор: Джонс Хопкинс, ISBN 978-0-8018-5414-9.
• Хорн, Роджер А.; Джонсон, Чарльз Р. (1985), Матричный анализ , Cambridge University Press, ISBN 978-0-521-38632-6. См. раздел 3.5. N  − 1
• Хаусхолдер, Олстон С. (1975), Теория матриц в численном анализе , Нью-Йорк: Dover Publications , MR  0378371.
• Лэй, Дэвид К.; Лэй, Стивен Р.; Макдональд, Джуди Дж. (2021), Линейная алгебра и ее приложения (шестое изд.), Pearson, ISBN 978-0-13-585125-8.
• Окунев, Павел; Джонсон, Чарльз Р. (1997), Необходимые и достаточные условия существования LU-факторизации произвольной матрицы , arXiv : math.NA/0506382.
• Пул, Дэвид (2006), Линейная алгебра: Современное введение (2-е изд.), Канада: Thomson Brooks/Cole, ISBN 978-0-534-99845-5.
• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), «Раздел 2.3», Numerical Recipes: The Art of Scientific Computing (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, Библиотека исходного кода математики
• Код Rust
• ЛУ в X10

Онлайн ресурсы

• Веб-приложение, описательно решающее системы линейных уравнений с помощью LU-разложения
• Калькулятор матриц с шагами, включая LU-разложение,
• Инструмент декомпозиции LU, uni-bonn.de
• Разложение LU Эда Пегга-младшего , The Wolfram Demonstrations Project , 2007.