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

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

Определения

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

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

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

Без надлежащего упорядочивания или перестановок в матрице факторизация может не осуществиться. Например, легко проверить (раскрывая матричное умножение), что . Если , то хотя бы один из и должно быть нулем, что означает, что либо L или же U является единственное число. Это невозможно, если А неособо (обратимо). Это процедурная проблема. Его можно удалить, просто изменив порядок строк А так что первый элемент переставленной матрицы отличен от нуля. Та же проблема на последующих этапах факторизации может быть устранена таким же образом; см. основную процедуру ниже.

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

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

куда L и U - снова нижняя и верхняя треугольные матрицы, и п это матрица перестановок, которая при умножении слева на А, меняет порядок строк А. Оказывается, все квадратные матрицы можно разложить на множители в таком виде:[2] и факторизация численно стабильна на практике.[3] Это делает разложение LUP полезным методом на практике.

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

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

куда L, U и п определены, как и раньше, и Q матрица перестановки, которая переупорядочивает столбцы А.[4]

Разложение LDU

An Разложение LDU является разложением вида

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

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

Пример

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

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

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

Подставляя эти значения в разложение LU выше, получаем


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

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

Любая квадратная матрица признает LUP и PLU факторизации.[2] Если является обратимый, то допускает LU (или же LDU) факторизация если и только если весь его ведущий руководитель несовершеннолетние ненулевые.[5] Если особая матрица ранга , то допускает LU факторизация, если первая ведущие главные миноры отличны от нуля, хотя обратное неверно.[6]

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

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

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

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

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

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

Алгоритмы

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

Когда факторизация LDU существует и уникальна, существует замкнутая (явная) формула для элементов L, D, и U через отношения определителей некоторых подматриц исходной матрицы А.[8] Особенно, , и для , это соотношение -я главная подматрица -я главная подматрица. Вычисление определителей вычислительно дорогой, поэтому эта явная формула на практике не используется.

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

Следующий алгоритм представляет собой модифицированную форму Гауссово исключение. Для вычисления разложения LU с использованием этого алгоритма требуется операции с плавающей запятой, игнорируя младшие члены. Частичное поворот добавляет только квадратичный член; это не относится к полному повороту.[9]

Учитывая N × N матрица , определять .

Мы исключаем элементы матрицы ниже главной диагонали в п-й столбец А(п − 1) добавив к я-й строке этой матрицы п-я строка умноженная на

Это можно сделать, умножив А(п − 1) слева с нижней треугольной матрицей

это N Икс N единичная матрица с ее п-й столбец заменен вектором

Мы установили

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

Обозначим верхнетреугольную матрицу А(N − 1) к U, и . Поскольку обратная нижнетреугольная матрица Lп снова является нижнетреугольной матрицей, а умножение двух нижнетреугольных матриц снова является нижнетреугольной матрицей, из этого следует, что L - нижнетреугольная матрица. Более того, видно, что

Мы получаем

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

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

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

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

Cormen et al.[10] описать рекурсивный алгоритм разложения LUP.

Учитывая матрицу А, позволять п1 матрица перестановок такая, что

,

куда , если в первом столбце А; или взять п1 в качестве единичной матрицы в противном случае. Теперь позвольте , если ; или же иначе. У нас есть

.

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

,

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

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

Можно найти приближение низкого ранга к разложению LU, используя рандомизированный алгоритм. Учитывая входную матрицу и желаемый низкий ранг , рандомизированный LU возвращает матрицы перестановок и нижняя / верхняя трапециевидные матрицы размера и соответственно такие, что с большой вероятностью , куда - константа, зависящая от параметров алгоритма и это -е сингулярное значение входной матрицы .[11]

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

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

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

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

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

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

Приложения

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

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

мы хотим решить уравнение для Икс, данный А и б. Предположим, мы уже получили LUP-разложение А такой, что , так .

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

  1. Сначала решаем уравнение за у.
  2. Во-вторых, решаем уравнение за Икс.

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

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

Стоимость решения системы линейных уравнений составляет примерно операции с плавающей запятой, если матрица имеет размер . Это делает его вдвое быстрее, чем алгоритмы, основанные на QR-разложение, что стоит около операции с плавающей запятой, когда Размышления домохозяина используются. По этой причине обычно предпочтительнее разложение LU.[13]

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

При решении систем уравнений б обычно рассматривается как вектор с длиной, равной высоте матрицы А. Однако при обращении матрицы вместо вектора б, имеем матрицу B, куда B является п-к-п матрица, так что мы пытаемся найти матрицу Икс (также п-к-п матрица):

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

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

Учитывая разложение LUP квадратной матрицы А, определитель А можно просто вычислить как

Второе уравнение следует из того факта, что определитель треугольной матрицы - это просто произведение ее диагональных элементов, и что определитель матрицы перестановок равен (−1)S куда S количество замен строк в разложении.

В случае разложения ЛУ с полным поворотом также равна правой части приведенного выше уравнения, если мы положим S быть общим количеством обменов строк и столбцов.

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

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

/ * INPUT: A - массив указателей на строки квадратной матрицы размерности N * Tol - небольшое число допуска для обнаружения отказа, когда матрица близка к вырожденной * ВЫХОД: Матрица A изменена, она содержит копию обеих матриц L-E и U как A = (L-E) + U, так что P * A = L * U. * Матрица перестановок хранится не как матрица, а в целочисленном векторе P размера N + 1  * содержащий индексы столбцов, в которых матрица перестановок имеет "1". Последний элемент P [N] = S + N,  * где S - количество обменов строк, необходимых для вычисления определителя, det (P) = (- 1) ^ S  */int LUPDecompose(двойной **А, int N, двойной Тол, int *п) {    int я, j, k, imax;     двойной maxA, *ptr, абсА;    за (я = 0; я <= N; я++)        п[я] = я; // Матрица единичной перестановки, P [N] инициализирована N    за (я = 0; я < N; я++) {        maxA = 0.0;        imax = я;        за (k = я; k < N; k++)            если ((абсА = фабрики(А[k][я])) > maxA) {                 maxA = абсА;                imax = k;            }        если (maxA < Тол) возвращаться 0; // сбой, матрица вырожденная        если (imax != я) {            // поворот P            j = п[я];            п[я] = п[imax];            п[imax] = j;            // поворот строк A            ptr = А[я];            А[я] = А[imax];            А[imax] = ptr;            // подсчет опорных точек, начиная с N (для определителя)            п[N]++;        }        за (j = я + 1; j < N; j++) {            А[j][я] /= А[я][я];            за (k = я + 1; k < N; k++)                А[j][k] -= А[j][я] * А[я][k];        }    }    возвращаться 1;  // разложение выполнено }/ * INPUT: A, P заполнены LUPDecompose; b - вектор rhs; N - размер * ВЫХОД: x - вектор решения A * x = b */пустота LUPSolve(двойной **А, int *п, двойной *б, int N, двойной *Икс) {    за (int я = 0; я < N; я++) {        Икс[я] = б[п[я]];        за (int k = 0; k < я; k++)            Икс[я] -= А[я][k] * Икс[k];    }    за (int я = N - 1; я >= 0; я--) {        за (int k = я + 1; k < N; k++)            Икс[я] -= А[я][k] * Икс[k];        Икс[я] /= А[я][я];    }}/ * INPUT: A, P заполнены LUPDecompose; N - размер * ВЫХОД: IA - это обратная исходная матрица */пустота LUPInvert(двойной **А, int *п, int N, двойной **Я) {      за (int j = 0; j < N; j++) {        за (int я = 0; я < N; я++) {            Я[я][j] = п[я] == j ? 1.0 : 0.0            за (int k = 0; k < я; k++)                Я[я][j] -= А[я][k] * Я[k][j];        }        за (int я = N - 1; я >= 0; я--) {            за (int k = я + 1; k < N; k++)                Я[я][j] -= А[я][k] * Я[k][j];            Я[я][j] /= А[я][я];        }    }}/ * INPUT: A, P заполнено LUPDecompose; N - размерность.  * ВЫХОД: функция возвращает определитель исходной матрицы */двойной LUP Детерминант(двойной **А, int *п, int N) {    двойной Det = А[0][0];    за (int я = 1; я < N; я++)        Det *= А[я][я];    возвращаться (п[N] - N) % 2 == 0 ? Det : -Det}

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

    общественный учебный класс SystemOfLinearEquations    {        общественный двойной[] РешитьUsingLU(двойной[,] матрица, двойной[] rightPart, int п)        {            // разложение матрицы            двойной[,] Лу = новый двойной[п, п];            двойной сумма = 0;            за (int я = 0; я < п; я++)            {                за (int j = я; j < п; j++)                {                    сумма = 0;                    за (int k = 0; k < я; k++)                        сумма += Лу[я, k] * Лу[k, j];                    Лу[я, j] = матрица[я, j] - сумма;                }                за (int j = я + 1; j < п; j++)                {                    сумма = 0;                    за (int k = 0; k < я; k++)                        сумма += Лу[j, k] * Лу[k, я];                    Лу[j, я] = (1 / Лу[я, я]) * (матрица[j, я] - сумма);                }            }                        // lu = L + U-I            // находим решение Ly = b            двойной[] у = новый двойной[п];            за (int я = 0; я < п; я++)            {                сумма = 0;                за (int k = 0; k < я; k++)                    сумма += Лу[я, k] * у[k];                у[я] = rightPart[я] - сумма;            }            // находим решение Ux = y            двойной[] Икс = новый двойной[п];            за (int я = п - 1; я >= 0; я--)            {                сумма = 0;                за (int k = я + 1; k < п; k++)                    сумма += Лу[я, k] * Икс[k];                Икс[я] = (1 / Лу[я, я]) * (у[я] - сумма);            }            возвращаться Икс;        }    }

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

функцияИкс =SolveLinearSystem(А, б)п = длина(б);    Икс = нули(п, 1);    у = нули(п, 1);    % разложения матрицы, метод Дулиттла    за я = 1: 1: п        за j = 1: 1: (я - 1)            альфа = А(я,j);            за к = 1: 1: (j - 1)                альфа = альфа - А(я,k)*А(k,j);            конецA (i, j) = альфа / A (j, j);        конецдля j = i: 1: n            альфа = А(я,j);            за к = 1: 1: (я - 1)                альфа = альфа - А(я,k)*А(k,j);            конецA (i, j) = альфа;        конецконец    % A = L + U-I    % найти решение Ly = b    за я = 1: 1: п        альфа = 0;        за к = 1: 1: я            альфа = альфа + А(я,k)*у(k);        конецy (i) = b (i) - альфа;    конец% найти решение Ux = y    за я = п: (- 1): 1        альфа = 0;        за к = (я + 1): 1: п            альфа = альфа + А(я,k)*Икс(k);        конецх (я) = (у (я) - альфа) / А (я, я);    конец конец

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

Примечания

  1. ^ Шварценберг-Черни, А. (1995). «О матричной факторизации и эффективном решении методом наименьших квадратов». Серия дополнений по астрономии и астрофизике. 110: 405. Bibcode:1995A и AS..110..405S.
  2. ^ а б Окунев и Джонсон (1997), Следствие 3.
  3. ^ Trefethen & Bau (1997), п. 166.
  4. ^ Trefethen & Bau (1997), п. 161.
  5. ^ а б Хорн и Джонсон (1985), Следствие 3.5.5
  6. ^ Хорн и Джонсон (1985), Теорема 3.5.2
  7. ^ Окунев и Джонсон (1997)
  8. ^ Домохозяин (1975)
  9. ^ Голуб и Ван Лоан (1996), п. 112, 119.
  10. ^ Кормен, Томас Х.; Лейзерсон, Чарльз Э.; Ривест, Рональд Л.; Штейн, Клиффорд (2001). Введение в алгоритмы. MIT Press и McGraw-Hill. ISBN  978-0-262-03293-3.
  11. ^ Шабат, Гиль; Шмуэли, Янив; Айзенбуд, Ярив; Авербух, Амир (2016). «Рандомизированное разложение LU». Прикладной и вычислительный гармонический анализ. 44 (2): 246–272. arXiv:1310.7202. Дои:10.1016 / j.acha.2016.04.006.
  12. ^ Банч и Хопкрофт (1974)
  13. ^ Trefethen & Bau (1997), п. 152.
  14. ^ Голуб и Ван Лоан (1996), п. 121

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

внешняя ссылка

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

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

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