В числовой линейной алгебре , то Якоби собственного значение алгоритм является итерационным методом для вычисления собственных значений и собственных векторов одного реальной симметричной матрицы (процесс , известный как диагонализация ). Он назван в честь Карла Густава Якоба Якоби , который первым предложил метод в 1846 году [1], но стал широко использоваться только в 1950-х годах с появлением компьютеров. [2]
Описание
Позволять - симметричная матрица, а - матрица вращения Гивенса . Потом:
симметричен и похож на.
Более того, есть записи:
где а также .
С ортогонален, а также имеют ту же норму Фробениуса (квадратный корень из суммы квадратов всех компонентов), однако мы можем выбрать такой, что , в таком случае имеет большую сумму квадратов по диагонали:
Установите его равным 0 и переставьте:
если
Чтобы оптимизировать этот эффект, S ij должен быть недиагональным элементом с наибольшим абсолютным значением, называемым точкой поворота .
Метод собственных значений Якоби многократно выполняет вращения, пока матрица не станет почти диагональной. Тогда элементы в диагонали являются приближениями (вещественных) собственных значений S .
Конвергенция
Если является стержневым элементом, то по определению для . Позволять обозначают сумму квадратов всех недиагональных элементов . С точно недиагональные элементы, имеем или же . Сейчас. Из этого следует или же , т.е. последовательность поворотов Якоби сходится хотя бы линейно с коэффициентом к диагональной матрице.
Номер Вращение Якоби называется разверткой; позволятьобозначим результат. Предыдущая оценка дает
- ,
т.е. последовательность разверток сходится хотя бы линейно с множителем ≈ .
Однако следующий результат Шёнхаге [3] дает локально квадратичную сходимость. Для этого пусть S имеет m различных собственных значений с кратностями и пусть d > 0 - наименьшее расстояние между двумя разными собственными значениями. Позвольте нам позвонить в ряд
Якоби вращает развертку Шенхаге. Если обозначает результат, тогда
- .
Таким образом, сходимость становится квадратичной, как только
Расходы
Каждое вращение Якоби может быть выполнено за O ( n ) шагов, если известен поворотный элемент p . Однако поиск p требует проверки всех N ≈ 1/2 n 2 недиагональных элемента. Мы также можем уменьшить эту сложность до O ( n ), если введем дополнительный массив индексов со свойством, что является показателем наибольшего элемента в строке I , ( я = 1, ..., п - 1) текущего S . Тогда индексы стержня ( k , l ) должны быть одной из пар. Также обновление индексного массива может быть выполнено в среднем за O ( n ) сложность: во-первых, максимальная запись в обновленных строках k и l может быть найдена за O ( n ) шагов. В других строках i меняются только записи в столбцах k и l . Перебирая эти строки, еслине является ни k, ни l , достаточно сравнить старый максимум на к новым записям и обновлению если необходимо. Еслидолжен быть равен k или l, а соответствующая запись уменьшилась во время обновления, максимум по строке i должен быть найден с нуля за O ( n ) сложность. Однако в среднем это происходит только один раз за оборот. Таким образом, каждое вращение имеет сложность в среднем O ( n ) и одну развертку O ( n 3 ), что эквивалентно одному умножению матриц. Кроме того,должен быть инициализирован до начала процесса, что может быть выполнено за n 2 шагов.
Обычно метод Якоби сходится с числовой точностью после небольшого количества разверток. Обратите внимание, что несколько собственных значений сокращают количество итераций, поскольку.
Алгоритм
Следующий алгоритм представляет собой описание метода Якоби в математической нотации. Он вычисляет вектор e, который содержит собственные значения, и матрицу E, которая содержит соответствующие собственные векторы, т. Е. - собственное значение, а столбец ортонормированный собственный вектор для , я = 1, ..., п .
процедура jacobi ( S ∈ R n × n ; out e ∈ R n ; out E ∈ R n × n ) var i , k , l , m , state ∈ N s , c , t , p , y , d , r ∈ R ind ∈ N n изменено ∈ L n функция maxind ( k ∈ N ) ∈ N ! индекс самого большого недиагонального элемента в строке k m : = k +1 для i : = k +2 to n do if │ S ki │> │ S km │ then m : = i endif endfor return m endfunc обновление процедуры ( k ∈ N ; t ∈ R )! обновить e k и его статус y : = e k ; e k : = y + t, если изменилось k и ( y = e k ), то изменилось k : = false; state : = state −1 elsif ( k не изменено ) и ( y ≠ e k ) затем изменилось k : = true; состояние : = состояние +1 endif endproc процедура rotate ( k , l , i , j ∈ N )! выполнять вращение S Ij , S кл ┌ ┐ ┌ ┐┌ ┐ │ S kl │ │ c - s ││ S kl │ │ │: = │ ││ │ │ S ij │ │ s c ││ S ij │ └ ┘ └ ┘└ ┘ ENDPROC ! init e, E и массивы ind, изменено E : = I ; state : = n для k : = от 1 до n do ind k : = maxind ( k ); e k : = S kk ; изменено k : = true endfor, а состояние ≠ 0 do ! следующий оборот m : = 1! найти индекс (K, L) поворотного р для к : = 2 до п -1 делать , если │ S к Ind к │> │ S м Ind м │ затем м : = к ENDIF ENDFOR к : = м ; l : = ind m ; p : = S kl ! вычислить c = cos φ, s = sin φ y : = ( e l - e k ) / 2; d : = │ y │ + √ ( p 2 + y 2 ) r : = √ ( p 2 + d 2 ); c : = d / r ; s : = p / r ; t : = p 2 / d, если y <0, то s : = - s ; t : = - t endif S kl : = 0,0; обновить ( k , - t ); обновление ( l , t ) ! Поворот строк и столбцов к и л для I : = 1 до K -1 сделать поворот ( я , к , я , л ) ENDFOR для I : = K +1 к л -1 сделать поворот ( к , я , я , л ) ENDFOR для я : = л + 1 , чтобы п сделать поворот ( к , я , л , я ) ENDFOR ! Поворот собственных векторов для I : = 1 до п делать ┌ ┐ ┌ ┐┌ ┐ │ E ik │ │ c - s ││ E ik │ │ │: = │ ││ │ │ E il │ │ s c ││ E il │ └ ┘ └ ┘└ ┘ ENDFOR ! строки k, l изменились, обновить строки ind k , ind l ind k : = maxind ( k ); ind l : = maxind ( l ) конец цикла
Заметки
1. Измененный логический массив содержит статус каждого собственного значения. Если числовое значение или же изменяется во время итерации, соответствующий компонент changed устанавливается в значение true , в противном случае - в false . Целочисленное состояние подсчитывает количество измененных компонентов, имеющих значение true . Итерация прекращается, как только state = 0. Это означает, что ни одно из приближенийнедавно изменил свое значение, и поэтому маловероятно, что это произойдет, если итерация продолжится. Здесь предполагается, что операции с плавающей запятой оптимально округляются до ближайшего числа с плавающей запятой.
2. Верхний треугольник матрицы S разрушается, а нижний треугольник и диагональ не изменяются. Таким образом, при необходимости можно восстановить S согласно
для k : = от 1 до n −1 делать ! восстановить матрицу S для l : = k +1 до n do S kl : = S lk endfor endfor
3. Собственные значения не обязательно в порядке убывания. Этого можно добиться с помощью простого алгоритма сортировки.
for k : = 1 to n −1 do m : = k for l : = k +1 to n do if e l > e m then m : = l endif endfor if k ≠ m then swap e m , e k swap E m , E k endif endfor
4. Алгоритм написан с использованием матричной записи (массивы на основе 1 вместо 0).
5. При реализации алгоритма часть, указанная в матричной записи, должна выполняться одновременно.
6. Эта реализация неправильно учитывает случай, когда одно измерение является независимым подпространством. Например, если задана диагональная матрица, вышеуказанная реализация никогда не завершится, поскольку ни одно из собственных значений не изменится. Следовательно, в реальных реализациях необходимо добавить дополнительную логику для учета этого случая.
Пример
Позволять
Затем jacobi производит следующие собственные значения и собственные векторы после 3 циклов (19 итераций):
Приложения для вещественных симметричных матриц
Когда собственные значения (и собственные векторы) симметричной матрицы известны, легко вычисляются следующие значения.
- Особые значения
- Сингулярные значения (квадратной) матрицы A являются квадратными корнями из (неотрицательных) собственных значений матрицы . В случае симметричной матрицы S имеем , следовательно, сингулярные значения S являются модулями собственных значений S
- 2-норма и спектральный радиус
- 2-норма матрицы A - это норма, основанная на евклидовой векторной норме, т. Е. Наибольшее значение когда x пробегает все векторы с . Это самое большое сингулярное значение A . В случае симметричной матрицы это наибольшее абсолютное значение ее собственных векторов и, следовательно, равно ее спектральному радиусу.
- Номер условия
- Число обусловленности невырожденной матрицы A определяется как . В случае симметричной матрицы это модуль отношения наибольшего и наименьшего собственного значения. Матрицы с большими числами обусловленности могут привести к численно нестабильным результатам: небольшое возмущение может привести к большим ошибкам. Матрицы Гильберта - самые известные плохо обусловленные матрицы. Например, матрица Гильберта четвертого порядка имеет условие 15514, а для порядка 8 - 2,7 × 10 8 .
- Классифицировать
- Матрица A имеет ранг r, если она имеет r столбцов, которые линейно независимы, а остальные столбцы линейно зависят от них. Эквивалентно, г есть размерность диапазона A . Кроме того, это количество ненулевых особых значений.
- В случае симметричной матрицы r - количество ненулевых собственных значений. К сожалению, из-за ошибок округления числовые приближения нулевых собственных значений могут не быть нулевыми (также может случиться, что численное приближение равно нулю, а истинное значение - нет). Таким образом, можно вычислить числовой ранг, только приняв решение, какое из собственных значений достаточно близко к нулю.
- Псевдообратный
- Псевдообратная матрица A - это единственная матрица для которых AX и XA симметричны и для которых AXA = A, XAX = X выполняется. Если A неособое, то ' .
- Когда вызывается процедура jacobi (S, e, E), то соотношение где Diag ( e ) обозначает диагональную матрицу с вектором e на диагонали. Позволять обозначим вектор, где заменяется на если и на 0, если равен (численно близок к) нулю. Поскольку матрица E ортогональна, отсюда следует, что псевдообратная матрица S имеет вид .
- Решение методом наименьших квадратов
- Если матрица A не имеет полного ранга, решения линейной системы Ax = b может не быть . Однако можно искать вектор x, для которого минимально. Решение . В случае симметричной матрицы S, как и раньше, имеем .
- Матрица экспоненциальная
- Из можно найти где exp e - вектор, где заменяется на . Таким же образом очевидным образом вычисляется f ( S ) для любой (аналитической) функции f .
- Линейные дифференциальные уравнения
- Дифференциальное уравнение x ' = Ax , x (0) = a имеет решение x ( t ) = exp ( t A ) a . Для симметричной матрицы S следует, что . Если является разложением a по собственным векторам S , то .
- Позволять - векторное пространство, натянутое на собственные векторы матрицы S, соответствующие отрицательному собственному значению, и аналогично для положительных собственных значений. Если тогда т.е. точка равновесия 0 притягивает x ( t ). Если тогда , т.е. 0 отталкивает x ( t ). а также называются стабильными и неустойчивые многообразия для S . Если a имеет компоненты в обоих многообразиях, то один компонент притягивается, а один компонент отталкивается. Следовательно, x ( t ) приближается в виде .
Обобщения
Метод Якоби был обобщен на комплексные эрмитовы матрицы , общие несимметричные вещественные и комплексные матрицы, а также на блочные матрицы.
Поскольку сингулярные значения вещественной матрицы являются квадратными корнями из собственных значений симметричной матрицы его также можно использовать для расчета этих значений. В этом случае метод модифицирован таким образом, что S не следует явно вычислять, что снижает опасность ошибок округления . Обратите внимание, что с участием .
Метод Якоби также хорошо подходит для параллелизма.
Рекомендации
- ^ Якоби, CGJ (1846). "Über ein leichtes Verfahren, die in der Theorie der Säkularstörungen vorkommenden Gleichungen numerisch aufzulösen" . Журнал Crelle (на немецком языке). 1846 (30): 51–94. DOI : 10.1515 / crll.1846.30.51 .
- ^ Голуб, Г.Х .; ван дер Ворст, HA (2000). «Вычисление собственных значений в 20 веке» . Журнал вычислительной и прикладной математики . 123 (1–2): 35–65. DOI : 10.1016 / S0377-0427 (00) 00413-1 .
- ^ Шёнхаге, А. (1964). "Zur quadratischen Konvergenz des Jacobi-Verfahrens". Numerische Mathematik (на немецком языке). 6 (1): 410–412. DOI : 10.1007 / BF01386091 . Руководство по ремонту 0174171 .
дальнейшее чтение
- Нажмите, WH; Теукольский, С.А. Феттерлинг, Вашингтон; Фланнери, Б.П. (2007), «Раздел 11.1. Преобразования Якоби симметричной матрицы» , Численные рецепты: Искусство научных вычислений (3-е изд.), Нью-Йорк: Cambridge University Press, ISBN 978-0-521-88068-8
- Рутисхаузер, Х. (1966). «Серия справочников по линейной алгебре: метод Якоби для вещественных симметричных матриц». Numerische Mathematik . 9 (1): 1–10. DOI : 10.1007 / BF02165223 . Руководство по ремонту 1553948 .
- Самех, AH (1971). «Об алгоритмах Якоби и подобных Якоби для параллельного компьютера» . Математика вычислений . 25 (115): 579–590. DOI : 10.1090 / s0025-5718-1971-0297131-6 . JSTOR 2005221 . Руководство по ремонту 0297131 .
- Шрофф, Гаутам М. (1991). «Параллельный алгоритм для собственных значений и собственных векторов общей комплексной матрицы». Numerische Mathematik . 58 (1): 779–805. CiteSeerX 10.1.1.134.3566 . DOI : 10.1007 / BF01385654 . Руководство по ремонту 1098865 .
- Веселич, К. (1979). «Об одном классе процедур типа Якоби для диагонализации произвольных вещественных матриц». Numerische Mathematik . 33 (2): 157–172. DOI : 10.1007 / BF01399551 . Руководство по ремонту 0549446 .
- Веселич, К .; Венцель, HJ (1979). «Квадратично сходящийся метод Якоби для вещественных матриц с комплексными собственными значениями». Numerische Mathematik . 33 (4): 425–435. DOI : 10.1007 / BF01399324 . Руководство по ремонту 0553351 .
Внешние ссылки
- Реализация в Matlab алгоритма Якоби, избегающего тригонометрических функций
- Реализация C ++ 11