stringtranslate.com

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

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

Заявление

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

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

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

Когда A — действительная матрица (следовательно, симметричная положительно определенная), факторизация может быть записана так, где L — действительная нижняя треугольная матрица с положительными диагональными элементами. [4] [5] [6]

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

Если эрмитова матрица A является только положительно полуопределенной, а не положительно определенной, то она все равно имеет разложение вида A = LL * , где диагональные элементы L могут быть равны нулю. [7] Разложение не обязательно должно быть уникальным, например: для любого θ . Однако, если ранг A равен r , то существует уникальный нижний треугольный L с ровно r положительными диагональными элементами и nr столбцами, содержащими все нули. [8]

В качестве альтернативы, разложение может быть сделано уникальным, когда выбор поворота фиксирован. Формально, если A является положительно полуопределенной матрицей n × n ранга r , то существует по крайней мере одна матрица перестановки P такая, что PAP T имеет уникальное разложение вида PAP T = LL * с , где L 1 является нижней треугольной матрицей r × r с положительной диагональю. [9]

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

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

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

По этой причине разложение LDL часто называют разложением Холецкого без квадратного корня . Для действительных матриц факторизация имеет вид A = LDL T и часто называется разложением LDLT (или разложением LDL T , или LDL′ ). Оно напоминает собственное разложение действительных симметричных матриц , A = QΛQ T , но на практике оно совсем другое, поскольку Λ и D не являются подобными матрицами .

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

Наоборот, учитывая классическое разложение Холецкого положительно определенной матрицы, если S — диагональная матрица, содержащая главную диагональ , то A можно разложить следующим образом (это изменяет масштаб каждого столбца, делая диагональные элементы равными 1),

Если A положительно определена , то все диагональные элементы D положительны. Для положительно полуопределенной A существует разложение , в котором число ненулевых элементов на диагонали D в точности равно рангу A. [11] Некоторые неопределенные матрицы, для которых не существует разложения Холецкого, имеют разложение LDL с отрицательными элементами в D : достаточно, чтобы первые n − 1 ведущих главных миноров A были невырожденными. [12]

Пример

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

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

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

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

Разложение Холецкого эквивалентно конкретному выбору сопряженных осей эллипсоида . [13] Подробно, пусть эллипсоид определен как , тогда по определению набор векторов является сопряженными осями эллипсоида тогда и только тогда . Тогда эллипсоид находится в точности там , где отображает базисный вектор , и является единичной сферой в n измерениях. То есть эллипсоид является линейным образом единичной сферы.

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

Разложение Холецкого соответствует выбору параллельности первой оси, нахождения в плоскости, охватываемой первыми двумя осями, и т. д. Это создает верхнетреугольную матрицу. Затем, есть , где нижнетреугольная.

Аналогично, анализ главных компонент соответствует выбору перпендикулярности. Тогда пусть и , и есть где — ортогональная матрица. Это тогда дает .

Приложения

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

Разложение Холецкого в основном используется для численного решения линейных уравнений . Если A симметрично и положительно определено, то можно решить, сначала вычислив разложение Холецкого , затем решив для y прямой подстановкой и, наконец, решив для x обратной подстановкой .

Альтернативный способ исключить извлечение квадратных корней при разложении — вычислить разложение ЛПНП , затем найти решение относительно y и, наконец, решить .

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

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

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

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

Нелинейные многомерные функции могут быть минимизированы по их параметрам с использованием вариантов метода Ньютона, называемых квазиньютоновскими методами. На итерации k поиск идет в направлении, определяемом решением для , где — направление шага, — градиент , а — приближение к матрице Гессе, сформированное путем повторения обновлений ранга 1 на каждой итерации. Две известные формулы обновления называются Дэвидон–Флетчер–Пауэлл (DFP) и Бройден–Флетчер–Гольдфарб–Шанно (BFGS). Потеря положительно определенного условия из-за ошибки округления исключается, если вместо обновления приближения к обратной матрице Гессе обновляется разложение Холецкого приближения самой матрицы Гессе. [14]

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

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

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

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

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

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

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

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

Вычисление

Существуют различные методы вычисления разложения Холецкого. Вычислительная сложность обычно используемых алгоритмов составляет O ( n 3 ) в общем случае. [ необходима цитата ] Все описанные ниже алгоритмы требуют около (1/3) n 3 FLOPs ( n 3 /6 умножений и столько же сложений) для действительных ароматов и (4/3) n 3 FLOPs для сложных ароматов, [16] где n — размер матрицы A . Следовательно, они имеют половину стоимости разложения LU , которое использует 2 n 3 /3 FLOPs (см. Trefethen and Bau 1997).

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

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

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

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

А (1)  := А .

На шаге i матрица A ( i ) имеет следующий вид: где I i −1 обозначает единичную матрицу размерности i − 1 .

Если матрица L i определяется как (обратите внимание, что a i,i > 0 , поскольку A ( i ) положительно определена), то A ( i ) можно записать как, где Обратите внимание, что b i b * i является внешним произведением , поэтому этот алгоритм называется версией внешнего произведения в (Голуб и Ван Лоан).

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

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

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

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

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

и поэтому следующие формулы для записей L :

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

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

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

для ( i = 0 ; i < dimensionSize ; i ++ ) { для ( j = 0 ; j <= i ; j ++ ) { float sum = 0 ; для ( k = 0 ; k < j ; k ++ ) sum += L [ i ][ k ] * L [ j ][ k ];                                   если ( i == j ) L [ i ][ j ] = sqrt ( A [ i ][ i ] - сумма ); иначе L [ i ][ j ] = ( 1,0 / L [ j ][ j ] * ( A [ i ][ j ] - сумма )); } }                   

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

сделать i = 1 , размер ( A , 1 ) L ( i , i ) = sqrt ( A ( i , i ) - скалярное_произведение ( 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 ) конец сделать                 

где conjgотносится к комплексному сопряжению элементов.

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

сделать i = 1 , размер ( A , 1 ) L ( i , i ) = sqrt ( A ( i , i ) - скалярное_произведение ( 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 ) конец сделать                 

где conjgотносится к комплексному сопряжению элементов.

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

Устойчивость вычислений

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Обновление рейтинга один

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

Вот функция [21], написанная на синтаксисе Matlab , которая реализует обновление ранга один:

функция  [L] = cholupdate ( L, x ) n = длина ( x ); для 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 ; если 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 ); конец конец конец                                                          

Обновление ранга n — это такое, при котором для матрицы обновляется разложение таким образом, что . Этого можно достичь, последовательно выполняя обновления ранга один для каждого из столбцов .

Понижение рейтинга на один

Понижение ранга на один похоже на обновление ранга на один, за исключением того, что сложение заменяется вычитанием: . Это работает только в том случае, если новая матрица по-прежнему положительно определена.

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

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

Если симметричная и положительно определенная матрица представлена ​​в блочной форме как

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

затем для новой матрицы , которая является такой же, но со вставкой новых строк и столбцов,

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

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

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

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

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

матрицы с удаленными строками и столбцами,

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

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

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

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

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

Если — положительно полуопределенная матрица , то последовательность состоит из положительно определенных матриц . (Это является непосредственным следствием, например, теоремы о спектральном отображении для полиномиального функционального исчисления.) Также в операторной норме . Из положительно определенного случая каждая имеет разложение Холецкого . По свойству операторной нормы,

Имеет место, поскольку снабженная операторной нормой, является алгеброй C*. Так же является ограниченным множеством в банаховом пространстве операторов, поэтому относительно компактным (потому что лежащее в основе векторное пространство конечномерно). Следовательно, у него есть сходящаяся подпоследовательность, также обозначаемая как , с пределом . Можно легко проверить, что это имеет требуемые свойства, т. е . , и является нижнетреугольным с неотрицательными диагональными элементами: для всех и ,

Следовательно, . Поскольку базовое векторное пространство конечномерно, все топологии на пространстве операторов эквивалентны. Так что стремится к по норме, значит стремится к по элементам. Это, в свою очередь, подразумевает, что, поскольку каждый является нижним треугольным с неотрицательными диагональными элементами, также является.

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

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

Обобщение

Факторизация Холецкого может быть обобщена [ требуется ссылка ] на (не обязательно конечные) матрицы с операторными записями. Пусть будет последовательностью гильбертовых пространств . Рассмотрим операторную матрицу

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

где каждый

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

существует , то существует нижняя треугольная операторная матрица L такая, что A = LL * . Можно также считать диагональные элементы L положительными.

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

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

Примечания

  1. ^ Бенуа (1924). «Примечание к методу разрешения нормальных уравнений, происшедшему с применением метода moindres carrés à un système d'équations linéaires en nombre inférieur à celui des inconnues (Procédé du Commandant Cholesky)». Бюллетень Géodesique (на французском языке). 2 : 66–67. дои : 10.1007/BF03031308.
  2. ^ ab Press, William H.; Saul A. Teukolsky; William T. Vetterling; Brian P. Flannery (1992). Numerical Recipes in C: The Art of Scientific Computing (второе изд.). Cambridge University England EPress. стр. 994. ISBN 0-521-43108-5. Получено 28.01.2009 .
  3. ^ Голуб и Ван Лоан (1996, стр. 143), Хорн и Джонсон (1985, стр. 407), Трефетен и Бау (1997, стр. 174).
  4. ^ Хорн и Джонсон (1985, стр. 407).
  5. ^ "матрицы - Диагонализация комплексной симметричной матрицы". MathOverflow . Получено 2020-01-25 .
  6. ^ Шабауэр, Ханнес; Пахер, Кристоф; Сандерленд, Эндрю Г.; Ганстерер, Вильфрид Н. (2010-05-01). «К параллельному решателю для обобщенных комплексных симметричных задач на собственные значения». Procedia Computer Science . ICCS 2010. 1 (1): 437–445. doi : 10.1016/j.procs.2010.04.047 . ISSN  1877-0509.
  7. ^ Голуб и Ван Лоан (1996, стр. 147).
  8. ^ Джентл, Джеймс Э. (1998). Численная линейная алгебра для приложений в статистике . Springer. стр. 94. ISBN 978-1-4612-0623-1.
  9. ^ Хайэм, Николас Дж. (1990). «Анализ разложения Холецкого полуопределенной матрицы». В Коксе, МГ; Хаммарлинге, СДж (ред.). Надежные численные вычисления . Оксфорд, Великобритания: Oxford University Press. стр. 161–185. ISBN 978-0-19-853564-5.
  10. ^ ab Krishnamoorthy, Aravindh; Menon, Deepak (2011). «Обращение матрицы с использованием разложения Холецкого». 1111 : 4144. arXiv : 1111.4144 . Bibcode :2011arXiv1111.4144K. {{cite journal}}: Цитировать журнал требует |journal=( помощь )
  11. ^ Так, Энтони Ман-Чо (2007). Подход полуопределенного программирования к проблеме реализации графа: теория, приложения и расширения (PDF) (PhD). Теорема 2.2.6.
  12. ^ Голуб и Ван Лоан (1996, Теорема 4.1.3)
  13. ^ Поуп, Стивен Б. «Алгоритмы для эллипсоидов». Отчет Корнелльского университета № FDA (2008): 08-01.
  14. ^ Арора, Джасбир Сингх (2004-06-02). Введение в оптимальное проектирование. Elsevier. ISBN 978-0-08-047025-2.
  15. ^ Документация Matlab randn. mathworks.com.
  16. ^ ?potrf Библиотека Intel® Math Kernel [1]
  17. ^ Фанг, Хоу-рен; О'Лири, Дайан П. (8 августа 2006 г.). «Модифицированные алгоритмы Холецкого: каталог с новыми подходами» (PDF) . {{cite journal}}: Цитировать журнал требует |journal=( помощь )
  18. ^ Уоткинс, Д. (1991). Основы матричных вычислений . Нью-Йорк: Wiley. стр. 84. ISBN 0-471-61414-9.
  19. ^ Нокедаль, Хорхе (2000). Численная оптимизация . Springer.
  20. ^ Фанг, Хоу-рен (24 августа 2007 г.). «Анализ блочных факторизаций LDLT для симметричных неопределенных матриц». {{cite journal}}: Цитировать журнал требует |journal=( помощь )
  21. ^ На основе: Stewart, GW (1998). Базовые разложения . Филадельфия: Soc. for Industrial and Applied Mathematics. ISBN 0-89871-414-1.
  22. ^ Осборн, М. (2010), Приложение B.

Ссылки

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

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

Информация

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

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

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