LU разложение

редактировать

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

Содержание

1 Определения
- 1.1 Факторизация LU с частичным поворотом
- 1.2 Факторизация LU с полным поворотом
- 1.3 LDU разложение
2 Пример
3 Существование и уникальность
- 3.1 Квадратные матрицы
- 3.2 Симметричные положительно определенные матрицы
- 3.3 Общие матрицы
4 Алгоритмы
- 4.1 Замкнутая формула
- 4.2 Использование Гаусса исключение
- 4.3 Через рекурсию
- 4.4 Рандомизированный алгоритм
- 4.5 Теоретическая сложность
- 4.6 Разложение по разреженной матрице
5 Приложения
- 5.1 Решение линейных уравнений
- 5.2 Инвертирование матрицы
- 5.3 Вычисление определителя
- 5.4 Примеры кода C
- 5.5 Примеры кода C #
- 5.6 Примеры кода MATLAB
6 См. Также
7 Примечания
8 Ссылки
9 Внешние ссылки

Определения

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

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

A = L U. {\ displaystyle A = LU.}

{\ displaystyle A = LU.}

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

[a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33] = [l 11 0 0 l 21 l 22 0 l 31 l 32 1 33] [u 11 u 12 u 13 0 u 22 u 23 0 0 u 33]. {\ 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} l_ {11} 0 0 \\ l_ {21} l_ {22} 0 \\ l_ {31} l_ {32} 1_ {33} \ end {bmatrix}} {\ begin {bmatrix} u_ {11} u_ {12} u_ {13} \\ 0 u_ {22} u_ {23} \\ 0 0 u_ {33} \ end {bmatrix}}.}

{\ 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} l_ {11} 0 0 \\ l_ {21} l_ {22} 0 \\ l_ {31} l_ {32} 1_ {33 } \ end {bmatrix}} {\ begin {bmatrix} u_ {11} u_ {12} u_ {13} \\ 0 u_ {22} u_ {23} \\ 0 0 u_ {33} \ end {bmatrix}}.}

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

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

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

PA = LU, {\ displaystyle PA = LU,}

{\ displaystyle PA = LU, }

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

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

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

PAQ = LU, {\ displaystyle PAQ = LU, }

{\ displaystyle PAQ = LU,}

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

разложение LDU

разложение LDU представляет собой разложение формы

A = LDU, {\ displaystyle A = LDU,}

{\ displaystyle A = LDU,}

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

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

Пример

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

[4 3 6 3] = [1 0 l 21 1] [u 11 u 12 0 u 22]. {\ displaystyle {\ begin {bmatrix} 4 3 \\ 6 3 \ end {bmatrix}} = {\ begin {bmatrix} 1 0 \\ l_ {21} 1 \ end {bmatrix}} {\ begin {bmatrix} u_ {11} u_ {12} \\ 0 u_ {22} \ end {bmatrix}}.}

{\ displaystyle {\ begin {bmatrix} 4 3 \\ 6 3 \ end {bmatrix }} = {\ begin {bmatrix} 1 0 \\ l_ {21} 1 \ end {bmatrix}} {\ begin {bmatrix} u_ {11} u_ {12} \\ 0 u_ {22} \ end {bmatrix}}. }

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

1 ⋅ u 11 + 0 ⋅ 0 = 4 1 ⋅ u 12 + 0 ⋅ u 22 = 3 l 21 ⋅ u 11 + 1 ⋅ 0 = 6 l 21 ⋅ u 12 + 1 ⋅ u 22 = 3. {\ Displaystyle {\ begin {align} 1 \ cdot u_ {11} +0 \ cdot 0 = 4 \\ 1 \ cdot u_ {12} +0 \ cdot u_ {22} = 3 \\ l_ { 21} \ cdot u_ {11} +1 \ cdot 0 = 6 \\ l_ {21} \ cdot u_ {12} +1 \ cdot u_ {22} = 3. \ end {align}}}

{\ displaystyle {\ begin {align} 1 \ cdot u_ {11} +0 \ cdot 0 = 4 \\ 1 \ cdot u_ {12} +0 \ cdot u_ {22} = 3 \\ l_ {21} \ cdot u_ {11} +1 \ cdot 0 = 6 \ \ l_ {21} \ cdot u_ {12} +1 \ cdot u_ {22} = 3. \ end {align}}}

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

l 21 = 1,5 u 11 = 4 u 12 = 3 u 22 = - 1,5 {\ displaystyle {\ begin {align} l_ {21} = 1.5 \\ u_ { 11} = 4 \\ u_ {12} = 3 \\ u_ {22} = - 1.5 \ end {align}}}

{\ displaystyle {\ begin {align} l_ {21} = 1.5 \\ u_ {11} = 4 \\ u_ {12} = 3 \\ u_ { 22} = - 1,5 \ конец {выровненный}}}

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

[4 3 6 3] = [1 0 1,5 1] [4 3 0 - 1,5]. {\ 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}}.}

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

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

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

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

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

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

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

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

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

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

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

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

Алгоритмы

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

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

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

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

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

Мы удалите элементы матрицы ниже главной диагонали в n-м столбце матрицы A, добавив к i-й строке этой матрицы n-ю строку, умноженную на

- li, n: = - ai, n (n - 1) an, n (n - 1), i = n + 1,…, N. {\ displaystyle -l_ {i, n}: = - {\ frac {a_ {i, n} ^ {(n-1)}} {a_ {n, n} ^ {(n-1)}}}, \ quad i = n + 1, \ ldots, N.}

{\ displaystyle -l_ {i, n}: = - {\ frac {a_ {i, n} ^ {(n-1)}} {a_ {n, n} ^ {(n-1)}}}, \ quad i = n + 1, \ ldots, N.}

Это можно сделать, умножив A слева на нижнюю треугольную матрицу

L n = (1 ⋱ 1 - ln + 1, n ⋮ ⋱ - l N, n 1), {\ displaystyle L_ {n} = {\ begin {pmatrix} 1 \\ \ ddots \\ 1 \\ - l_ {n + 1, n} \\ \ vdots \ ddots \\ - l_ {N, n} 1 \ end {pmatrix}},}

{\ displaystyle L_ {n} = {\ begin {pmatrix} 1 \\ \ ddots \\ 1 \\ - l_ {n + 1, n} \\ \ vdots \ ddots \\ - l_ {N, n} 1 \ end {pmatrix}},}

то есть единичная матрица N x N, в которой ее n-й столбец заменен вектором $(0 ⋯ 0 1 - ln + 1, n ⋯ - l N, n) T. {\ displaystyle {\ begin {pmatrix} 0 \ cdots 0 1 -l_ {n + 1, n} \ cdots -l_ {N, n} \ end {pmatrix}} ^ {T}.}$ ${\ displaystyle {\ begin {pmatrix} 0 \ cdots 0 1 -l_ {n + 1, n} \ cdots -l_ {N, n} \ end {pmatrix}} ^ {T}.}$

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

A (n): = L n A (n - 1). {\ displaystyle A ^ {(n)}: = L_ {n} A ^ {(n-1)}.}

A ^ {(n)}: = L_n A ^ {(n-1)}.

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

A = L 1 - 1 L 1 A (0) = L 1 - 1 A (1) = L 1 - 1 L 2 - 1 L 2 A (1) = L 1 - 1 L 2 - 1 A (2) = ⋯ = L 1 - 1… LN - 1 - 1 A (N - 1). {\ Displaystyle A = L_ {1} ^ {- 1} L_ {1} A ^ {(0)} = L_ {1} ^ {- 1} A ^ {(1)} = L_ {1} ^ {- 1} L_ {2} ^ {- 1} L_ {2} A ^ {(1)} = L_ {1} ^ {- 1} L_ {2} ^ {- 1} A ^ {(2)} = \ dotsb = L_ {1} ^ {- 1} \ ldots L_ {N-1} ^ {- 1} A ^ {(N-1)}.}

{\ displaystyle A = L_ {1} ^ {- 1} L_ {1} A ^ {(0)} = L_ {1} ^ {- 1} A ^ {(1)} = L_ {1} ^ {- 1} L_ {2} ^ {- 1} L_ {2} A ^ {(1)} = L_ {1} ^ {- 1} L_ {2} ^ {- 1} A ^ {(2)} = \ dotsb = L_ {1} ^ {- 1} \ ldots L_ {N-1} ^ {- 1} A ^ {(N-1)}.}

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

L = (1 l 2, 1 ⋱ ⋱ 1 ⋮… ln + 1, n 1 ⋮ ⋱ ⋱ l N, 1… l N, n… l N, N - 1 1). {\ displaystyle L = {\ begin {pmatrix} 1 \\ l_ {2,1} \ ddots \\ \ ddots 1 \\\ vdots \ dots l_ {n + 1, n} 1 \\ \ vdots \ ddots \ ddots \\ l_ {N, 1} \ dots l_ {N, n} \ dots l_ {N, N-1} 1 \ end {pmatrix}}.}

{\ displaystyle L = {\ begin {p matrix} 1 \\ l_ {2,1} \ ddots \\ \ ddots 1 \\\ vdots \ dots l_ {n + 1, n} 1 \\ \ vdots \ ddots \ ddots \ \ l_ {N, 1} \ dots l_ {N, n} \ dots l_ {N, N-1} 1 \ end {pmatrix}}.}

Получаем $А = LU. {\ displaystyle A = LU.}$ ${\ displaystyle A = LU.}$

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

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

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

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

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

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

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

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

$P 1 A = (1 0 cv I n - 1) (aw T 0 A ′ - cvw T) {\ 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 ^ {T} \\\ hline \\ 0 A'-cvw ^ {T} \\ \ end {array}} \ right)}$ $P_{1}A=\left({\begin{array}{c|ccc}10\\\hline \\cvI_{n-1}\\\end{array}}\right)\left({\begin{array}{c|c}aw^{T}\\\hline \\0A'-cvw^{T}\\\end{array}}\right)$ .

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

$(1 0 0 P ′) P 1 A = (1 0 cv ′ L ′) (aw T 0 U ′) {\ 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 ^ {T} \\\ hline \\ 0 U' \\ \ end {array}} \ right)}$ $\left({\begin{array}{c|ccc}10\\\hline \\0P'\\\end{array}}\right)P_{1}A=\left({\begin{array}{c|ccc}10\\\hline \\cv'L'\\\end{array}}\right)\left({\begin{array}{c|ccc}aw^{T}\\\hline \\0U'\\\end{array}}\right)$ ,

, который является LUP-разложением A.

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

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

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

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

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

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

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

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

приложений

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

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

A x = b, {\ displaystyle Ax = b,}

Ax = b,

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

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

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

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

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

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

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

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

AX = ЛЮКС = B. {\ displaystyle AX = LUX = B.}

{\ displaystyle AX = LUX = B.}

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

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

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

det (A) = det (P - 1) det (L) det (U) = ( - 1) S (∏ i = 1 nlii) (∏ i = 1 nuii). {\ Displaystyle \ Det (A) = \ Det (P ^ {- 1}) \ Det (L) \ Det (U) = (- 1) ^ {S} \ left (\ prod _ {i = 1} ^ {n} l_ {ii} \ right) \ left (\ prod _ {i = 1} ^ {n} u_ {ii} \ right).}

{\ displaystyle \ det (A) = \ det (P ^ {- 1}) \ 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 с полным поворотом, $det (A) {\ textstyle \ det (A)}$ ${\ textstyle \ det (A)}$ также равен правой части приведенного выше уравнения, если мы позволим S быть общим количеством обменов строк и столбцов.

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

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

/ * INPUT: A - массив указателей на строки квадратной матрицы, имеющей размерность N * Tol - малое число допуска для обнаружения отказа, когда матрица близка к вырождению * OUTPUT: Matrix 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; двойной maxA, * ptr, absA; для (i = 0; i <= N; i++) P[i] = i; //Unit permutation matrix, P[N] initialized with N for (i = 0; i < N; i++) { maxA = 0.0; imax = i; for (k = i; k < N; k++) if ((absA = fabs(A[k][i]))>maxA) {maxA = absA; imax = k; } if (maxA < Tol) return 0; //failure, matrix is degenerate if (imax != i) { //pivoting P j = P[i]; P[i] = P[imax]; P[imax] = j; //pivoting rows of A ptr = A[i]; A[i] = A[imax]; A[imax] = ptr; //counting pivots starting from N (for determinant) P[N]++; } for (j = i + 1; j < N; j++) { A[j][i] /= A[i][i]; for (k = i + 1; k < N; k++) A[j][k] -= A[j][i] * A[i][k]; } } return 1; //decomposition done } /* INPUT: A,P filled in LUPDecompose; b - rhs vector; N - dimension * OUTPUT: x - solution vector of A*x=b */ void LUPSolve(double **A, int *P, double *b, int N, double *x) { for (int i = 0; i < N; i++) { x[i] = b[P[i]]; for (int k = 0; k < i; k++) x[i] -= A[i][k] * x[k]; } for (int i = N - 1; i>= 0; i--) {for (int k = i + 1; k < N; k++) x[i] -= A[i][k] * x[k]; x[i] /= A[i][i]; } } /* INPUT: A,P filled in LUPDecompose; N - dimension * OUTPUT: IA is the inverse of the initial matrix */ void LUPInvert(double **A, int *P, int N, double **IA) { for (int j = 0; j < N; j++) { for (int i = 0; i < N; i++) { IA[i][j] = P[i] == j ? 1.0 : 0.0 for (int k = 0; k < i; k++) IA[i][j] -= A[i][k] * IA[k][j]; } for (int i = N - 1; i>= 0; i--) {for (int k = i + 1; k < N; k++) IA[i][j] -= A[i][k] * IA[k][j]; IA[i][j] /= A[i][i]; } } } /* INPUT: A,P filled in LUPDecompose; N - dimension. * OUTPUT: Function returns the determinant of the initial matrix */ double LUPDeterminant(double **A, int *P, int N) { double det = A[0][0]; for (int i = 1; i < N; i++) det *= A[i][i]; return (P[N] - N) % 2 == 0 ? det : -det }

C # примеры кода

открытый класс 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; for (int k = 0; k < i; k++) sum += lu[j, k] * lu[k, i]; lu[j, i] = (1 / lu[i, i]) * (matrix[j, i] - sum); } } // lu = L+U-I // find solution of Ly = b double y = new double[n]; for (int i = 0; i < n; i++) { sum = 0; for (int k = 0; k < i; k++) sum += lu[i, k] * y[k]; y[i] = rightPart[i] - sum; } // find solution of Ux = y double x = new double[n]; for (int i = n - 1; i>= 0; i--) {sum = 0; for (int k = i + 1; k < n; k++) sum += lu[i, k] * x[k]; x[i] = (1 / lu[i, i]) * (y[i] - sum); } return x; } }

Примеры кода MATLAB

function x = SolveLinearSystem (A, b) n = длина (b); x = нули (n, 1); y = нули (n, 1);% разложение матрицы, метод Дулиттла для i = 1: 1: n для j = 1: 1: (i - 1) альфа = A (i, j); для k = 1: 1: (j - 1) alpha = alpha - A (i, k) * A (k, j); конец A (i, j) = alpha / A (j, j); конец для j = i: 1: n alpha = A (i, j); для k = 1: 1: (i - 1) alpha = alpha - A (i, k) * A (k, j); end A (i, j) = alpha; end end% A = L + UI% найти решение Ly = b для i = 1: 1: n alpha = 0; для k = 1: 1 : i alpha = alpha + A (i, k) * y (k); end y (i) = b (i) - alpha; end% найти решение Ux = y для i = n: (- 1): 1 альфа = 0; для k = (i + 1): 1: n альфа = альфа + A (i, к) * х (к); конец x (i) = (y (i) - альфа) / A (i, i); end end

См. также

Примечания

Ссылки

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

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

Ссылки

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

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

LAPACK - это набор подпрограмм FORTRAN для решения задач плотной линейной алгебры.
ALGLIB включает частичный порт LAPACK для C ++, C #, Delphi и т. Д.
Код C ++, профессор Дж. Лумис, Университет Дейтона
Код C, Математическая исходная библиотека
LU в X10
[1] Лю в C ++ Анил Педгаонкар
Рандомизированный код LU MATLAB

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