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

В линейной алгебре или разложение Холецкого факторизация Холецкого (произносится / ʃ ə ˈ l ɛ s k i / shə- LES -kee ) — это разложение эрмитовой определенной положительно матрицы в произведение нижней треугольной матрицы и ее сопряженной transpose , что полезно для эффективных численных решений, например, моделирования Монте-Карло . Он был открыт Андре-Луи Холеским для реальных матриц и посмертно опубликован в 1924 году. Когда это применимо, разложение Холецкого примерно в два раза эффективнее LU-разложения для решения систем линейных уравнений .

Заявление

Разложение Холецкого эрмитовой положительно определенной матрицы A представляет собой разложение вида

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

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

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

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

$\mathbf {A} =\mathbf {LL} ^{\mathsf {T}},$

где L — действительная нижне-треугольная матрица с положительными диагональными элементами.

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

Если эрмитова матрица A является только положительно полуопределенной, а не положительно определенной, то она все равно имеет разложение формы A = LL * , где диагональные элементы L могут быть равны нулю. Декомпозиция не обязательно должна быть уникальной, например:

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

для любого θ . Однако если ранг A равен r , то существует единственный нижний треугольник L с ровно r положительными диагональными элементами и n - r столбцами, содержащими все нули.

Альтернативно, декомпозиция может быть сделана уникальной, если фиксирован поворотный выбор. Формально, если A — положительно-полуопределенная матрица размера n × n ранга r , то существует хотя бы одна матрица перестановок P такая, что PAP ^Т имеет уникальное разложение вида PAP ^Т = ЛЛ ^* с $\mathbf {L} ={\begin{bmatrix}\mathbf {L} _{1}&0\\\mathbf {L} _{2}&0\end{bmatrix}}$, где L ₁ — нижняя треугольная матрица размера r × r с положительной диагональю.

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

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

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

где L — младшая единичная треугольная (однотреугольная) матрица, а D — диагональная матрица. То есть диагональные элементы L должны быть равны 1 за счет введения дополнительной диагональной матрицы D в разложение. Основное преимущество заключается в том, что разложение ЛПНП можно вычислять и использовать по существу с помощью тех же алгоритмов, но без извлечения квадратных корней.

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

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

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

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

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

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

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

Пример

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

${\begin{aligned}{\begin{pmatrix}4&12&-16\\12&37&-43\\-16&-43&98\\\end{pmatrix}}={\begin{pmatrix}2&0&0\\6&1&0\\-8&5&3\\\end{pmatrix}}{\begin{pmatrix}2&6&-8\\0&1&5\\0&0&3\\\end{pmatrix}}.\end{aligned}}$

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

${\begin{aligned}{\begin{pmatrix}4&12&-16\\12&37&-43\\-16&-43&98\\\end{pmatrix}}&={\begin{pmatrix}1&0&0\\3&1&0\\-4&5&1\\\end{pmatrix}}{\begin{pmatrix}4&0&0\\0&1&0\\0&0&9\\\end{pmatrix}}{\begin{pmatrix}1&3&-4\\0&1&5\\0&0&1\\\end{pmatrix}}.\end{aligned}}$

Геометрическая интерпретация

См. Также: Преобразование отбеливания.

Эллипс — это линейное изображение единичного круга. Два вектора $v_{1},v_{2}$ — сопряженные оси эллипса, выбранные так, что $v_{1}$ быть параллельным первой оси, $v_{2}$ находиться внутри плоскости, охватываемой первыми двумя осями.

Разложение Холецкого эквивалентно определенному выбору осей эллипсоида сопряженных . Подробно, пусть эллипсоид определяется как $y^{T}Ay=1$, то по определению набор векторов $v_{1},...,v_{n}$ являются сопряженными осями эллипсоида тогда и только тогда, когда $v_{i}^{T}Av_{j}=\delta _{ij}$. Тогда эллипсоид в точности $$\left\{\sum _{i}x_{i}v_{i}:x^{T}x=1\right\}=f(\mathbb {S} ^{n})$$где $f$ отображает базисный вектор $e_{i}\mapsto v_{i}$, и $\mathbb {S} ^{n}$ - единичная сфера в n измерениях. То есть эллипсоид представляет собой линейное изображение единичной сферы.

Определите матрицу $V:=[v_{1}|v_{2}|\cdots |v_{n}]$, затем $v_{i}^{T}Av_{j}=\delta _{ij}$ эквивалентно $V^{T}AV=I$. Разный выбор сопряженных осей соответствует разным разложениям.

Разложение Холецкого соответствует выбору $v_{1}$ быть параллельным первой оси, $v_{2}$ находиться внутри плоскости, охватываемой первыми двумя осями, и так далее. Это делает $V$ верхнетреугольная матрица. Тогда есть $A=LL^{T}$, где $L=(V^{-1})^{T}$ является нижне-треугольным.

Аналогично, анализ главных компонент соответствует выбору $v_{1},...,v_{n}$ быть перпендикулярным. Тогда пусть $\lambda =1/\|v_{i}\|^{2}$ и $\Sigma =\mathrm {diag} (\lambda _{1},...,\lambda _{n})$, и есть $V=U\Sigma ^{-1/2}$ где $U$ является ортогональной матрицей. Тогда это дает $A=U\Sigma U^{T}$.

Приложения

Численное решение системы линейных уравнений

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

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

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

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

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

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

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

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

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

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

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

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

Инверсия матрицы

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

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

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

Вычисление

Существуют различные методы расчета разложения Холецкого. Вычислительная сложность часто используемых алгоритмов составляет O ( n ³ ) в общем. Все алгоритмы, описанные ниже, включают около (1/3) n ³ ФЛОПЫ ( n ³ /6 умножений и столько же сложений) для реальных ароматов и (4/3) n ³ FLOP для сложных вкусов, где n размер матрицы A. — Следовательно, у них есть половина стоимости LU-разложения , которое использует 2 n ³ /3 флопа (см. Трефетен и Бау, 1997).

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

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

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

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

А ⁽¹⁾ : А. =

На шаге i матрица A ^{( я )} имеет следующую форму:

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

где I _{i −1} обозначает единичную матрицу размерности i − 1.

Если матрица L _i определяется формулой

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

(заметим, что a _i,i > 0, поскольку A ^{( я )} положительно определен), тогда А ^{( я )} можно записать как

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

где

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

Обратите внимание, что b _i b * _i — внешний продукт этот алгоритм называется версией внешнего продукта , поэтому в (Golub & Van Loan) .

Это повторяется для i от 1 до n . После n шагов A ^{( п +1)} = I , и, следовательно, искомая нижняя треугольная матрица L вычисляется как

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

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

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

Если уравнение

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

выписано, получается следующее:

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

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

$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.$

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

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

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

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

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

• Алгоритм Холеского-Банакевича начинается с верхнего левого угла матрицы L и продолжает вычисление матрицы построчно.

for (i = 0; i < dimensionSize; i++) {
    for (j = 0; j <= i; j++) {
        float sum = 0;
        for (k = 0; k < j; k++)
            sum += L[i][k] * L[j][k];

        if (i == j)
            L[i][j] = sqrt(A[i][i] - sum);
        else
            L[i][j] = (1.0 / L[j][j] * (A[i][j] - sum));
    }
}

Приведенный выше алгоритм можно кратко выразить как объединение скалярного произведения и умножения матриц в векторизованных языках программирования, таких как Фортран, следующим образом:

do i = 1, size(A,1)
    L(i,i) = sqrt(A(i,i) - dot_product(L(i,1:i-1), L(i,1:i-1)))
    L(i+1:,i) = (A(i+1:,i) - matmul(conjg(L(i,1:i-1)), L(i+1:,1:i-1))) / L(i,i)
end do

где conjg относится к комплексно-сопряженным элементам.

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

  for (j = 0; j < dimensionSize; j++) {
      float sum = 0;
      for (k = 0; k < j; k++) {
          sum += L[j][k] * L[j][k];
      }
      L[j][j] = sqrt(A[j][j] - sum);
  
      for (i = j + 1; i < dimensionSize; i++) {
          sum = 0;
          for (k = 0; k < j; k++) {
              sum += L[i][k] * L[j][k];
          }
          L[i][j] = (1.0 / L[j][j] * (A[i][j] - sum));
      }
  }
  

Приведенный выше алгоритм можно кратко выразить как объединение скалярного произведения и умножения матриц в векторизованных языках программирования, таких как Фортран, следующим образом:

do i = 1, size(A,1)
    L(i,i) = sqrt(A(i,i) - dot_product(L(1:i-1,i), L(1:i-1,i)))
    L(i,i+1:) = (A(i,i+1:) - matmul(conjg(L(1:i-1,i)), L(1:i-1,i+1:))) / L(i,i)
end do

где conjg относится к комплексно-сопряженным элементам.

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

Стабильность вычислений

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

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

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

Здесь ||·|| ₂ — матричная 2-норма , c _n — небольшая константа, зависящая от n , а ε обозначает единичное округление .

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

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

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

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

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

$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 вещественны, если A вещественно.

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

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

$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.$

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

Блочный вариант

* при использовании с неопределенными матрицами Известно, что факторизация LDL без тщательного поворота нестабильна; в частности, элементы факторизации могут расти произвольно. Возможным улучшением является выполнение факторизации блочных подматриц, обычно 2 × 2:

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

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

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

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

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

Обновление декомпозиции

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

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

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

Вот функция написанный в синтаксисе Matlab , который реализует обновление первого ранга:

function [L] = cholupdate(L, x)
    n = length(x);
    for k = 1:n
        r = sqrt(L(k, k)^2 + x(k)^2);
        c = r / L(k, k);
        s = x(k) / L(k, k);
        L(k, k) = r;
        if k < n
            L((k+1):n, k) = (L((k+1):n, k) + s * x((k+1):n)) / c;
            x((k+1):n) = c * x((k+1):n) - s * L((k+1):n, k);
        end
    end
end

Обновление ранга n — это обновление матрицы $\mathbf {M}$ обновляют разложение так, что ${\tilde {\mathbf {A} }}=\mathbf {A} +\mathbf {M} \mathbf {M} ^{*}$. Этого можно достичь, последовательно выполняя обновления первого ранга для каждого из столбцов таблицы. $\mathbf {M}$.

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

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

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

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

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

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

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

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

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

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

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

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

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

${\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}&=\mathrm {chol} \left(\mathbf {A} _{22}-\mathbf {S} _{12}^{\mathrm {T} }\mathbf {S} _{12}\right),\\\mathbf {S} _{23}&=\mathbf {S} _{22}^{\mathrm {T} }\setminus \left(\mathbf {A} _{23}-\mathbf {S} _{12}^{\mathrm {T} }\mathbf {S} _{13}\right),\\\mathbf {S} _{33}&=\mathrm {chol} \left(\mathbf {L} _{33}^{\mathrm {T} }\mathbf {L} _{33}-\mathbf {S} _{23}^{\mathrm {T} }\mathbf {S} _{23}\right).\end{aligned}}$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

$\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 \,.$

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

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

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

Обобщение

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

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

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

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

где каждый

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

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

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

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

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

• Язык программирования C : Научная библиотека GNU предоставляет несколько реализаций разложения Холецкого.
• Система компьютерной алгебры Maxima : функция cholesky вычисляет разложение Холецкого.
• Система численных вычислений GNU Octave предоставляет несколько функций для расчета, обновления и применения разложения Холецкого.
• Библиотека LAPACK обеспечивает высокопроизводительную реализацию разложения Холецкого, доступ к которой возможен из Fortran , C и большинства языков.
• В Python функция cholesky из numpy.linalg модуль выполняет разложение Холецкого.
• В Матлабе ​ chol функция дает разложение Холецкого. Обратите внимание, что chol по умолчанию использует верхний треугольный коэффициент входной матрицы, т.е. он вычисляет $A=R^{*}R$ где $R$ имеет верхнюю треугольную форму. Вместо этого можно передать флаг, чтобы использовать нижний треугольный коэффициент.
• В Р , chol функция дает разложение Холецкого.
• В Юлии , cholesky функция от LinearAlgebra стандартная библиотека дает разложение Холецкого.
• В системе Mathematica функция « CholeskyDecomposition" можно применить к матрице.
• В C++ несколько библиотек линейной алгебры поддерживают это разложение:
  • Armadillo (библиотека C++) предоставляет команду chol выполнить разложение Холецкого.
  • Библиотека Эйгена предоставляет факторизации Холецкого как для разреженных, так и для плотных матриц.
  • В ROOT- пакете TDecompChol класс доступен.

• В Аналитике функция Decompose дает разложение Холецкого.
• Библиотека Apache Commons Math имеет реализацию , которую можно использовать в Java, Scala и любом другом языке JVM.

См. также

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

Ссылки

• Деренёвский, Дариуш; Кубале, Марек (2004). «Параллельная факторизация матриц Холецкого и ранжирование графов». 5-я Международная конференция по параллельной обработке и прикладной математике (PDF) . Конспекты лекций по информатике. Том. 3019. Шпрингер-Верлаг. стр. 985–992. дои : 10.1007/978-3-540-24669-5_127 . ISBN  978-3-540-21946-0 . Архивировано из оригинала (PDF) 16 июля 2011 г.
• Голуб, Джин Х .; Ван Лоан, Чарльз Ф. (1996). Матричные вычисления (3-е изд.). Балтимор: Джонс Хопкинс. ISBN  978-0-8018-5414-9 .
• Хорн, Роджер А.; Джонсон, Чарльз Р. (1985). Матричный анализ . Издательство Кембриджского университета. ISBN  0-521-38632-2 .
• С. Дж. Жюльер и Дж. К. Ульман. « Общий метод аппроксимации нелинейных преобразований вероятностных распределений ».
• С. Дж. Жюльер и Дж. К. Ульманн, « Новое расширение фильтра Калмана для нелинейных систем », в Proc. AeroSense: 11-й Международный. Симп. Аэрокосмическое/оборонное зондирование, моделирование и управление, 1997, стр. 182–193.
• Трефетен, Ллойд Н .; Бау, Дэвид (1997). Численная линейная алгебра . Филадельфия: Общество промышленной и прикладной математики. ISBN  978-0-89871-361-9 .
• Осборн, Майкл (2010). Байесовские гауссовские процессы для последовательного прогнозирования, оптимизации и квадратуры (PDF) (диссертация). Оксфордский университет.
• Рушель, Жоау Паулу Тараскони, степень бакалавра « Параллельные реализации разложения Холецкого на процессорах и графических процессорах », Федеральный университет Риу-Гранди-ду-Сул, Институт информатики, 2016, стр. 29-30.

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

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

• О численном разрешении систем линейных уравнений , рукопись Холеского 1910 года, онлайн и проанализированная на BibNum (на французском и английском языках) [для английского языка нажмите «Для загрузки»]

Информация

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

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

• LAPACK — это набор подпрограмм FORTRAN для решения задач плотной линейной алгебры (DPOTRF, DPOTRF2, подробная информация о производительности ).
• ALGLIB включает частичный порт LAPACK на C++, C#, Delphi, Visual Basic и т. д. (spdmatrixcholesky, hpdmatrixcholesky).
• libflame — это библиотека C с функциональностью LAPACK.
• Заметки и видео о высокопроизводительной реализации факторизации Холецкого в Техасском университете в Остине.
• Cholesky: TBB + Threads + SSE — книга, объясняющая реализацию CF с TBB, потоками и SSE (на испанском языке).
• библиотека «Ceres Solver» от Google.
• Процедуры разложения ЛПНП в Matlab.
• Armadillo — пакет линейной алгебры C++.
• Rosetta Code — сайт по программированию хрестоматии. по теме страницы .
• AlgoWiki — открытая энциклопедия свойств алгоритмов и особенностей их реализации по теме страницы.
• Intel® oneAPI Math Kernel Library Оптимизированная Intel математическая библиотека для числовых вычислений ?potrf , ?potrs

Использование матрицы в моделировании

• Генерация коррелированных случайных величин и случайных процессов , Мартин Хо, Колумбийский университет

Онлайн калькуляторы

• Онлайн-калькулятор матриц Выполняет разложение матриц Холецкого онлайн.