Алгоритм собственных значений Якоби

В числовой линейной алгебре алгоритм собственных значений Якоби представляет собой итеративный метод вычисления собственных значений и собственных векторов вещественной симметричной матрицы (процесс, известный как диагонализация ) . Он назван в честь Карла Густава Якоба Якоби , который впервые предложил этот метод в 1846 году, ^[1] но стал широко использоваться только в 1950-х годах с появлением компьютеров. ^[2]

Описание

Пусть – симметричная матрица, а – матрица вращения Гивенса . Затем: ${\displaystyle S}$ ${\displaystyle G=G(i,j,\theta )}$

${\displaystyle S'=GSG^{\top }\,}$

симметричен и подобен . ${\displaystyle S}$

Кроме того, есть записи: ${\displaystyle S^{\prime }}$

${\displaystyle {\begin{aligned}S'_{ii}&=c^{2}\,S_{ii}-2\,sc\,S_{ij}+s^{2}\,S_{jj}\\S'_{jj}&=s^{2}\,S_{ii}+2sc\,S_{ij}+c^{2}\,S_{jj}\\S'_{ij}&=S'_{ji}=(c^{2}-s^{2})\,S_{ij}+sc\,(S_{ii}-S_{jj})\\S'_{ik}&=S'_{ki}=c\,S_{ik}-s\,S_{jk}&k\neq i,j\\S'_{jk}&=S'_{kj}=s\,S_{ik}+c\,S_{jk}&k\neq i,j\\S'_{kl}&=S_{kl}&k,l\neq i,j\end{aligned}}}$

где и . ${\displaystyle s=\sin(\theta )}$ ${\displaystyle c=\cos(\theta )}$

Поскольку он ортогонален и имеет ту же норму Фробениуса (сумма квадратного корня квадратов всех компонентов), однако мы можем выбрать такой, что , и в этом случае имеет большую сумму квадратов на диагонали: ${\displaystyle G}$ ${\displaystyle S}$ ${\displaystyle S^{\prime }}$ ${\displaystyle ||\cdot ||_{F}}$ ${\displaystyle \theta }$ ${\displaystyle S_{ij}^{\prime }=0}$ ${\displaystyle S^{\prime }}$

${\displaystyle S'_{ij}=\cos(2\theta )S_{ij}+{\tfrac {1}{2}}\sin(2\theta )(S_{ii}-S_{jj})}$

Установите это значение равным 0 и переставьте:

${\displaystyle \tan(2\theta )={\frac {2S_{ij}}{S_{jj}-S_{ii}}}}$

если ${\displaystyle S_{jj}=S_{ii}}$

${\displaystyle \theta ={\frac {\pi }{4}}}$

Чтобы оптимизировать этот эффект, Sij должен быть _{недиагональным} элементом с наибольшим абсолютным значением, называемым опорной точкой .

Метод собственных значений Якоби неоднократно выполняет повороты , пока матрица не станет почти диагональной. Тогда элементы на диагонали являются аппроксимациями (действительных) собственных значений S .

Конвергенция

Если это опорный элемент, то по определению для . Пусть обозначает сумму квадратов всех недиагональных элементов . Поскольку имеет ровно недиагональные элементы, мы имеем или . Сейчас . Это подразумевает или ; то есть последовательность вращений Якоби сходится по крайней мере линейно с коэффициентом к диагональной матрице. ${\displaystyle p=S_{kl}}$ ${\displaystyle |S_{ij}|\leq |p|}$ ${\displaystyle 1\leq i,j\leq n,i\neq j}$ ${\displaystyle \Gamma (S)^{2}}$ ${\displaystyle S}$ ${\displaystyle S}$ ${\displaystyle 2N:=n(n-1)}$ ${\displaystyle p^{2}\leq \Gamma (S)^{2}\leq 2Np^{2}}$ ${\displaystyle 2p^{2}\geq \Gamma (S)^{2}/N}$ ${\displaystyle \Gamma (S^{J})^{2}=\Gamma (S)^{2}-2p^{2}}$ ${\displaystyle \Gamma (S^{J})^{2}\leq (1-1/N)\Gamma (S)^{2}}$ ${\displaystyle \Gamma (S^{J})\leq (1-1/N)^{1/2}\Gamma (S)}$ ${\displaystyle (1-1/N)^{1/2}}$

Ряд вращений Якоби называется разверткой; обозначим результат. Предыдущая оценка дает ${\displaystyle N}$ ${\displaystyle S^{\sigma }}$

${\displaystyle \Gamma (S^{\sigma })\leq \left(1-{\frac {1}{N}}\right)^{N/2}\Gamma (S)}$;

т. е. последовательность прогонов сходится по крайней мере линейно с множителем ≈ . ${\displaystyle e^{1/2}}$

Однако следующий результат Шенхаге ^[3] дает локально квадратичную сходимость. Для этого пусть S имеет m различных собственных значений с кратностями и пусть d > 0 — наименьшее расстояние между двумя разными собственными значениями. Давайте назовем несколько ${\displaystyle \lambda _{1},...,\lambda _{m}}$ ${\displaystyle \nu _{1},...,\nu _{m}}$

${\displaystyle N_{S}:={\frac {n(n-1)}{2}}-\sum _{\mu =1}^{m}{\frac {1}{2}}\nu _{\mu }(\nu _{\mu }-1)\leq N}$

Якоби вращает ход Шёнхаге. Если обозначает результат, то ${\displaystyle S^{s}}$

${\displaystyle \Gamma (S^{s})\leq {\sqrt {{\frac {n}{2}}-1}}\left({\frac {\gamma ^{2}}{d-2\gamma }}\right),\quad \gamma :=\Gamma (S)}$ .

Таким образом, сходимость становится квадратичной, как только ${\displaystyle \Gamma (S)<{\frac {d}{2+{\sqrt {{\frac {n}{2}}-1}}}}}$

Расходы

Каждое вращение Якоби можно выполнить за O( n ) шагов, если известен поворотный элемент p . Однако поиск p требует проверки всех N  ≈ 1/2 n ² недиагональных элемента. Мы также можем уменьшить это до сложности O( n ), если введем дополнительный массив индексов со свойством, которое является индексом самого большого элемента в строке i ( i = 1, ..., n  − 1) текущего S. . Тогда индексы пивота ( k , l ) должны быть одной из пар . Также обновление индексного массива может быть выполнено со средней сложностью O( n ): во-первых, максимальная запись в обновленных строках k и l может быть найдена за O( n ) шагов. В остальных строках i изменяются только записи в столбцах k и l . При циклическом просмотре этих строк, если нет ни k , ни l , достаточно сравнить старый максимум at с новыми записями и при необходимости обновить. Если должно быть равно k или l и соответствующая запись уменьшилась во время обновления, максимум по строке i необходимо найти с нуля со сложностью O( n ). Однако это произойдет в среднем только один раз за оборот. Таким образом, каждое вращение имеет сложность среднего случая O( n ) и одну прогонку O( n ³ ), что эквивалентно одному умножению матрицы. Кроме того , перед запуском процесса необходимо инициализировать его , что можно сделать за n ² шагов. ${\displaystyle m_{1},\,\dots \,,\,m_{n-1}}$ ${\displaystyle m_{i}}$ ${\displaystyle (i,m_{i})}$ ${\displaystyle m_{i}}$ ${\displaystyle m_{i}}$ ${\displaystyle m_{i}}$ ${\displaystyle m_{i}}$ ${\displaystyle m_{i}}$

Обычно метод Якоби сходится с точностью до числа после небольшого количества проходов. Обратите внимание, что несколько собственных значений уменьшают количество итераций, поскольку . ${\displaystyle N_{S}<N}$

Алгоритм

Следующий алгоритм представляет собой описание метода Якоби в математических обозначениях. Он вычисляет вектор e , содержащий собственные значения, и матрицу E , содержащую соответствующие собственные векторы; то есть является собственным значением, а столбец — ортонормированным собственным вектором для , i = 1, ..., n . ${\displaystyle e_{i}}$ ${\displaystyle E_{i}}$ ${\displaystyle e_{i}}$

процедура jacobi( S ∈ R n × n ; out  e ∈ R n ; out  E ∈ R n × n ) вар  я , k , л , м , состояние ∈ N  s , c , т , п , y , d , р ∈ R  ind ∈ N n  изменено ∈ L n функция maxind( k ∈ N ) ∈ N  ! индекс наибольшего недиагонального элемента в строке k  m  := k +1 for  i  := k +2 to  n  do  if │ S ki │ > │ S km │ then  m  := i  endif  endfor  return  m  endfunc обновление процедуры ( k ∈ N ; т ∈ R ) ! обновить e k и его статус  y  := e k ; e k  := y + t,  если  изменилось k и ( y = e k ) , то  изменилось k  := false; состояние  := состояние −1 elsif (не изменено k ) и ( y ≠ e k ) затем  изменено k  := true; состояние  := состояние +1 endif  endproc процедура Rotate( k , l , i , j ∈ N ) ! выполнить вращение Sij , S kl ┌  ┐ ┌ ┐┌  ┐ │ S кл │ │ c − s ││ S кл │ │  │ := │ ││  │ │ S ij │ │ s  c ││ S ij │ └  ┘ └ ┘└  ┘ endproc  ! init e, E и массивы ind, изменены  E  := I ; состояние  := n  for  k  := от 1 до  n  do  ind k  := maxind( k ); е к  := S kk ; изменено k  := true endfor  while  state ≠0 do  ! следующий оборот  m  := 1 ! найти индекс (k,l) опорной точки p  для  k  := 2 до  n −1 do  if │ S k  ind k │ > │ S m  ind m │ then  m  := k  endif  endfor  k  := m ; л  := инд м ; п  := S кл  ! вычислить c = cos φ, s = sin φ  y  := ( e l − e k )/2; d  := │ y │+√( п 2 + y 2 ) р  := √( п 2 + d 2 ); с  := д / р ; с  := п / р ; t  := p 2 / d,  если  y <0 , то  s  := − s ; t  := − t  endif  S kl  := 0,0; обновление( к ,- т ); обновление( л , т )  ! повернуть строки и столбцы k и l  для  i  := 1 до  k −1 do Rotate( i , k , i , l ) endfor  for  i  := k +1 до  l −1 do Rotate( k , i , i , l ) endfor  for  i  := l +1 to  n  do Rotate( k , i , l , i ) endfor  ! повернуть собственные векторы  для  i  := 1 в  n  do ┌  ┐ ┌ ┐┌  ┐ │ E ik │ │ c − s ││ E ik │ │  │ := │ ││  │ │ Э иль │ │ s  c ││ Э иль │ └  ┘ └ ┘└  ┘ конец для  ! обновить все потенциально измененные ind i  for  i  := 1 до  n  do  ind i  := maxind( i ) endfor  цикла endproc

Примечания

1. Измененный логический массив сохраняет статус каждого собственного значения. Если числовое значение или изменяется во время итерации, соответствующему компоненту изменено присваивается значение true , в противном случае — значение false . Целочисленное состояние подсчитывает количество компонентов изменения , имеющих значение true . Итерация прекращается, как только состояние = 0. Это означает, что ни одно из приближений в последнее время не изменило свое значение, и поэтому маловероятно, что это произойдет, если итерация продолжится. Здесь предполагается, что операции с плавающей запятой оптимально округляются до ближайшего числа с плавающей запятой. ${\displaystyle e_{k}}$ ${\displaystyle e_{l}}$ ${\displaystyle e_{1},\,...\,,e_{n}}$

2. Верхний треугольник матрицы S разрушается, а нижний треугольник и диагональ остаются неизменными. Таким образом, при необходимости можно восстановить S согласно

для  k  := 1 до  n −1 do  ! восстановить матрицу 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  если  k ≠ m  то поменять местами e m , e k поменять местами E m , E k  endif endfor

4. Алгоритм написан с использованием матричной записи (массивы с отсчетом от 1 вместо 0).

5. При реализации алгоритма часть, заданная с помощью матричной записи, должна выполняться одновременно.

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

Пример

Позволять ${\displaystyle S={\begin{pmatrix}4&-30&60&-35\\-30&300&-675&420\\60&-675&1620&-1050\\-35&420&-1050&700\end{pmatrix}}}$

Затем jacobi выдает следующие собственные значения и собственные векторы после 3 проходов (19 итераций):

${\displaystyle e_{1}=2585.25381092892231}$

${\displaystyle E_{1}={\begin{pmatrix}0.0291933231647860588\\-0.328712055763188997\\0.791411145833126331\\-0.514552749997152907\end{pmatrix}}}$

${\displaystyle e_{2}=37.1014913651276582}$

${\displaystyle E_{2}={\begin{pmatrix}-0.179186290535454826\\0.741917790628453435\\-0.100228136947192199\\-0.638282528193614892\end{pmatrix}}}$

${\displaystyle e_{3}=1.4780548447781369}$

${\displaystyle E_{3}={\begin{pmatrix}-0.582075699497237650\\0.370502185067093058\\0.509578634501799626\\0.514048272222164294\end{pmatrix}}}$

${\displaystyle e_{4}=0.1666428611718905}$

${\displaystyle E_{4}={\begin{pmatrix}0.792608291163763585\\0.451923120901599794\\0.322416398581824992\\0.252161169688241933\end{pmatrix}}}$

Приложения для вещественных симметричных матриц

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

Сингулярные значения

Сингулярные значения (квадратной) матрицы являются квадратными корнями (неотрицательных) собственных значений матрицы . В случае симметричной матрицы мы имеем , следовательно, сингулярные значения являются абсолютными значениями собственных значений ${\displaystyle A}$ ${\displaystyle A^{T}A}$ ${\displaystyle S}$ ${\displaystyle S^{T}S=S^{2}}$ ${\displaystyle S}$ ${\displaystyle S}$

2-норма и спектральный радиус

2-норма матрицы A — это норма, основанная на евклидовой векторной норме; то есть наибольшее значение, когда x проходит через все векторы с . Это самое большое единственное значение . В случае симметричной матрицы это наибольшее абсолютное значение ее собственных векторов и, следовательно, равное ее спектральному радиусу. ${\displaystyle \|Ax\|_{2}}$ ${\displaystyle \|x\|_{2}=1}$ ${\displaystyle A}$

Номер условия

Число обусловленности неособой матрицы определяется как . В случае симметричной матрицы это абсолютное значение частного наибольшего и наименьшего собственного значения. Матрицы с большими числами обусловленности могут привести к численно нестабильным результатам: небольшое возмущение может привести к большим ошибкам. Матрицы Гильберта — самые известные плохо обусловленные матрицы. Например, матрица Гильберта четвертого порядка имеет условие 15514, а для 8-го порядка — 2,7×10 ⁸ . ${\displaystyle A}$ ${\displaystyle {\mbox{cond}}(A)=\|A\|_{2}\|A^{-1}\|_{2}}$

Классифицировать

Матрица имеет ранг , если в ней есть столбцы, которые линейно независимы, а остальные столбцы линейно зависят от них. Эквивалентно, это размер диапазона  . Кроме того, это количество ненулевых сингулярных значений. ${\displaystyle A}$ ${\displaystyle r}$ ${\displaystyle r}$ ${\displaystyle r}$ ${\displaystyle A}$

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

Псевдообратный

Псевдообратная матрица — это единственная матрица, для которой и симметричны и для которой выполняется. Если неособо, то . ${\displaystyle A}$ ${\displaystyle X=A^{+}}$ ${\displaystyle AX}$ ${\displaystyle XA}$ ${\displaystyle AXA=A,XAX=X}$ ${\displaystyle A}$ ${\displaystyle A^{+}=A^{-1}}$

При вызове процедуры jacobi (S, e, E) выполняется соотношение, где Diag( e ) обозначает диагональную матрицу с вектором e на диагонали. Обозначим вектор, где заменяется на if и на 0, если (численно близко) к нулю. Поскольку матрица E ортогональна, отсюда следует, что псевдообратная матрица S определяется выражением . ${\displaystyle S=E^{T}{\mbox{Diag}}(e)E}$ ${\displaystyle e^{+}}$ ${\displaystyle e_{i}}$ ${\displaystyle 1/e_{i}}$ ${\displaystyle e_{i}\leq 0}$ ${\displaystyle e_{i}}$ ${\displaystyle S^{+}=E^{T}{\mbox{Diag}}(e^{+})E}$

Решение методом наименьших квадратов

Если матрица не имеет полного ранга, решения линейной системы может не быть . Однако можно найти вектор x, для которого он минимален. Решение такое . В случае симметричной матрицы S, как и раньше, имеем . ${\displaystyle A}$ ${\displaystyle Ax=b}$ ${\displaystyle \|Ax-b\|_{2}}$ ${\displaystyle x=A^{+}b}$ ${\displaystyle x=S^{+}b=E^{T}{\mbox{Diag}}(e^{+})Eb}$

Матричная экспонента

Отсюда можно найти , где exp  — вектор, где заменяется на . Точно так же очевидным образом можно вычислить любую (аналитическую) функцию . ${\displaystyle S=E^{T}{\mbox{Diag}}(e)E}$ ${\displaystyle \exp S=E^{T}{\mbox{Diag}}(\exp e)E}$ ${\displaystyle e}$ ${\displaystyle e_{i}}$ ${\displaystyle \exp e_{i}}$ ${\displaystyle f(S)}$ ${\displaystyle f}$

Линейные дифференциальные уравнения

Дифференциальное уравнение имеет решение . Для симметричной матрицы отсюда следует, что . Если – разложение по собственным векторам , то . ${\displaystyle x'=Ax,x(0)=a}$ ${\displaystyle x(t)=\exp(tA)}$ ${\displaystyle S}$ ${\displaystyle x(t)=E^{T}{\mbox{Diag}}(\exp te)Ea}$ ${\displaystyle a=\sum _{i=1}^{n}a_{i}E_{i}}$ ${\displaystyle a}$ ${\displaystyle S}$ ${\displaystyle x(t)=\sum _{i=1}^{n}a_{i}\exp(te_{i})E_{i}}$

Пусть - векторное пространство, натянутое на собственные векторы, которые соответствуют отрицательному собственному значению и аналогично положительным собственным значениям. Если тогда ; то есть точка равновесия 0 притягивает к . Если тогда ; то есть 0 отталкивает . и называются устойчивыми и неустойчивыми многообразиями для . Если компоненты есть в обоих многообразиях, то один компонент притягивается, а другой отталкивается. Следовательно подходит как . ${\displaystyle W^{s}}$ ${\displaystyle S}$ ${\displaystyle W^{u}}$ ${\displaystyle a\in W^{s}}$ ${\displaystyle {\mbox{lim}}_{t\rightarrow \infty }x(t)=0}$ ${\displaystyle x(t)}$ ${\displaystyle a\in W^{u}}$ ${\displaystyle {\mbox{lim}}_{t\rightarrow \infty }x(t)=\infty }$ ${\displaystyle x(t)}$ ${\displaystyle W^{s}}$ ${\displaystyle W^{u}}$ ${\displaystyle S}$ ${\displaystyle a}$ ${\displaystyle x(t)}$ ${\displaystyle W^{u}}$ ${\displaystyle t\to \infty }$

Реализация Юлии

Следующий код представляет собой прямую реализацию математического описания алгоритма собственных значений Якоби на языке программирования Julia .

используя LinearAlgebra , Тест  функция find_pivot ( Sprime ) n = размер ( Sprime , 1 ) Pivot_i = Pivot_j = 0 Pivot = 0,0             for j = 1 : n for i = 1 : ( j - 1 ) if abs ( Sprime [ i , j ]) > Pivot Pivot_i = i Pivot_j = J Pivot = ABS ( Sprime [ i , J ]) end end end                        return ( Pivot_i , Pivot_j , Pivot ) конец   # на практике не следует явно создавать экземпляр матричной функции вращения Гивенса .​​​​​​​​​​​​​​​​​​​​​​​ ] = потому что ( θ ) грамм [ я , j ] = грех ( θ ) грамм [ j , я ] = - грех ( θ ) вернуть G конец                 # S — симметричная матрица размера n на n n = 4 sqrtS = randn ( n , n ); S = sqrtS * sqrtS ' ;      # наибольший разрешенный недиагональный элемент U' * S * U # где U — собственные векторы tol = 1e-14  Sprime = копия ( S ) U = Матрица { Float64 } ( I , ( n , n ))    в то время как true ( Pivot_i , Pivot_j , Pivot ) = find_pivot ( Sprime )       если точка поворота < tol конец разрыва      θ = atan ( 2 * Спрайм [ Pivot_i , Pivot_j ] / ( Спрайм [ Pivot_j , Pivot_j ] - Спрайм [ Pivot_i , Pivot_i ] )) / 2        G = заданная_матрица_вращения ( n , Pivot_i , Pivot_j , θ )   # обновляем Sprime и U Sprime .= G '* Sprime * G U .= U * G end        # Sprime теперь (почти) диагональная матрица # извлекаем собственные значения λ = diag ( Sprime )  # сортируем собственные значения (и соответствующие собственные векторы U) по возрастанию значений i = sortperm ( λ ) λ = λ [ i ] U = U [ : , i ]      # S должно быть равно U *diagm(λ) * U' @test S ≈ U * diagm ( λ ) * U '       

Обобщения

Метод Якоби был обобщен на комплексные эрмитовы матрицы , общие несимметричные вещественные и комплексные матрицы, а также на блочные матрицы.

Поскольку сингулярные значения вещественной матрицы являются квадратными корнями собственных значений симметричной матрицы, их также можно использовать для расчета этих значений. Для этого случая метод модифицируется таким образом, что S не нужно вычислять явно, что снижает опасность ошибок округления . Обратите внимание, что с . ${\displaystyle S=A^{T}A}$ ${\displaystyle JSJ^{T}=JA^{T}AJ^{T}=JA^{T}J^{T}JAJ^{T}=B^{T}B}$ ${\displaystyle B\,:=JAJ^{T}}$

Метод Якоби также хорошо подходит для параллелизма.

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

1. Якоби, CGJ (1846). «Über ein leichtes Verfahren, die in der Theorie der Säkularstörungen vorkommenden Gleichungen numerisch aufzulösen». Журнал Крелля (на немецком языке). 1846 (30): 51–94. дои : 10.1515/crll.1846.30.51. S2CID  199546177.
2. Голуб, Г.Х.; ван дер Ворст, HA (2000). «Вычисление собственных значений в 20 веке». Журнал вычислительной и прикладной математики . 123 (1–2): 35–65. дои : 10.1016/S0377-0427(00)00413-1 .
3. Шёнхаге, А. (1964). «Квадратная конвергенция Якоби-Верфаренса». Numerische Mathematik (на немецком языке). 6 (1): 410–412. дои : 10.1007/BF01386091. MR  0174171. S2CID  118301078.

дальнейшее чтение

• Пресс, WH; Теукольский, С.А.; Феттерлинг, WT; Фланнери, Б.П. (2007), «Раздел 11.1. Преобразования Якоби симметричной матрицы», Численные рецепты: искусство научных вычислений (3-е изд.), Нью-Йорк: издательство Кембриджского университета, ISBN 978-0-521-88068-8
• Рутисхаузер, Х. (1966). «Серия справочников по линейной алгебре: метод Якоби для вещественных симметричных матриц». Нумерическая математика . 9 (1): 1–10. дои : 10.1007/BF02165223. MR  1553948. S2CID  120520713.
• Самех, А.Х. (1971). «О Якоби и якобиподобных алгоритмах для параллельного компьютера». Математика вычислений . 25 (115): 579–590. дои : 10.1090/s0025-5718-1971-0297131-6 . JSTOR  2005221. МР  0297131.
• Шрофф, Гаутам М. (1991). «Параллельный алгоритм для собственных значений и собственных векторов общей комплексной матрицы». Числовая математика . 58 (1): 779–805. CiteSeerX  10.1.1.134.3566 . дои : 10.1007/BF01385654. MR  1098865. S2CID  13904356.
• Веселич, К. (1979). «Об одном классе якобиподобных процедур диагонализации произвольных вещественных матриц». Нумерическая математика . 33 (2): 157–172. дои : 10.1007/BF01399551. MR  0549446. S2CID  119919630.
• Веселич, К.; Венцель, HJ (1979). «Квадратично сходящийся метод Якоби для вещественных матриц с комплексными собственными значениями». Нумерическая математика . 33 (4): 425–435. дои : 10.1007/BF01399324. MR  0553351. S2CID  119554420.
• Юсеф Саад: «Возвращаясь к (блочному) методу вращения подпространства Якоби для симметричной проблемы собственных значений», Numerical Algorithms, vol.92 (2023), стр.917-944. https://doi.org/10.1007/s11075-022-01377-w.

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

• Реализация в Matlab алгоритма Якоби, позволяющая избежать тригонометрических функций
• реализация С++11