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

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

Описание

Позволять $S$ быть симметричной матрицей, и $G=G(i,j,\theta )$ быть матрицей вращения Гивенса . Затем:

$S'=G^{\top }SG\,$

симметричен подобен и $S$.

Более того, $S^{\prime }$ есть записи:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Расходы

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

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

Алгоритм

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

процедура  jacobi(  S  ∈  R n × n ;  out   e  ∈  R н ;  out   E  ∈  R n × n )   вар      я  ,  k  ,  л  ,  м  ,  состояние  ∈  N      s  ,  c  ,  т  ,  п  ,  y  ,  d  ,  р  ∈  R      ind  ∈  N н     изменилось  ∈  L н   функция  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    процедура  update(  k  ∈  N  ;  t  ∈  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  ) !  выполнить вращение S  ij  , S  kl     ┌    ┐ ┌ ┐┌    ┐    │  S  кл  │ │  c  −  s  ││  S  кл  │    │    │ := │ ││    │    │  S  ij  │ │  s     c  ││  S  ij  │    └    ┘ └ ┘└    ┘   конечный процесс   !  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  │+√(  p 2 +  и 2 )     р  := √(  р 2 +  д 2 );  с  :=  д  /  р  ;  с  :=  п  /  р  ;  т  :=  р 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. Измененный сохраняет статус каждого собственного значения. Если числовое значение $e_{k}$ или $e_{l}$ изменяется во время итерации, соответствующий компонент изменения устанавливается в значение true , в противном случае — в значение false . Целочисленное состояние подсчитывает количество компонентов изменения , имеющих значение true . Итерация прекращается, как только состояние = 0. Это означает, что ни одно из приближений $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  до   n  −1  do      m  :=  k      for   l  :=  k  +1  до   n   do          if   e  l  >  e  m   then              m  :=  l          endif      endfor      если   k  ≠  m   тогда         поменять  e  m  местами  и  k         поменять  E  m  ,  E  k      endif  endfor 

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

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

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

Пример

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

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

$e_{1}=2585.25381092892231$

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

$e_{2}=37.1014913651276582$

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

$e_{3}=1.4780548447781369$

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

$e_{4}=0.1666428611718905$

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

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

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

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

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

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

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

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

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

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

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

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

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

Псевдообратная матрица $A$ это уникальная матрица $X=A^{+}$ для чего $AX$ и $XA$ симметричны и для которых $AXA=A,XAX=X$ держит. Если $A$ неособа, то $A^{+}=A^{-1}$.

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

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

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

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

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

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

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

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

Реализация Джулии

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

используя   LinearAlgebra  ,   тестовая  функция   find_pivot  (  Sprime  )      n   =   size  (  Sprime  ,  1  )      Pivot_i   =   Pivot_j   =   0      Pivot   =   0,0      для   j   =   1  :  n          для   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  )  end  # на практике не следует явно создавать экземпляр матричной  функции   вращения Гивенса Gives_rotation_matrix  (  n  ,  i  ,  j  ,  θ  )      G   знак равно   Матрица  {  Float64  }(  я  ,(  n  ,  n  ))      грамм  [  я  ,  я  ]   знак равно   грамм  [  j  ,  j  ]   знак   грех  равно потому что (  θ  )      грамм  [  я  ,  j  ]   знак   равно  (  θ  )      грамм  [  j  ,  i  ]   =   -  sin  (  θ  )      return   G  end  # S — симметричная матрица размера n на n  n   =   4  sqrtS   =   randn  (  n  ,  n  );  S   =   sqrtS  *  sqrtS  '  ;  # наибольший разрешенный недиагональный элемент U' * S * U  # где U — собственные векторы  tol   =   1e-14  Sprime   =   copy  (  S  )  U   =   Matrix  {  Float64  }(  I  ,(  n  ,  n  ))  while   true      (  Pivot_i  ,   Pivot_j  ,   Pivot  )   =   find_pivot  (  Sprime  ),      если   точка вращения   <   tol,          разрыва      конец      θ   =   atan  (  2  *  Sprime  [  Pivot_i  ,  Pivot_j  ]  /  (  Sprime  [  Pivot_j  ,  Pivot_j  ]   -   Sprime  [  Pivot_i  ,  Pivot_i  ]   ))   /   2      G   =   Gives_rotation_matrix  (  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=A^{T}A$ его также можно использовать для расчета этих значений. Для этого случая метод модифицируется таким образом, что S не нужно вычислять явно, что снижает опасность ошибок округления . Обратите внимание, что $JSJ^{T}=JA^{T}AJ^{T}=JA^{T}J^{T}JAJ^{T}=B^{T}B$ с $B\,:=JAJ^{T}$ .

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

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

• Пресс, WH; Теукольский, С.А.; Феттерлинг, WT; Фланнери, Б.П. (2007), «Раздел 11.1. Преобразования Якоби симметричной матрицы» , Численные рецепты: искусство научных вычислений (3-е изд.), Нью-Йорк: Издательство Кембриджского университета, ISBN  978-0-521-88068-8
• Рутисхаузер, Х. (1966). «Серия справочников по линейной алгебре: метод Якоби для вещественных симметричных матриц». Нумерическая математика . 9 (1): 1–10. дои : 10.1007/BF02165223 . МР   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 . МР   1098865 . S2CID   13904356 .
• Веселич, К. (1979). «Об одном классе якобиподобных процедур диагонализации произвольных вещественных матриц». Нумерическая математика . 33 (2): 157–172. дои : 10.1007/BF01399551 . МР   0549446 . S2CID   119919630 .
• Веселич, К.; Венцель, HJ (1979). «Квадратично сходящийся метод Якоби для вещественных матриц с комплексными собственными значениями». Нумерическая математика . 33 (4): 425–435. дои : 10.1007/BF01399324 . МР   0553351 . S2CID   119554420 .
• Юсеф Саад: «Возвращаясь к (блочному) методу вращения подпространства Якоби для симметричной проблемы собственных значений», Numerical Algorithms, vol.92 (2023), стр.917-944. https://doi.org/10.1007/s11075-022-01377-w .

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

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