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

В числовой линейной алгебре алгоритм собственных значений Якоби представляет собой итерационный метод вычисления собственных значений и собственных векторов вещественной диагонализация симметричной матрицы (процесс, известный как ) . Он назван в честь Карла Густава Якоба Якоби , впервые предложившего метод в 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 2 н ² недиагональные элементы. Мы также можем уменьшить эту сложность до 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, ..., п .

procedure jacobi(S ∈ Rn×n; out e ∈ Rn; out E ∈ Rn×n)
  var
    i, k, l, m, state ∈ N
    s, c, t, p, y, d, r ∈ R
    ind ∈ Nn
    changed ∈ Ln

  function maxind(k ∈ N) ∈ N ! index of largest off-diagonal element in row k
    m := k+1
    for i := k+2 to n do
      if │Ski│ > │Skm│ then m := i endif
    endfor
    return m
  endfunc

  procedure update(k ∈ N; t ∈ R) ! update ek and its status
    y := ek; ek := y+t
    if changedk and (y=ek) then changedk := false; state := state−1
    elsif (not changedk) and (y≠ek) then changedk := true; state := state+1
    endif
  endproc

  procedure rotate(k,l,i,j ∈ N) ! perform rotation of Sij, Skl
    ┌   ┐    ┌     ┐┌   ┐
    │Skl│    │c  −s││Skl│
    │   │ := │     ││   │
    │Sij│    │s   c││Sij│
    └   ┘    └     ┘└   ┘
  endproc

  ! init e, E, and arrays ind, changed
  E := I; state := n
  for k := 1 to n do indk := maxind(k); ek := Skk; changedk := true endfor
  while state≠0 do ! next rotation
    m := 1 ! find index (k,l) of pivot p
    for k := 2 to n−1 do
      if │Sk indk│ > │Sm indm│ then m := k endif
    endfor
    k := m; l := indm; p := Skl
    ! calculate c = cos φ, s = sin φ
    y := (el−ek)/2; d := │y│+√(p2+y2)
    r := √(p2+d2); c := d/r; s := p/r; t := p2/d
    if y<0 then s := −s; t := −t endif
    Skl := 0.0; update(k,−t); update(l,t)
    ! rotate rows and columns k and l
    for i := 1 to k−1 do rotate(i,k,i,l) endfor
    for i := k+1 to l−1 do rotate(k,i,i,l) endfor
    for i := l+1 to n do rotate(k,i,l,i) endfor
    ! rotate eigenvectors
    for i := 1 to n do
      ┌   ┐    ┌     ┐┌   ┐
      │Eik│    │c  −s││Eik│
      │   │ := │     ││   │
      │Eil│    │s   c││Eil│
      └   ┘    └     ┘└   ┘
    endfor
    ! update all potentially changed indi
    for i := 1 to n do indi := maxind(i) endfor
  loop
endproc

Примечания

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

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

for k := 1 to n−1 do ! restore matrix S
    for l := k+1 to n do
        Skl := Slk
    endfor
endfor

3. Собственные значения не обязательно расположены в порядке убывания. Этого можно добиться с помощью простого алгоритма сортировки.

for k := 1 to n−1 do
    m := k
    for l := k+1 to n do
        if el > em then
            m := l
        endif
    endfor
    if k ≠ m then
        swap em,ek
        swap Em,Ek
    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 .

using LinearAlgebra, Test

function find_pivot(Sprime)
    n = size(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)
end

# in practice one should not instantiate explicitly the Givens rotation matrix
function givens_rotation_matrix(n,i,j,θ)
    G = Matrix{Float64}(I,(n,n))
    G[i,i] = G[j,j] = cos(θ)
    G[i,j] = sin(θ)
    G[j,i] = -sin(θ)
    return G
end

# S is a symmetric n by n matrix
n = 4
sqrtS = randn(n,n);
S = sqrtS*sqrtS';

# the largest allowed off-diagonal element of U' * S * U
# where U are the eigenvectors
tol = 1e-14

Sprime = copy(S)
U = Matrix{Float64}(I,(n,n))

while true
    (pivot_i, pivot_j, pivot) = find_pivot(Sprime)

    if pivot < tol
        break
    end

    θ = atan(2*Sprime[pivot_i,pivot_j]/(Sprime[pivot_j,pivot_j] - Sprime[pivot_i,pivot_i] )) / 2

    G = givens_rotation_matrix(n,pivot_i,pivot_j,θ)

    # update Sprime and U
    Sprime .= G'*Sprime*G
    U .= U * G
end

# Sprime is now (almost) a diagonal matrix
# extract eigenvalues
λ = diag(Sprime)

# sort eigenvalues (and corresponding eigenvectors U) by increasing values
i = sortperm(λ)
λ = λ[i]
U = U[:,i]

# S should be equal to 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