Разложение Холецкого - Cholesky decomposition

В линейная алгебра, то Разложение Холецкого или же Факторизация Холецкого (произносится /ʃə.ˈлɛs.kя/) это разложение из Эрмитский, положительно определенная матрица в продукт нижняя треугольная матрица и это сопряженный транспонировать, что полезно для эффективных численных решений, например, Моделирование Монте-Карло. Это было обнаружено Андре-Луи Холески для реальных матриц. Когда это применимо, разложение Холецкого примерно в два раза эффективнее, чем LU разложение для решения системы линейных уравнений.^[1]

Заявление

Разложение Холецкого Эрмитский положительно определенная матрица А, является разложением вида

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

куда L это нижняя треугольная матрица с действительными и положительными диагональными элементами, и L* обозначает сопряженный транспонировать из L. Каждая эрмитова положительно определенная матрица (а значит, и каждая вещественнозначная симметричная положительно определенная матрица) имеет уникальное разложение Холецкого.^[2]

Обратное утверждение тривиально: если А можно записать как LL* для некоторых обратимых L, нижний треугольник или иное, то А является эрмитовым и положительно определенным.

Когда А является вещественной матрицей (следовательно, симметричной положительно определенной), факторизация может быть записана

А = LL^Т,

куда L - вещественная нижнетреугольная матрица с положительными диагональными элементами.^[3][4][5]

Положительные полуопределенные матрицы

Если эрмитова матрица А является только положительно полуопределенным, а не положительно определенным, то он все еще имеет разложение вида А = LL* где диагональные элементы L могут быть равны нулю.^[6]Разложение не обязательно должно быть уникальным, например:

${\displaystyle {\begin{pmatrix}0&0\\0&1\end{pmatrix}}=\mathbf {L} \mathbf {L} ^{*},\quad \quad \mathbf {L} ={\begin{pmatrix}0&0\\\cos \theta &\sin \theta \end{pmatrix}}.}$

Однако если ранг А является р, то существует единственный нижний треугольник L с точно р положительные диагональные элементы и н-р столбцы, содержащие все нули.^[7]

В качестве альтернативы декомпозиция может быть сделана уникальной, если фиксирован выбор поворота. Формально, если А является п × п положительно полуопределенная матрица ранга р, то существует хотя бы одна матрица перестановок п такой, что П А П^Т имеет единственное разложение вида П А П^Т = L L^* с${\displaystyle \mathbf {L} ={\begin{bmatrix}\mathbf {L} _{1}&0\\\mathbf {L} _{2}&0\end{bmatrix}}}$,куда L₁ является р × р нижнетреугольная матрица с положительной диагональю.^[8]

Разложение ЛПНП

Близким вариантом классического разложения Холецкого является разложение ЛПНП,

${\displaystyle \mathbf {A} =\mathbf {LDL} ^{*},}$

куда L это нижний блок треугольный (унитреугольный) матрица и D это диагональ матрица, то есть диагональные элементы L должны быть равны 1 за счет введения дополнительной диагональной матрицы D Основное преимущество состоит в том, что разложение LDL можно вычислить и использовать по существу с теми же алгоритмами, но без извлечения квадратных корней.^[9]

По этой причине разложение ЛПНП часто называют без квадратного корня Холецкий разложение. Для вещественных матриц факторизация имеет вид А = ЛПНП^Т и часто упоминается как Разложение LDLT (или ЛПНП^Т разложение, или ЛПНП '). Это тесно связано с собственное разложение вещественных симметричных матриц, А = QΛQ^Т.

Разложение ЛПНП связано с классическим разложением Холецкого вида LL* следующее:

${\displaystyle \mathbf {A} =\mathbf {LDL} ^{*}=\mathbf {L} \mathbf {D} ^{1/2}(\mathbf {D} ^{1/2})^{*}\mathbf {L} ^{*}=\mathbf {L} \mathbf {D} ^{1/2}(\mathbf {L} \mathbf {D} ^{1/2})^{*}.}$

Наоборот, учитывая классическое разложение Холецкого ${\displaystyle \mathbf {A} =\mathbf {C} \mathbf {C} ^{*}}$ положительно определенной матрицы, если S - диагональная матрица, содержащая главную диагональ ${\displaystyle \mathbf {C} }$, затем А можно разложить как ${\displaystyle \mathbf {L} \mathbf {D} \mathbf {L} ^{*}}$ куда

${\displaystyle \mathbf {L} =\mathbf {C} \mathbf {S} ^{-1}}$ (это изменяет масштаб каждого столбца, чтобы сделать диагональные элементы равными 1),

${\displaystyle \mathbf {D} =\mathbf {S} ^{2}.}$

Если А положительно определена, то диагональные элементы D все положительны. для положительных полуопределенных А, ${\displaystyle \mathbf {L} \mathbf {D} \mathbf {L} ^{*}}$ разложение существует там, где количество ненулевых элементов на диагонали D это в точности ранг А.^[10]Некоторые неопределенные матрицы, для которых не существует разложения Холецкого, имеют разложение LDL с отрицательными элементами в D: достаточно, чтобы первая п-1 ведущие основные несовершеннолетние из А неособые.^[11]

Пример

Вот разложение Холецкого симметричной вещественной матрицы:

${\displaystyle {\begin{aligned}\left({\begin{array}{*{3}{r}}4&12&-16\\12&37&-43\\-16&-43&98\\\end{array}}\right)=\left({\begin{array}{*{3}{r}}2&0&0\\6&1&0\\-8&5&3\\\end{array}}\right)\left({\begin{array}{*{3}{r}}2&6&-8\\0&1&5\\0&0&3\\\end{array}}\right).\end{aligned}}}$

А вот его ЛПНП^Т разложение:

${\displaystyle {\begin{aligned}\left({\begin{array}{*{3}{r}}4&12&-16\\12&37&-43\\-16&-43&98\\\end{array}}\right)&=\left({\begin{array}{*{3}{r}}1&0&0\\3&1&0\\-4&5&1\\\end{array}}\right)\left({\begin{array}{*{3}{r}}4&0&0\\0&1&0\\0&0&9\\\end{array}}\right)\left({\begin{array}{*{3}{r}}1&3&-4\\0&1&5\\0&0&1\\\end{array}}\right).\end{aligned}}}$

Приложения

Разложение Холецкого в основном используется для численного решения линейные уравнения ${\displaystyle \mathbf {Ax} =\mathbf {b} }$. Если А симметрично и положительно определено, то мы можем решить ${\displaystyle \mathbf {Ax} =\mathbf {b} }$ сначала вычислив разложение Холецкого ${\displaystyle \mathbf {A} =\mathbf {LL} ^{\mathrm {*} }}$, затем решая ${\displaystyle \mathbf {Ly} =\mathbf {b} }$ за у к форвардная замена, и наконец решение ${\displaystyle \mathbf {L^{*}x} =\mathbf {y} }$ за Икс к обратная замена.

Альтернативный способ избавиться от извлечения квадратного корня из ${\displaystyle \mathbf {LL} ^{\mathrm {*} }}$ разложение заключается в вычислении разложения Холецкого ${\displaystyle \mathbf {A} =\mathbf {LDL} ^{\mathrm {*} }}$, затем решая ${\displaystyle \mathbf {Ly} =\mathbf {b} }$ за у, и наконец решение ${\displaystyle \mathbf {DL} ^{\mathrm {*} }\mathbf {x} =\mathbf {y} }$.

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

Линейный метод наименьших квадратов

Системы формы Топор = б с А симметричные и положительно определенные возникают в приложениях довольно часто. Например, нормальные уравнения в линейный метод наименьших квадратов проблемы имеют такую ​​форму. Также может случиться, что матрица А исходит от энергетического функционала, который должен быть положительным с физических соображений; это часто случается при численном решении уравнения в частных производных.

Нелинейная оптимизация

Нелинейные многомерные функции могут быть минимизированы по их параметрам, используя варианты Метод Ньютона называется квазиньютон методы. На итерации k поиск идет в направлении ${\displaystyle p_{k}}$ определяется путем решения ${\displaystyle B_{k}p_{k}}$ = ${\displaystyle -g_{k}}$ за ${\displaystyle p_{k}}$, куда ${\displaystyle p_{k}}$ это направление шага, ${\displaystyle g_{k}}$ это градиент, и ${\displaystyle B_{k}}$ является приближением к Матрица Гессе формируется путем повторения обновлений ранга 1 на каждой итерации. Две известные формулы обновления называются Дэвидон – Флетчер – Пауэлл (DFP) и Бройден – Флетчер – Гольдфарб – Шанно (BFGS). Потери положительно-определенного условия из-за ошибки округления можно избежать, если вместо обновления приближения к обратному к гессиану обновлять разложение Холецкого аппроксимации самой матрицы Гессе.^[12]

Моделирование Монте-Карло

Разложение Холецкого обычно используется в Метод Монте-Карло для моделирования систем с множеством коррелированных переменных. В ковариационная матрица раскладывается, чтобы получить нижнетреугольную L. Применяя это к вектору некоррелированных выборок ты производит образец вектора Лу с ковариационными свойствами моделируемой системы.^[13]

Следующий упрощенный пример показывает экономию, которую можно получить из разложения Холецкого: предположим, цель состоит в том, чтобы сгенерировать две коррелированные нормальные переменные. ${\displaystyle x_{1}}$ и ${\displaystyle x_{2}}$ с заданным коэффициентом корреляции ${\displaystyle \rho }$. Для этого необходимо сначала сгенерировать две некоррелированные гауссовские случайные величины. ${\displaystyle z_{1}}$ и ${\displaystyle z_{2}}$, что можно сделать с помощью Преобразование Бокса – Мюллера. С учетом необходимого коэффициента корреляции ${\displaystyle \rho }$, коррелированные нормальные переменные могут быть получены с помощью преобразований ${\displaystyle x_{1}=z_{1}}$ и ${\displaystyle x_{2}=\rho z_{1}+{\sqrt {1-\rho ^{2}}}z_{2}}$.

Фильтры Калмана

Фильтры Калмана без запаха обычно используют разложение Холецкого для выбора набора так называемых сигма-точек. Фильтр Калмана отслеживает среднее состояние системы как вектор Икс длины N и ковариация как N × N матрица п. Матрица п всегда положительно полуопределен и может быть разложен на LL^Т. Столбцы L можно складывать и вычитать из среднего Икс чтобы сформировать набор из 2N векторы, называемые сигма точки. Эти сигма-точки полностью отражают среднее значение и ковариацию состояния системы.

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

Явный обратный эрмитовой матрицы можно вычислить с помощью разложения Холецкого аналогично решению линейных систем с использованием ${\displaystyle n^{3}}$ операции (${\displaystyle {\tfrac {1}{2}}n^{3}}$ умножения).^[9] Вся инверсия может быть даже эффективно выполнена на месте.

Неэрмитова матрица B также можно инвертировать, используя следующее тождество, где BB* всегда будет эрмитским:

${\displaystyle \mathbf {B} ^{-1}=\mathbf {B} ^{*}(\mathbf {BB} ^{*})^{-1}.}$

Вычисление

Существуют различные методы вычисления разложения Холецкого. Вычислительная сложность обычно используемых алгоритмов составляет О(п³) в целом.^{[нужна цитата ]} Все описанные ниже алгоритмы включают около п³/3 FLOPs (п³/ 6 умножений и такое же количество сложений), где п размер матрицы А. Следовательно, они получают половину стоимости LU разложение, который использует 2п³/ 3 FLOP (см. Trefethen and Bau 1997).

Какой из приведенных ниже алгоритмов быстрее, зависит от деталей реализации. Как правило, первый алгоритм будет немного медленнее, потому что он получает доступ к данным менее регулярно.

Алгоритм Холецкого

В Алгоритм Холецкого, используется для вычисления матрицы разложения L, это модифицированная версия Гауссово исключение.

Рекурсивный алгоритм начинается с я : = 1 и

А⁽¹⁾ := А.

На шаге я, матрица А^(я) имеет следующий вид:

${\displaystyle \mathbf {A} ^{(i)}={\begin{pmatrix}\mathbf {I} _{i-1}&0&0\\0&a_{i,i}&\mathbf {b} _{i}^{*}\\0&\mathbf {b} _{i}&\mathbf {B} ^{(i)}\end{pmatrix}},}$

куда я_я−1 обозначает единичная матрица измерения я − 1.

Если теперь определить матрицу L_я к

${\displaystyle \mathbf {L} _{i}:={\begin{pmatrix}\mathbf {I} _{i-1}&0&0\\0&{\sqrt {a_{i,i}}}&0\\0&{\frac {1}{\sqrt {a_{i,i}}}}\mathbf {b} _{i}&\mathbf {I} _{n-i}\end{pmatrix}},}$

тогда мы можем написать А^(я) в качестве

${\displaystyle \mathbf {A} ^{(i)}=\mathbf {L} _{i}\mathbf {A} ^{(i+1)}\mathbf {L} _{i}^{*}}$

куда

${\displaystyle \mathbf {A} ^{(i+1)}={\begin{pmatrix}\mathbf {I} _{i-1}&0&0\\0&1&0\\0&0&\mathbf {B} ^{(i)}-{\frac {1}{a_{i,i}}}\mathbf {b} _{i}\mathbf {b} _{i}^{*}\end{pmatrix}}.}$

Обратите внимание, что б_я б*_я является внешний продукт, поэтому этот алгоритм называется внешняя версия продукта в (Голуб и Ван Лоан).

Мы повторяем это для я от 1 до п. После п шаги, мы получаем А^(п+1) = я. Следовательно, нижнетреугольная матрица L мы ищем рассчитывается как

${\displaystyle \mathbf {L} :=\mathbf {L} _{1}\mathbf {L} _{2}\dots \mathbf {L} _{n}.}$

Алгоритмы Холецкого – Банахевича и Холецкого – Краута.

Шаблон доступа (белый) и шаблон записи (желтый) для локального алгоритма Холецкого-Банахевича на матрице 5 × 5

Если мы выпишем уравнение

${\displaystyle {\begin{aligned}\mathbf {A} =\mathbf {LL} ^{T}&={\begin{pmatrix}L_{11}&0&0\\L_{21}&L_{22}&0\\L_{31}&L_{32}&L_{33}\\\end{pmatrix}}{\begin{pmatrix}L_{11}&L_{21}&L_{31}\\0&L_{22}&L_{32}\\0&0&L_{33}\end{pmatrix}}\\[8pt]&={\begin{pmatrix}L_{11}^{2}&&({\text{symmetric}})\\L_{21}L_{11}&L_{21}^{2}+L_{22}^{2}&\\L_{31}L_{11}&L_{31}L_{21}+L_{32}L_{22}&L_{31}^{2}+L_{32}^{2}+L_{33}^{2}\end{pmatrix}},\end{aligned}}}$

получаем следующее:

${\displaystyle {\begin{aligned}\mathbf {L} ={\begin{pmatrix}{\sqrt {A_{11}}}&0&0\\A_{21}/L_{11}&{\sqrt {A_{22}-L_{21}^{2}}}&0\\A_{31}/L_{11}&\left(A_{32}-L_{31}L_{21}\right)/L_{22}&{\sqrt {A_{33}-L_{31}^{2}-L_{32}^{2}}}\end{pmatrix}}\end{aligned}}}$

и поэтому следующие формулы для элементов L:

${\displaystyle L_{j,j}=(\pm ){\sqrt {A_{j,j}-\sum _{k=1}^{j-1}L_{j,k}^{2}}},}$

${\displaystyle L_{i,j}={\frac {1}{L_{j,j}}}\left(A_{i,j}-\sum _{k=1}^{j-1}L_{i,k}L_{j,k}\right)\quad {\text{for }}i>j.}$

Для комплексных и вещественных матриц разрешены несущественные изменения знаков диагональных и связанных недиагональных элементов. Выражение под квадратный корень всегда положительно, если А является реальным и положительно определенным.

Для комплексной эрмитовой матрицы применяется следующая формула:

${\displaystyle L_{j,j}={\sqrt {A_{j,j}-\sum _{k=1}^{j-1}L_{j,k}L_{j,k}^{*}}},}$

${\displaystyle L_{i,j}={\frac {1}{L_{j,j}}}\left(A_{i,j}-\sum _{k=1}^{j-1}L_{i,k}L_{j,k}^{*}\right)\quad {\text{for }}i>j.}$

Итак, мы можем вычислить (я, j) запись, если нам известны записи слева и выше. Вычисления обычно располагаются в одном из следующих порядков:

• В Алгоритм Холецкого – Банахевича начинается с верхнего левого угла матрицы L и переходит к вычислению матрицы построчно.
• В Алгоритм Холецкого – Краута начинается с верхнего левого угла матрицы L и переходит к вычислению матрицы столбец за столбцом.

Любой из вариантов доступа позволяет при желании выполнять все вычисления на месте.

Стабильность расчета

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

Теперь предположим, что применимо разложение Холецкого. Как было сказано выше, алгоритм будет вдвое быстрее. Кроме того, нет поворот необходимо, и погрешность всегда будет небольшой. В частности, если мы хотим решить Топор = б, и у обозначает вычисленное решение, то у решает возмущенную систему (А + E)у = б, куда

${\displaystyle \|\mathbf {E} \|_{2}\leq c_{n}\varepsilon \|\mathbf {A} \|_{2}.}$

Здесь || · ||₂ это матрица 2-норма, c_п небольшая константа, зависящая от п, а ε обозначает округление единицы.

Одна из проблем разложения Холецкого, о которой следует помнить, - это использование квадратных корней. Если факторизуемая матрица положительно определена, как требуется, числа под квадратными корнями всегда положительны. в точной арифметике. К сожалению, числа могут стать отрицательными из-за ошибки округления, и в этом случае алгоритм не может продолжаться. Однако это может произойти только в том случае, если матрица очень плохо подготовлена. Один из способов решить эту проблему - добавить матрицу диагональной коррекции к разлагаемой матрице, чтобы попытаться обеспечить положительную определенность.^[14] Хотя это может снизить точность разложения, это может быть очень выгодно по другим причинам; например, при выполнении Метод Ньютона в оптимизации добавление диагональной матрицы может улучшить стабильность, когда она далека от оптимальной.

Разложение ЛПНП

Альтернативная форма, устраняющая необходимость извлекать квадратные корни при А симметрично, является симметричной неопределенной факторизацией^[15]

${\displaystyle {\begin{aligned}\mathbf {A} =\mathbf {LDL} ^{\mathrm {T} }&={\begin{pmatrix}1&0&0\\L_{21}&1&0\\L_{31}&L_{32}&1\\\end{pmatrix}}{\begin{pmatrix}D_{1}&0&0\\0&D_{2}&0\\0&0&D_{3}\\\end{pmatrix}}{\begin{pmatrix}1&L_{21}&L_{31}\\0&1&L_{32}\\0&0&1\\\end{pmatrix}}\\[8pt]&={\begin{pmatrix}D_{1}&&(\mathrm {symmetric} )\\L_{21}D_{1}&L_{21}^{2}D_{1}+D_{2}&\\L_{31}D_{1}&L_{31}L_{21}D_{1}+L_{32}D_{2}&L_{31}^{2}D_{1}+L_{32}^{2}D_{2}+D_{3}.\end{pmatrix}}.\end{aligned}}}$

Следующие рекурсивные отношения применяются для записей D и L:

${\displaystyle D_{j}=A_{jj}-\sum _{k=1}^{j-1}L_{jk}^{2}D_{k},}$

${\displaystyle L_{ij}={\frac {1}{D_{j}}}\left(A_{ij}-\sum _{k=1}^{j-1}L_{ik}L_{jk}D_{k}\right)\quad {\text{for }}i>j.}$

Это работает до тех пор, пока сгенерированные диагональные элементы в D оставайся ненулевым. Тогда разложение единственное. D и L реальны, если А реально.

Для сложной эрмитовой матрицы А, применяется следующая формула:

${\displaystyle D_{j}=A_{jj}-\sum _{k=1}^{j-1}L_{jk}L_{jk}^{*}D_{k},}$

${\displaystyle L_{ij}={\frac {1}{D_{j}}}\left(A_{ij}-\sum _{k=1}^{j-1}L_{ik}L_{jk}^{*}D_{k}\right)\quad {\text{for }}i>j.}$

Опять же, шаблон доступа позволяет при желании выполнять все вычисления на месте.

Вариант блока

При использовании с неопределенными матрицами ЛПНП* известно, что факторизация нестабильна без осторожного поворота;^[16] в частности, элементы факторизации могут расти произвольно. Возможное улучшение - выполнить факторизацию блочных подматриц, обычно 2 × 2:^[17]

${\displaystyle {\begin{aligned}\mathbf {A} =\mathbf {LDL} ^{\mathrm {T} }&={\begin{pmatrix}\mathbf {I} &0&0\\\mathbf {L} _{21}&\mathbf {I} &0\\\mathbf {L} _{31}&\mathbf {L} _{32}&\mathbf {I} \\\end{pmatrix}}{\begin{pmatrix}\mathbf {D} _{1}&0&0\\0&\mathbf {D} _{2}&0\\0&0&\mathbf {D} _{3}\\\end{pmatrix}}{\begin{pmatrix}\mathbf {I} &\mathbf {L} _{21}^{\mathrm {T} }&\mathbf {L} _{31}^{\mathrm {T} }\\0&\mathbf {I} &\mathbf {L} _{32}^{\mathrm {T} }\\0&0&\mathbf {I} \\\end{pmatrix}}\\[8pt]&={\begin{pmatrix}\mathbf {D} _{1}&&(\mathrm {symmetric} )\\\mathbf {L} _{21}\mathbf {D} _{1}&\mathbf {L} _{21}\mathbf {D} _{1}\mathbf {L} _{21}^{\mathrm {T} }+\mathbf {D} _{2}&\\\mathbf {L} _{31}\mathbf {D} _{1}&\mathbf {L} _{31}\mathbf {D} _{1}\mathbf {L} _{21}^{\mathrm {T} }+\mathbf {L} _{32}\mathbf {D} _{2}&\mathbf {L} _{31}\mathbf {D} _{1}\mathbf {L} _{31}^{\mathrm {T} }+\mathbf {L} _{32}\mathbf {D} _{2}\mathbf {L} _{32}^{\mathrm {T} }+\mathbf {D} _{3}\end{pmatrix}},\end{aligned}}}$

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

${\displaystyle \mathbf {D} _{j}=\mathbf {A} _{jj}-\sum _{k=1}^{j-1}\mathbf {L} _{jk}\mathbf {D} _{k}\mathbf {L} _{jk}^{\mathrm {T} },}$

${\displaystyle \mathbf {L} _{ij}=\left(\mathbf {A} _{ij}-\sum _{k=1}^{j-1}\mathbf {L} _{ik}\mathbf {D} _{k}\mathbf {L} _{jk}^{\mathrm {T} }\right)\mathbf {D} _{j}^{-1}.}$

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

Обновление разложения

На практике часто возникает задача обновить разложение Холецкого. Более подробно, разложение Холецкого уже вычислено ${\displaystyle \mathbf {A} =\mathbf {L} \mathbf {L} ^{*}}$ какой-то матрицы ${\displaystyle \mathbf {A} }$, то заменяют матрицу ${\displaystyle \mathbf {A} }$ каким-то образом в другую матрицу, скажем ${\displaystyle {\tilde {\mathbf {A} }}}$, и нужно вычислить разложение Холецкого обновленной матрицы: ${\displaystyle {\tilde {\mathbf {A} }}={\tilde {\mathbf {L} }}{\tilde {\mathbf {L} }}^{*}}$. Вопрос теперь в том, можно ли использовать разложение Холецкого ${\displaystyle \mathbf {A} }$ который был вычислен ранее для вычисления разложения Холецкого ${\displaystyle {\tilde {\mathbf {A} }}}$.

Обновление первого ранга

Конкретный случай, когда обновленная матрица ${\displaystyle {\tilde {\mathbf {A} }}}$ связана с матрицей ${\displaystyle \mathbf {A} }$ к ${\displaystyle {\tilde {\mathbf {A} }}=\mathbf {A} +\mathbf {x} \mathbf {x} ^{*}}$, известен как ранговое обновление.

Вот небольшая функция^[18] написано в Matlab синтаксис, реализующий обновление первого ранга:

функция[L] =cholupdate(L, х)п = длина(Икс);    за k = 1:п        р = sqrt(L(k, k)^2 + Икс(k)^2);        c = р / L(k, k);        s = Икс(k) / L(k, k);        L(k, k) = р;        если k < п            L((k+1):п, k) = (L((k+1):п, k) + s * Икс((k+1):п)) / c;            Икс((k+1):п) = c * Икс((k+1):п) - s * L((k+1):п, k);        конец    конецконец

Понижение первого ранга

А понижение первого ранга похоже на обновление первого ранга, за исключением того, что добавление заменяется вычитанием: ${\displaystyle {\tilde {\mathbf {A} }}=\mathbf {A} -\mathbf {x} \mathbf {x} ^{*}}$. Это работает, только если новая матрица ${\displaystyle {\tilde {\mathbf {A} }}}$ все еще положительно определен.

Код для обновления первого ранга, показанный выше, можно легко адаптировать для понижения ранга один: нужно просто заменить два добавления в назначении на р и L ((к + 1): п, к) путем вычитания.

Добавление и удаление строк и столбцов

Если у нас есть симметричная и положительно определенная матрица ${\displaystyle \mathbf {A} }$ представлен в виде блока как

${\displaystyle \mathbf {A} ={\begin{pmatrix}\mathbf {A} _{11}&\mathbf {A} _{13}\\\mathbf {A} _{13}^{\mathrm {T} }&\mathbf {A} _{33}\\\end{pmatrix}}}$

и его верхний фактор Холецкого

${\displaystyle \mathbf {L} ={\begin{pmatrix}\mathbf {L} _{11}&\mathbf {L} _{13}\\0&\mathbf {L} _{33}\\\end{pmatrix}},}$

затем для новой матрицы ${\displaystyle {\tilde {\mathbf {A} }}}$, что совпадает с ${\displaystyle \mathbf {A} }$ но с добавлением новых строк и столбцов,

${\displaystyle {\begin{aligned}{\tilde {\mathbf {A} }}&={\begin{pmatrix}\mathbf {A} _{11}&\mathbf {A} _{12}&\mathbf {A} _{13}\\\mathbf {A} _{12}^{\mathrm {T} }&\mathbf {A} _{22}&\mathbf {A} _{23}\\\mathbf {A} _{13}^{\mathrm {T} }&\mathbf {A} _{23}^{\mathrm {T} }&\mathbf {A} _{33}\\\end{pmatrix}}\end{aligned}}}$

нас интересует факторизация Холецкого ${\displaystyle {\tilde {\mathbf {A} }}}$, который мы называем ${\displaystyle {\tilde {\mathbf {S} }}}$, без прямого вычисления всего разложения.

${\displaystyle {\begin{aligned}{\tilde {\mathbf {S} }}&={\begin{pmatrix}\mathbf {S} _{11}&\mathbf {S} _{12}&\mathbf {S} _{13}\\0&\mathbf {S} _{22}&\mathbf {S} _{23}\\0&0&\mathbf {S} _{33}\\\end{pmatrix}}.\end{aligned}}}$

Письмо ${\displaystyle \mathbf {A} \setminus \mathbf {b} }$ для решения ${\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} }$, который легко найти для треугольных матриц, и ${\displaystyle {\text{chol}}(\mathbf {M} )}$ для разложения Холецкого ${\displaystyle \mathbf {M} }$, можно найти следующие отношения:

${\displaystyle {\begin{aligned}\mathbf {S} _{11}&=\mathbf {L} _{11},\\\mathbf {S} _{12}&=\mathbf {L} _{11}^{\mathrm {T} }\setminus \mathbf {A} _{12},\\\mathbf {S} _{13}&=\mathbf {L} _{13},\\\mathbf {S} _{22}&={\text{chol}}(\mathbf {A} _{22}-\mathbf {S} _{12}^{\mathrm {T} }\mathbf {S} _{12}),\\\mathbf {S} _{23}&=\mathbf {S} _{22}^{\mathrm {T} }\setminus (\mathbf {A} _{23}-\mathbf {S} _{12}^{\mathrm {T} }\mathbf {S} _{13}),\\\mathbf {S} _{33}&={\text{chol}}(\mathbf {L} _{33}^{\mathrm {T} }\mathbf {L} _{33}-\mathbf {S} _{23}^{\mathrm {T} }\mathbf {S} _{23}).\\\end{aligned}}}$

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

${\displaystyle {\begin{aligned}{\tilde {\mathbf {A} }}&={\begin{pmatrix}\mathbf {A} _{11}&\mathbf {A} _{12}&\mathbf {A} _{13}\\\mathbf {A} _{12}^{\mathrm {T} }&\mathbf {A} _{22}&\mathbf {A} _{23}\\\mathbf {A} _{13}^{\mathrm {T} }&\mathbf {A} _{23}^{\mathrm {T} }&\mathbf {A} _{33}\\\end{pmatrix}}\end{aligned}}}$

с известным разложением Холецкого

${\displaystyle {\begin{aligned}{\tilde {\mathbf {S} }}&={\begin{pmatrix}\mathbf {S} _{11}&\mathbf {S} _{12}&\mathbf {S} _{13}\\0&\mathbf {S} _{22}&\mathbf {S} _{23}\\0&0&\mathbf {S} _{33}\\\end{pmatrix}}\end{aligned}}}$

и желаете определить фактор Холецкого

${\displaystyle {\begin{aligned}\mathbf {L} &={\begin{pmatrix}\mathbf {L} _{11}&\mathbf {L} _{13}\\0&\mathbf {L} _{33}\\\end{pmatrix}}\end{aligned}}}$

матрицы ${\displaystyle \mathbf {A} }$ с удаленными строками и столбцами,

${\displaystyle {\begin{aligned}\mathbf {A} &={\begin{pmatrix}\mathbf {A} _{11}&\mathbf {A} _{13}\\\mathbf {A} _{13}^{\mathrm {T} }&\mathbf {A} _{33}\\\end{pmatrix}},\end{aligned}}}$

дает следующие правила:

${\displaystyle {\begin{aligned}\mathbf {L} _{11}&=\mathbf {S} _{11},\\\mathbf {L} _{13}&=\mathbf {S} _{13},\\\mathbf {L} _{33}&={\text{chol}}(\mathbf {S} _{33}+\mathbf {S} _{23}^{\mathrm {T} }\mathbf {S} _{23}).\end{aligned}}}$

Обратите внимание, что все приведенные выше уравнения, которые включают поиск разложения Холецкого новой матрицы, имеют вид ${\displaystyle {\tilde {\mathbf {A} }}=\mathbf {A} \pm \mathbf {x} \mathbf {x} ^{*}}$, что позволяет эффективно рассчитывать их с помощью процедур обновления и понижения, подробно описанных в предыдущем разделе.^[19]

Доказательство для положительно полуопределенных матриц

Доказательство ограничивающим аргументом

Приведенные выше алгоритмы показывают, что любая положительно определенная матрица ${\displaystyle \mathbf {A} }$ имеет разложение Холецкого. Этот результат можно распространить на положительный полуопределенный случай ограничивающим аргументом. Аргумент не является полностью конструктивным, т.е. не дает явных численных алгоритмов для вычисления факторов Холецкого.

Если ${\displaystyle \mathbf {A} }$ является ${\displaystyle n\times n}$ положительная полуопределенная матрица, то последовательность ${\displaystyle \left(\mathbf {A} _{k}\right)_{k}:=\left(\mathbf {A} +{\frac {1}{k}}\mathbf {I} _{n}\right)_{k}}$ состоит из положительно определенные матрицы. (Это является непосредственным следствием, например, теоремы о спектральном отображении для полиномиального функционального исчисления.) Кроме того,

${\displaystyle \mathbf {A} _{k}\rightarrow \mathbf {A} \quad {\text{for}}\quad k\rightarrow \infty }$

в норма оператора. Из положительно определенного случая каждый ${\displaystyle \mathbf {A} _{k}}$ имеет разложение Холецкого ${\displaystyle \mathbf {A} _{k}=\mathbf {L} _{k}\mathbf {L} _{k}^{*}}$. По свойству операторной нормы

${\displaystyle \|\mathbf {L} _{k}\|^{2}\geq \|\mathbf {L} _{k}\mathbf {L} _{k}^{*}\|=\|\mathbf {A} _{k}\|\,.}$

Так ${\displaystyle \left(\mathbf {L} _{k}\right)_{k}}$ является ограниченным множеством в Банахово пространство операторов, поэтому относительно компактный (поскольку лежащее в основе векторное пространство конечномерно). Следовательно, он имеет сходящуюся подпоследовательность, также обозначаемую ${\displaystyle \left(\mathbf {L} _{k}\right)_{k}}$, с лимитом ${\displaystyle \mathbf {L} }$. Легко проверить, что это ${\displaystyle \mathbf {L} }$ имеет желаемые свойства, т.е. ${\displaystyle \mathbf {A} =\mathbf {L} \mathbf {L} ^{*}}$, и ${\displaystyle \mathbf {L} }$ является нижним треугольником с неотрицательными диагональными элементами: для всех ${\displaystyle x}$ и ${\displaystyle y}$,

${\displaystyle \langle \mathbf {A} x,y\rangle =\left\langle \lim \mathbf {A} _{k}x,y\right\rangle =\langle \lim \mathbf {L} _{k}\mathbf {L} _{k}^{*}x,y\rangle =\langle \mathbf {L} \mathbf {L} ^{*}x,y\rangle \,.}$

Следовательно, ${\displaystyle \mathbf {A} =\mathbf {L} \mathbf {L} ^{*}}$. Поскольку основное векторное пространство конечномерно, все топологии в пространстве операторов эквивалентны. Так ${\displaystyle \left(\mathbf {L} _{k}\right)_{k}}$ как правило ${\displaystyle \mathbf {L} }$ в норме означает ${\displaystyle \left(\mathbf {L} _{k}\right)_{k}}$ как правило ${\displaystyle \mathbf {L} }$ входной. Это, в свою очередь, означает, что, поскольку каждый ${\displaystyle \mathbf {L} _{k}}$ является нижним треугольником с неотрицательными диагональными элементами, ${\displaystyle \mathbf {L} }$ это также.

Доказательство QR-разложением

Позволять ${\displaystyle \mathbf {A} }$ быть положительный полуопределенный Эрмитова матрица. Тогда его можно записать как произведение матрица квадратного корня, ${\displaystyle \mathbf {A} =\mathbf {B} \mathbf {B} ^{*}}$. Сейчас же QR-разложение может быть применен к ${\displaystyle \mathbf {B} ^{*}}$, в результате чего ${\displaystyle \mathbf {B} ^{*}=\mathbf {Q} \mathbf {R} }$, куда ${\displaystyle \mathbf {Q} }$ унитарен и ${\displaystyle \mathbf {R} }$ верхнетреугольный. Подставляя разложение в исходное равенство, получаем ${\displaystyle A=\mathbf {B} \mathbf {B} ^{*}=(\mathbf {QR} )^{*}\mathbf {QR} =\mathbf {R} ^{*}\mathbf {Q} ^{*}\mathbf {QR} =\mathbf {R} ^{*}\mathbf {R} }$. Параметр ${\displaystyle \mathbf {L} =\mathbf {R} ^{*}}$ завершает доказательство.

Обобщение

Факторизация Холецкого может быть обобщена^{[нужна цитата ]} в (не обязательно конечные) матрицы с операторными элементами. Позволять ${\displaystyle \{{\mathcal {H}}_{n}\}}$ быть последовательностью Гильбертовы пространства. Рассмотрим операторную матрицу

${\displaystyle \mathbf {A} ={\begin{bmatrix}\mathbf {A} _{11}&\mathbf {A} _{12}&\mathbf {A} _{13}&\;\\\mathbf {A} _{12}^{*}&\mathbf {A} _{22}&\mathbf {A} _{23}&\;\\\mathbf {A} _{13}^{*}&\mathbf {A} _{23}^{*}&\mathbf {A} _{33}&\;\\\;&\;&\;&\ddots \end{bmatrix}}}$

действуя на прямую сумму

${\displaystyle {\mathcal {H}}=\bigoplus _{n}{\mathcal {H}}_{n},}$

где каждый

${\displaystyle \mathbf {A} _{ij}:{\mathcal {H}}_{j}\rightarrow {\mathcal {H}}_{i}}$

это ограниченный оператор. Если А положительно (полуопределено) в том смысле, что для всех конечных k и для любого

${\displaystyle h\in \bigoplus _{n=1}^{k}{\mathcal {H}}_{k},}$

у нас есть ${\displaystyle \langle h,\mathbf {A} h\rangle \geq 0}$, то существует нижнетреугольная операторная матрица L такой, что А = LL*. Можно также взять диагональные записи L быть позитивным.

Реализации в библиотеках программирования

• Язык программирования C: the Научная библиотека GNU предоставляет несколько реализаций разложения Холецкого.
• Максима система компьютерной алгебры: функция холецкий вычисляет разложение Холецкого.
• GNU Octave Система численных вычислений предоставляет несколько функций для вычисления, обновления и применения разложения Холецкого.
• В ЛАПАК Библиотека обеспечивает высокопроизводительную реализацию разложения Холецкого, к которой можно получить доступ из Fortran, C и большинства языков.
• В Python, функция cholesky из модуля numpy.linalg выполняет разложение Холецкого.
• В Matlab и р, функция "chol" дает разложение Холецкого ..
• В Юля, функция "cholesky" из стандартной библиотеки LinearAlgebra дает разложение Холецкого.
• В Mathematica, функция "CholeskyDecomposition" может быть применена к матрице.
• В C ++, команда "чоль" из библиотеки броненосцев выполняет разложение Холецкого. В Библиотека Eigen предоставляет факторизации Холецкого как для разреженных, так и для плотных матриц.
• в КОРЕНЬ пакет доступен класс TDecompChol.
• В Аналитика, функция Decompose дает разложение Холецкого.
• В В библиотеке Apache Commons Math есть реализация который можно использовать в Java, Scala и любом другом языке JVM.

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

• Ранг цикла
• Неполная факторизация Холецкого
• Разложение матрицы
• Алгоритм минимальной степени
• Квадратный корень из матрицы
• Закон инерции Сильвестра
• Символическое разложение Холецкого

Примечания

1. Press, William H .; Саул А. Теукольский; Уильям Т. Веттерлинг; Брайан П. Фланнери (1992). Числовые рецепты на языке C: искусство научных вычислений (второе изд.). Кембриджский университет, Англия, EPress. п.994. ISBN  0-521-43108-5. Получено 2009-01-28.
2. Голуб и Ван Лоан (1996), п. 143), Хорн и Джонсон (1985), п. 407), Trefethen & Bau (1997)., п. 174).
3. Хорн и Джонсон (1985), п. 407).
4. «Матрицы - Диагонализация комплексной симметричной матрицы». MathOverflow. Получено 2020-01-25.
5. Шабауэр, Ханнес; Пачер, Кристоф; Сандерленд, Эндрю Дж .; Ганстерер, Уилфрид Н. (01.05.2010). «К параллельному решателю для обобщенных сложных симметричных задач на собственные значения». Процедуры информатики. ICCS 2010. 1 (1): 437–445. Дои:10.1016 / j.procs.2010.04.047. ISSN  1877-0509.
6. Голуб и Ван Лоан (1996), п. 147).
7. Нежный, Джеймс Э. (1998). Численная линейная алгебра для приложений в статистике. Springer. п. 94. ISBN  978-1-4612-0623-1.
8. Хайэм, Николас Дж. (1990). "Анализ разложения Холецкого полуопределенной матрицы". In Cox, M. G .; Хаммарлинг, С. Дж. (Ред.). Надежные численные вычисления. Оксфорд, Великобритания: Издательство Оксфордского университета. С. 161–185. ISBN  978-0-19-853564-5.
9. Кришнамурти, Аравинд; Менон, Дипак (2011). «Инверсия матриц с помощью разложения Холецкого». 1111: 4144. arXiv:1111.4144. Bibcode:2011arXiv1111.4144K.
10. Итак, Энтони Ман-Чо (2007). Подход полуопределенного программирования к проблеме реализации графа: теория, приложения и расширения (PDF) (Кандидат наук). Теорема 2.2.6.
11. Голуб и Ван Лоан (1996), Теорема 4.1.3)
12. Арора, Дж. Введение в оптимальный дизайн (2004), стр. 327. https://books.google.com/books?id=9FbwVe577xwC&pg=PA327
13. Документация Matlab randn. mathworks.com.
14. Фанг, Хаврен; О'Лири, Дайан П. (8 августа 2006 г.). «Модифицированные алгоритмы Холецкого: каталог с новыми подходами» (PDF).
15. Уоткинс, Д. (1991). Основы матричных вычислений. Нью-Йорк: Вили. п.84. ISBN  0-471-61414-9.
16. Нокедаль, Хорхе (2000). Численная оптимизация. Springer.
17. Фанг, Хаврен (24 августа 2007 г.). «Анализ блочных LDLT-факторизаций для симметричных неопределенных матриц».
18. По материалам: Стюарт, Г. В. (1998). Основные разложения. Филадельфия: Soc. по промышленной и прикладной математике. ISBN  0-89871-414-1.
19. Осборн, М. (2010), Приложение B.

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

• Дерениовский, Дариуш; Кубале, Марек (2004). «Факторизация Холецкого параллельных матриц и ранжирование графов». 5-я Международная конференция по параллельной обработке и прикладной математике (PDF). Конспект лекций по информатике. 3019. Springer-Verlag. С. 985–992. Дои:10.1007/978-3-540-24669-5_127. ISBN  978-3-540-21946-0. Архивировано из оригинал (PDF) на 2011-07-16.
• Голуб, Джин Х.; Ван Лоан, Чарльз Ф. (1996). Матричные вычисления (3-е изд.). Балтимор: Джонс Хопкинс. ISBN  978-0-8018-5414-9.
• Хорн, Роджер А .; Джонсон, Чарльз Р. (1985). Матричный анализ. Издательство Кембриджского университета. ISBN  0-521-38632-2.
• С. Дж. Юльер и Дж. К. Ульманн. "Общий метод аппроксимации нелинейных преобразований вероятностных распределений. ".
• С. Дж. Джулиер и Дж. К. Ульманн "Новое расширение фильтра Калмана на нелинейные системы ", в Proc. AeroSense: 11th Int. Symp. Aerospace / Defense Sensing, Simulation and Controls, 1997, pp. 182–193.
• Trefethen, Ллойд Н.; Бау, Дэвид (1997). Числовая линейная алгебра. Филадельфия: Общество промышленной и прикладной математики. ISBN  978-0-89871-361-9.
• Осборн, Майкл (2010). Байесовские гауссовские процессы для последовательного прогнозирования, оптимизации и квадратуры (PDF) (Тезис). Оксфордский университет.

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

История науки

• Sur la résolution numérique des systèmes d'équations linéaires, Рукопись Холецкого 1910 г., онлайн и проанализирована BibNum (на французском и английском языках) [для английского нажмите "Зарядное устройство"]

Информация

• «Факторизация Холецкого», Энциклопедия математики, EMS Press, 2001 [1994]
• «Холецкий разложение». PlanetMath.
• Разложение Холецкого, Краткое руководство по анализу данных
• Разложение Холецкого на www.math-linux.com
• Разложение Холецкого стало проще по науке меандертальца

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

• ЛАПАК представляет собой набор подпрограмм FORTRAN для решения задач плотной линейной алгебры
• АЛГЛИБ включает частичный перенос LAPACK на C ++, C #, Delphi, Visual Basic и т. д.
• libflame это библиотека C с функциональностью LAPACK.
• Заметки и видео по высокопроизводительной реализации факторизации Холецкого в Техасском университете в Остине.
• Холецкий: TBB + Threads + SSE это книга, объясняющая реализацию CF с TBB, потоками и SSE (на испанском языке).
• библиотека "Ceres Solver" от Google.
• Разложение ЛПНП подпрограммы в Matlab.
• Армадилло это пакет линейной алгебры C ++

Use of the matrix in simulation

• Generating Correlated Random Variables and Stochastic Processes, Martin Haugh, Колумбийский университет

Online calculators

• Online Matrix Calculator Performs Cholesky decomposition of matrices online.