Algorytm wartości własnych Jacobiego

W numerycznej algebrze liniowej algorytm wartości własnej Jacobiego jest iteracyjną metodą obliczania wartości własnych i wektorów własnych rzeczywistej macierzy symetrycznej (proces znany jako diagonalizacja ). Jej nazwa pochodzi od Carla Gustava Jacoba Jacobiego , który jako pierwszy zaproponował tę metodę w 1846 roku, ale stała się szeroko stosowana dopiero w latach pięćdziesiątych XX wieku wraz z pojawieniem się komputerów.

Opis

Niech $_$$macierzą$ i _ _ _ Następnie:

${\ Displaystyle S '= GSG ^ {\ wierzchołek} \,}$

jest symetryczny i podobny do ${\ displaystyle S}$ .

Ponadto ${\ displaystyle S ^ {\ pierwsza}}$ ma wpisy: S ′ {\ Displaystyle S ^ {\ pierwsza}}

${\ Displaystyle {\ rozpocząć {wyrównane} 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{wyrównane}}}$

gdzie $)}$${\ Displaystyle c = \ cos (\ teta)}$ i do .

Ponieważ $ortogonalny$ , i $mają$$Frobeniusa$ samą _ ${\ Displaystyle ||\ cdot ||_ {F}}$ (pierwiastkowa suma kwadratów wszystkich składników), jednak możemy wybrać ${\ displaystyle S_ {ij}^{\prime}=0}$$jot$ że , w takim przypadku ${\ Displaystyle S ^ {\ liczba pierwsza}}$ ma większą sumę kwadratów na przekątnej:

${\ Displaystyle S'_ {ij} = \ sałata (2 \ teta) S_ {ij}+{\tfrac {1}{2}}\sin(2\theta )(S_{ii}-S_{jj})}$

Ustaw to na 0 i zmień kolejność:

${\ Displaystyle \ dębnik (2 \ teta) = {\ Frac {2S_ {ij}} {S_ {jj} -S_ {ii}}} }$

jeśli ${\ Displaystyle S_ {jj} = S_ {ii}}$

${\ Displaystyle \ teta = {\ Frac {\ pi} {4}}}$

Aby zoptymalizować ten efekt, S _ij powinien być elementem pozadiagonalnym o największej wartości bezwzględnej, zwanym sworzniem .

Metoda wartości własnej Jacobiego wielokrotnie wykonuje obroty , aż macierz stanie się prawie diagonalna. Wtedy elementy na przekątnej są przybliżeniami (rzeczywistych) wartości własnych S .

Konwergencja

Jeśli ${\ displaystyle p = S_ {kl}}$ jest elementem obrotowym, to z definicji ${\ Displaystyle | S_ {ij} | \ równoważnik | p |}$ dla ${\ Displaystyle 1 \ równoważnik ja, j \ równoważnik n, i \ neq j}$ . Niech ${\ Displaystyle \ Gamma (S) ^ {2}}$ oznaczają sumę kwadratów wszystkich wpisów $przekątną$ . Ponieważ ${\ Displaystyle p ^ {2} \ równoważnik \ Gamma (S) ^ {2} \ równoważnik 2Np ^ {2}}$ dokładnie $Γ$ poza $przekątną$ lub ${\ Displaystyle 2p ^ {2} \ geq \ Gamma (S) ^ {2} / N}$ . teraz ${\ Displaystyle \ Gamma (S ^ {J}) ^ {2} = \ Gamma (S) ^ {2} -2 p ^ {2}}$ . To implikuje ${\ Displaystyle \ Gamma (S ^ {J}) ^ {2} \ równoważnik (1-1 / N) \ Gamma (S) ^ {2}}$ lub ${\ Displaystyle \ Gamma (S ^ {J}) \ równoważnik (1-1 / N) ^ {1/2} \ Gamma (S)}$ , czyli sekwencja obrotów Jacobiego jest zbieżna przynajmniej liniowo przez współczynnik ${\ Displaystyle (1-1 / N) ^ {1/2}}$ do przekątnej matrycy.

Liczba $obrotów$ nazywana jest przeciągnięciem; niech ${\ Displaystyle S ^ {\ sigma}}$ oznacza wynik. Poprzednie oszacowanie daje plony

${\ Displaystyle \ Gamma (S ^ {\ sigma}) \ równoważnik \ lewo (1- {\ Frac {1} {N}} \ prawo)^{N/2}\Gamma (S)}$ ,

tj. sekwencja przemiatań zbiega się co najmniej liniowo ze współczynnikiem ≈ ${\ displaystyle e ^ {1/2}}$ .

Jednak następujący wynik Schönhage daje lokalnie kwadratową zbieżność. W tym celu niech S ma m różnych wartości własnych ${\ Displaystyle \ lambda _ {1}, ..., \ lambda _ {m}}$ krotnościami ${\ Displaystyle \ nu _ {1}, ..., \ nu _ {m}}$ niech d > 0 będzie najmniejszą odległością dwóch różnych wartości własnych. Zadzwońmy po numer

${\ Displaystyle N_ {S}: = {\ Frac {n (n-1)} 2}}-\suma _{\mu =1}^{m}{\frac {1}{2}}\nu _{\mu }(\nu _{\mu }-1)\równoważnik N}$

Obroty Jacobiego i zamiatanie Schönhage. Jeśli ${\ displaystyle S ^ {s}}$ oznacza wynik wtedy

${\ Displaystyle \ Gamma (S ^ {s}) \ równoważnik {\ sqrt {{\ Frac {n} {2}}-1}}\left({\frac {\gamma ^{2}}{d-2\gamma }}\right),\quad \gamma :=\Gamma (S)}$ .

Zatem zbieżność staje się kwadratowa, gdy tylko ${\ Displaystyle \ Gamma (S) <{\ Frac {d} {2 + {\ sqrt {{\ Frac {n} {2} }-1}}}}}$

Koszt

  Każdy obrót Jacobiego można wykonać w krokach O ( n ), gdy znany jest element obrotu p . Jednak poszukiwanie p wymaga sprawdzenia wszystkich elementów poza przekątną N ≈ 1 / 2 n ^{2 .} Możemy to również zredukować do złożoności O ( n ), jeśli wprowadzimy dodatkową tablicę indeksów ${\ Displaystyle m_ {1}, \, \ kropki \,, \, m_ {n-1} }$ z właściwością ${\ Displaystyle m_ {i}}$ jest indeksem największego elementu w rzędzie ja , ( ja = 1, ..., n - 1) bieżącego S . Wtedy indeksy obrotu ( k , l ) muszą być jedną z par ${\ Displaystyle (i, m_ {i})}$ . Również aktualizacja tablicy indeksów może być wykonana w O( n ) złożoności średniego przypadku: Po pierwsze, maksymalny wpis w zaktualizowanych wierszach k i l można znaleźć w O ( n ) kroki. W pozostałych wierszach i zmieniają się tylko wpisy w kolumnach k i l . Zapętlając te wiersze, jeśli ${\ displaystyle m_ {i}$ jest ani k, ani $stare$ wystarczy porównać maksimum z wpisami i zaktualizować $} styl wyświetlania m_{i}}$ w razie potrzeby. Jeśli powinno być równe k lub . ${\ displaystyle m_ {i}}$ l i odpowiadający mu wpis zmniejszyły się podczas aktualizacji, maksymalne przekroczenie rzędu i należy znaleźć od zera w złożoności O( n ). Jednak stanie się to średnio tylko raz na obrót. ^obrót ma O( n ) i jedno przemiatanie O( n3 ) złożoność przypadku średniego, co odpowiada jednemu mnożeniu macierzy. Dodatkowo przed rozpoczęciem procesu należy zainicjować proces, co można $wykonać$ n ² .

Zazwyczaj metoda Jacobiego jest zbieżna w zakresie precyzji numerycznej po niewielkiej liczbie przeciągnięć. Zauważ, że wiele wartości własnych zmniejsza liczbę iteracji od ${\ displaystyle N_ {S}<N}$ .

Algorytm

Poniższy algorytm jest opisem metody Jacobiego w notacji matematycznej. Oblicza wektor e , który zawiera wartości $}$ i macierz E , która zawiera odpowiednie wektory własne, tj. ${i}$ własną, a kolumna ortonormalnym wektorem własnym dla mi ja {\ displaystyle $\ Displaystyle e_ {i}}$ , ja = 1, ..., n .

  
    
     procedura  jacobi(  S  ∈  R  n  ×  n  ;  out  e  ∈  R  n  ;  out  E  ∈  R  n  ×  n  )  var  ja  ,  k  ,  l  ,  m  ,  stan  ∈  N  s  ,  do  ,  t  ,  p  ,  y  ,  re  ,  r  ∈  R 
    
    

  
       
        
    
     ind  ∈  N  n  zmieniony  ∈  L  n  funkcja  maksind(  k  ∈  N  ) ∈  N  !  indeks największego elementu poza przekątną w rzędzie k  m  :=  k  +1  for  i  :=  k  +2  to  n  do  if  │  S  ki  │ > │  S  km  │  then  m  :=  i  endif  endfor  return  
  

  
    
       m  funkcja końcowa aktualizacja   procedury  (  k  ∈  N  ;  t  ∈  R  ) !  aktualizacja e  k  i jej status  y  :=  e  k  ;  e  k  :=  y  +  t  jeśli  zmieniono  k  i (  y  =  e  k  )  to  zmieniono  k  := fałsz;  stan  :=  stan  −1  elsif  (nie  
  

       zmieniono  k  ) i (  y  ≠  e  k  )  następnie  zmieniono  k  := prawda;  stan  :=  stan  +1  endif  procedura  endproc  obrót(  k  ,  l  ,  i  ,  j  ∈  N  ) !  │  wykonać obrót Sij  , S  kl  ┌  ┐ ┌ ┐┌  ┐ │  S  kl  │  c  −  s            
  
       ││  S  kl  │ │  │ := │ ││  │ │  S  ij  │ │  s  c  ││  Si  ij  │ └  ┘ └ ┘└  ┘  endproc  !  init e, E i tablice ind, zmienione  E  :=  I  ;  stan  :=  n  dla  k  := 1  do  n  do  ind  k  := maxind(  k  );  e  k  :=  S  kk 
   
    
      
        
    
     ;  zmieniono  k  := true  endfor  while  state  ≠0  do  !  następny obrót  m  := 1 !  znajdź indeks (k,l) obrotu p  dla  k  := 2  do  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  !  oblicz c = cos φ, s = sin φ  y  := (  mi  l  −  mi  k  )/2;  re  := │  y  │+√(  p  2  +  y  2  )  r  := √(  p  2  +  re  2  );  do  :=  re  /  r  ;  s  :=  p  /  r 
       
    
       ;  t  :=  p  2  /  re  jeśli  y  <0  to  s  := -  s  ;  t  := −  t  endif  S  kl  := 0,0; aktualizacja(   k  ,−  t  ); aktualizacja (   l  ,  t  )!  obróć wiersze i kolumny k i l  dla  i  := 1  do  k  −1  obróć  (  i  ,  k  ,  i  , 
      
       
       l  )  koniec  dla  i  :=  k  +1  do  l  −1  obróć  (  k  ,  ja  ,  ja  ,  l  )  koniecdla  i  :=  l  +1  do  n  obróć  (  k  ,  ja  ,  l  ,  ja  )  koniec  dla  !  obróć wektory własne  dla  i  := 1  do  n                 
         zrobić  ┌  ┐ ┌ ┐┌  ┐ │  E  ik  │ │  c  -  s  ││  E  ik  │ │  │ := │ ││  │ │  E  il  │ │  s  c  ││  E  il  │ └  ┘ └ ┘└  ┘  koniec dla  !  zaktualizuj wszystkie potencjalnie zmienione ind  i  dla  i  := 1  do  n  do  ind  i  := maxind ( 
  
 i  )  endfor  pętli  endproc 

Notatki

Zmieniona tablica logiczna zawiera status każdej wartości własnej. Jeśli wartość liczbowa $na$ się podczas iteracji, odpowiedni składnik zmieniony $ustawiany$ true , w przeciwnym razie . Stan całkowitoliczbowy zlicza liczbę elementów zmienionych , które mają wartość true . Iteracja zatrzymuje się, gdy tylko stan = 0. Oznacza to, że żadne z przybliżeń ${\ Displaystyle e_ {1}, \,... \,, e_ {n}}$ ostatnio zmienił swoją wartość i dlatego jest mało prawdopodobne, że tak się stanie, jeśli iteracja będzie kontynuowana. Tutaj zakłada się, że operacje zmiennoprzecinkowe są optymalnie zaokrąglane do najbliższej liczby zmiennoprzecinkowej.

2. Górny trójkąt macierzy S zostaje zniszczony, podczas gdy dolny trójkąt i przekątna pozostają niezmienione. W ten sposób możliwe jest przywrócenie S w razie potrzeby zgodnie z

  
       
        
    
 dla  k  := 1  do  n  −1  zrób  !  przywróć macierz S  dla  l  :=  k  +1  do  n  do  S  kl  :=  S  lk  endfor  endfor 

3. Wartości własne niekoniecznie są w kolejności malejącej. Można to osiągnąć za pomocą prostego algorytmu sortowania.

  
    
       
          
            
        
    
       dla  k  := 1  do  n  −1  do  m  :=  k  dla  l  :=  k  +1  do  n  do  if  e  l  >  e  m  then  m  :  =  l  endif  endfor  if  k  ≠  m  to  zamień  e  m  ,  ek  zamień  E  m  ,  E  k 
    
 endif  endfor 

4. Algorytm jest napisany przy użyciu notacji macierzowej (tablice oparte na 1 zamiast 0).

5. Podczas implementacji algorytmu część określona za pomocą notacji macierzowej musi być wykonywana jednocześnie.

6. Ta implementacja nie uwzględnia poprawnie przypadku, w którym jeden wymiar jest niezależną podprzestrzenią. Na przykład, jeśli otrzymamy macierz diagonalną, powyższa implementacja nigdy się nie zakończy, ponieważ żadna z wartości własnych nie ulegnie zmianie. Dlatego w rzeczywistych implementacjach należy dodać dodatkową logikę, aby uwzględnić ten przypadek.

Przykład

Niech ${\ Displaystyle S = {\ rozpocząć {pmatrix} 4 i -30 i 60 i -35 \\ -30 i 300 i - 675&420\\60&-675&1620&-1050\\-35&420&-1050&700\koniec{pmacierz}}}$

Następnie jacobi tworzy następujące wartości własne i wektory własne po 3 przeciągnięciach (19 iteracjach):

${\ Displaystyle e_ {1} = 2585,25381092892231}$

${\ Displaystyle E_ {1} = {\ rozpocząć {pmatrix} 0,029193 3231647860588\\-0,328712055763188997\\0,791411145833126331\\-0,514552749997152907\end{pmatrix}}}$

${\ Displaystyle e_ {2} = 37,1014913651276582}$

$-0,100228136947192199\\-0,638282528193614892\end{pmatrix}}}$

${\ Displaystyle e_ {3} = 1,4780548447781369}$

${\ Displaystyle E_ {3} = {\ rozpocząć {pmatrix} -0,5820756 99497237650\\0.370502185067093058\\0.509578634501799626\\0.514048272222164294\end{pmatrix}}}$

${\ Displaystyle e_ {4} = 0,1666428611718905}$

${\ Displaystyle E_ {4} = {\ rozpocząć {pmatrix} 0,792608291 163763585\\0.451923120901599794\\0.322416398581824992\\0.252161169688241933\end{pmatrix}}}$

Zastosowania rzeczywistych macierzy symetrycznych

Gdy znane są wartości własne (i wektory własne) macierzy symetrycznej, można łatwo obliczyć następujące wartości.

Wartości osobliwe

Wartości osobliwe macierzy (kwadratowej) A to pierwiastki kwadratowe z (nieujemnych) wartości własnych ${\ Displaystyle A ^ {T} A}$ . W przypadku macierzy symetrycznej $S$ mamy , stąd wartości S 2-

osobliwe

norma i promień widma

2-normy macierzy A $x \|_{2}=1}$ normą opartą na wektorze euklidesowym, tj. największej wartości, $gdy$ x przechodzi przez wszystkie wektory z . Jest to największa wartość osobliwa A . W przypadku macierzy symetrycznej jest to największa wartość bezwzględna jej wektorów własnych, a więc równa jej promieniu widmowemu.

Numer warunku

Numer warunku nieosobliwej macierzy A jest zdefiniowany jako ${\ Displaystyle {\ mbox {warunek}} (A) = \| A \ |_ {2} \| A ^ {-1} \|_ {2 }}$ . W przypadku macierzy symetrycznej jest to wartość bezwzględna ilorazu największej i najmniejszej wartości własnej. Macierze z dużymi liczbami warunków mogą powodować numerycznie niestabilne wyniki: małe zaburzenia mogą skutkować dużymi błędami. Macierze Hilberta są najbardziej znanymi macierzami źle uwarunkowanymi. Na przykład macierz Hilberta czwartego rzędu ma warunek 15514, podczas gdy dla rzędu 8 jest to 2,7 × 10 ⁸ .

Ranga

Macierz A ma rangę r , jeśli ma r kolumn, które są liniowo niezależne, podczas gdy pozostałe kolumny są od nich liniowo zależne. Równoważnie r jest wymiarem zakresu A . Ponadto jest to liczba niezerowych wartości osobliwych.

W przypadku macierzy symetrycznej r jest liczbą niezerowych wartości własnych. Niestety z powodu błędów zaokrągleń przybliżenia liczbowe zerowych wartości własnych mogą nie być zerowe (może się również zdarzyć, że przybliżenie liczbowe wynosi zero, podczas gdy wartość prawdziwa nie jest). liczbową można obliczyć tylko podejmując decyzję, które z wartości własnych są wystarczająco bliskie zeru.

Pseudoodwrotność

$AX$ macierzy A macierz dla i XA _ są symetryczne i dla których zachodzi AXA = A, XAX = X. Jeśli A nie jest liczbą pojedynczą, to ' ${\ Displaystyle A ^ {+} = A ^ {- 1}}$ . Gdy wywoływana jest procedura jacobi ( S ,

$= E ^ {T} {\ mbox {Diag}} (e)$ , E ), to relacja zachodzi tam, gdzie ( e ) oznacza macierz diagonalną z wektorem e na przekątnej. Niech ${\ Displaystyle e ^ {+}}$ oznacza wektor, w którym ${\ Displaystyle e_ {i}}$ jest zastępowany przez ${\ Displaystyle 1/e_ {i}}$ jeśli ${\ Displaystyle e_ {i}} } \ równoważnik 0}$ i przez 0, jeśli ${\ displaystyle e_ {i}}$ jest (liczbowo blisko) zera. Ponieważ macierz E jest ortogonalna, wynika z tego że pseudoodwrotność S jest dana wzorem ${\ Displaystyle S ^ {+} = E ^ {T} {\ mbox {Diag}} (e ^ {+}) E}$ .

Rozwiązanie metodą najmniejszych kwadratów

Jeśli macierz A nie ma pełnego rzędu, może nie istnieć rozwiązanie układu liniowego Ax = b . Można $,$ szukać minimalny Rozwiązaniem jest ${\ Displaystyle x = A ^ {+} b}$ . W przypadku symetrycznej macierzy S jak poprzednio, mamy ${\ Displaystyle x = S ^ {+} b = E ^ {T} {\ mbox {Diag}} (e ^ {+}) Eb}$ .

Macierz wykładnicza

Z ${\ Displaystyle S = E ^ {T} {\ mbox {Diag}} (e) E}$ znajduje się ${\ Displaystyle \ exp S = E ^ {T} {\ mbox {Diag}} (\ exp e) E}$ gdzie exp e jest wektorem, w którym ${\ Displaystyle e_ {i}}$ jest zastępowane przez ${\ Displaystyle \ exp e_ {i}}$ . W ten sam sposób f ( S ) można obliczyć w oczywisty sposób dla dowolnej (analitycznej) funkcji f .

Liniowe równania różniczkowe

Równanie różniczkowe x' = Ax , x (0) = a ma rozwiązanie x ( t ) = exp( t A ) a . Dla macierzy symetrycznej S wynika z tego, że ${\ Displaystyle x (t) = E ^ {T} {\ mbox {Diag}} (\ exp te) Ea}$ . Jeśli ${\ Displaystyle a = \ suma _ {i = 1} ^ {n} a_ {i} E_ {i}} jest$ rozwinięciem a przez wektory własne S , a następnie ${\ Displaystyle x (t) = \ suma _ {i = 1} ^ {n} a_ {i} \ exp (te_ {i}) E_ {i}}$ .

Niech będzie przestrzenią wektorową rozpiętą przez wektory własne S. $\ displaystyle W ^ {s}}$ które odpowiadają ujemnej wartości własnej i $.$ dodatnich wartości własnych za ${\ Displaystyle a \ w W ^ {s}},$$lim}} _ {t \ \ infty} x (t) = 0}$ ) tj. punkt równowagi 0 jest atrakcyjny dla x ( t ). Jeśli ${\ Displaystyle a \ w W ^ {u}},$ to ${\ Displaystyle {\ mbox {lim}} _ {t \ \ infty} x (t) = \ infty}$ , tj. 0 jest odpychające dla x ( t ). ${\ Displaystyle W ^ {s}$ } $i$ nazywane i niestabilnymi rozmaitościami dla S . Jeśli a ma składniki w obu rozmaitościach, to jeden składnik jest przyciągany, a drugi odpychany. Stąd x ( t ) zbliża się $\$ t $infty}$ .

Uogólnienia

Metoda Jacobiego została uogólniona na zespolone macierze hermitowskie , ogólne niesymetryczne macierze rzeczywiste i zespolone oraz macierze blokowe.

$można$ własnych macierzy symetrycznej, wykorzystać do obliczenia tych wartości. W tym przypadku metoda jest modyfikowana w taki sposób, że S nie musi być jawnie obliczane, co zmniejsza niebezpieczeństwo błędów zaokrągleń . Zauważ, że ${\ Displaystyle JSJ ^ {T} = JA ^ {T} AJ ^ {T} = JA ^ {T} J ^ {T} JAJ ^ {T} = B ^ {T$ z ${\ Displaystyle B \,: = JAJ ^ {T}}$ .

Metoda Jacobiego jest również dobrze przystosowana do równoległości.

Dalsza lektura

•   Prasa, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), „Sekcja 11.1. Transformacje Jacobiego macierzy symetrycznej” , Przepisy numeryczne: The Art of Scientific Computing (wyd. 3), New York: Cambridge University Press, ISBN 978-0-521-88068-8
•   Rutishauser, H. (1966). „Seria podręczników algebra liniowa: metoda Jacobiego dla rzeczywistych macierzy symetrycznych”. Matematyka numeryczna . 9 (1): 1–10. doi : 10.1007/BF02165223 . MR 1553948 .
•    Sameh, AH (1971). „O algorytmach Jacobiego i podobnych do Jacobiego dla komputera równoległego” . Matematyka obliczeń . 25 (115): 579–590. doi : 10.1090/s0025-5718-1971-0297131-6 . JSTOR 2005221 . MR 0297131 .
•    Shroff, Gautam M. (1991). „Algorytm równoległy dla wartości własnych i wektorów własnych ogólnej macierzy zespolonej”. Matematyka numeryczna . 58 (1): 779–805. CiteSeerX 10.1.1.134.3566 . doi : 10.1007/BF01385654 . MR 1098865 .
•   Veselić, K. (1979). „O klasie procedur podobnych do Jacobiego do diagonalizacji dowolnych macierzy rzeczywistych”. Matematyka numeryczna . 33 (2): 157–172. doi : 10.1007/BF01399551 . MR 0549446 .
•   Veselić, K.; Wenzel, HJ (1979). „Kwadratowo zbieżna metoda podobna do Jacobiego dla rzeczywistych macierzy o złożonych wartościach własnych”. Matematyka numeryczna . 33 (4): 425–435. doi : 10.1007/BF01399324 . MR 0553351 .
• Yousef Saad: „Revisiting the (block) Jacobi metoda rotacji podprzestrzeni dla symetrycznego problemu wartości własnej”, Numerical Algorithms, tom 92 (2023), s. 917-944. https://doi.org/10.1007/s11075-022-01377-w .

Linki zewnętrzne

• Implementacja w Matlabie algorytmu Jacobiego, który unika funkcji trygonometrycznych
• Implementacja C++11