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

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

Описание

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

{\ Displaystyle S '= GSG ^ {\ top} \,}

симметричен и похож на . ${\ displaystyle S}$

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

{\ displaystyle {\ begin {align} 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 {выровнено}}}

где и . ${\ Displaystyle s = \ грех (\ тета)}$ ${\ Displaystyle с = \ соз (\ тета)}$

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

{\ Displaystyle S '_ {ij} = \ соз (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}}}

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

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

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

Если - это элемент поворота, то по определению для . Позвольте обозначить сумму квадратов всех недиагональных элементов . Поскольку имеет ровно недиагональные элементы, мы имеем или . Сейчас . Отсюда следует или , т. Е. Последовательность поворотов Якоби сходится по крайней мере линейно с множителем к диагональной матрице. ${\ displaystyle p = S_ {kl}}$ ${\ Displaystyle | S_ {ij} | \ leq | p |}$ ${\ Displaystyle 1 \ Leq я, j \ leq п, я \ neq j}$ ${\ displaystyle \ Gamma (S) ^ {2}}$ ${\ displaystyle S}$ ${\ displaystyle S}$ ${\ Displaystyle 2N: = п (п-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 \ Гамма (S ^ {J}) \ Leq (1-1 / N) ^ {1/2} \ Гамма (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}}$

Однако следующий результат Шёнхаге дает локально квадратичную сходимость. Для этого пусть 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 м_ {1}, \, \ точки \ ,, \, м_ {п-1}}$ ${\ displaystyle m_ {i}}$ ${\ Displaystyle (я, м_ {я})}$ ${\ 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}}$

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

  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 │S_ki│ > │S_km│ then m := i endif
    endfor
    return m
  endfunc

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

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

  ! init e, E, and arrays ind, changed
  E := I; state := n
  for k := 1 to n do ind_k := maxind(k); e_k := S_kk; changed_k := 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 │S_{k ind_k}│ > │S_{m ind_m}│ then m := k endif
    endfor
    k := m; l := ind_m; p := S_kl
    ! calculate c = cos φ, s = sin φ
    y := (e_l−e_k)/2; d := │y│+√(p²+y²)
    r := √(p²+d²); c := d/r; s := p/r; t := p²/d
    if y<0 then s := −s; t := −t endif
    S_kl := 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
      ┌ ┐    ┌     ┐┌ ┐
      │E_ik│    │c  −s││E_ik│
      │ │ := │     ││ │
      │E_il│    │s   c││E_il│
      └ ┘    └     ┘└ ┘
    endfor
    ! rows k, l have changed, update rows ind_k, ind_l
    ind_k := maxind(k); ind_l := maxind(l)
  loop
endproc

Заметки

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

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

for k := 1 to n−1 do ! restore matrix S
    for l := k+1 to 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. Эта реализация неправильно учитывает случай, когда одно измерение является независимым подпространством. Например, если задана диагональная матрица, вышеуказанная реализация никогда не завершится, поскольку ни одно из собственных значений не изменится. Следовательно, в реальных реализациях необходимо добавить дополнительную логику для учета этого случая.

Пример

Позволять ${\ 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}}}$

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

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

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

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

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

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

Псевдообратный: Псевдообратная матрица A - это единственная матрица, для которой AX и XA симметричны и для которой выполняется AXA = A, XAX = X. Если A неособо, то ' . ${\ displaystyle X = A ^ {+}}$ ${\ Displaystyle А ^ {+} = А ^ {- 1}}$; Когда вызывается процедура jacobi (S, e, E), то выполняется соотношение, где Diag ( e ) обозначает диагональную матрицу с вектором e на диагонали. Пусть обозначаться вектор , где заменяется , если и 0 , если есть (численно близка к нулю). Поскольку матрица E ортогональна, псевдообратная матрица S имеет вид . ${\ Displaystyle S = E ^ {T} {\ mbox {Diag}} (е) E}$ ${\ displaystyle e ^ {+}}$ ${\ displaystyle e_ {i}}$ ${\ displaystyle 1 / e_ {i}}$ ${\ Displaystyle е_ {я} \ leq 0}$ ${\ displaystyle e_ {i}}$ ${\ Displaystyle S ^ {+} = E ^ {T} {\ mbox {Diag}} (e ^ {+}) E}$

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

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

Линейные дифференциальные уравнения: Дифференциальное уравнение x ' = Ax , x (0) = a имеет решение x ( t ) = exp ( t A ) a . Для симметричной матрицы S отсюда следует, что . Если есть разложение a по собственным векторам S , то . ${\ Displaystyle х (т) = E ^ {T} {\ mbox {Diag}} (\ exp te) Ea}$ ${\ Displaystyle а = \ сумма _ {я = 1} ^ {п} а_ {я} Е_ {я}}$ ${\ displaystyle x (t) = \ sum _ {i = 1} ^ {n} a_ {i} \ exp (te_ {i}) E_ {i}}$; Позвольте быть векторным пространством, натянутым на собственные векторы S, которые соответствуют отрицательному собственному значению и аналогично положительным собственным значениям. Если тогда, т.е. точка равновесия 0 привлекательна для x ( t ). Если тогда , то есть 0 отталкивает x ( t ). и называются устойчивые и неустойчивые многообразия для S . Если a имеет компоненты в обоих многообразиях, то один компонент притягивается, а один компонент отталкивается. Следовательно, x ( t ) приближается как . ${\ Displaystyle W ^ {s}}$ ${\ displaystyle W ^ {u}}$ ${\ displaystyle a \ in W ^ {s}}$ ${\ Displaystyle {\ mbox {lim}} _ {т \ \ infty} х (т) = 0}$ ${\ displaystyle a \ in W ^ {u}}$ ${\ Displaystyle {\ mbox {lim}} _ {т \ \ infty} х (т) = \ infty}$ ${\ Displaystyle W ^ {s}}$ ${\ displaystyle W ^ {u}}$ ${\ displaystyle W ^ {u}}$ ${\ Displaystyle т \ к \ infty}$

Обобщения

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

Поскольку сингулярные значения реальной матрицы являются квадратными корнями из собственных значений симметричной матрицы, их также можно использовать для вычисления этих значений. В этом случае метод модифицируется таким образом, что 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}}$

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

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

Нажмите, 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 .

Languages

In other projects