Jacobi-módszer

A Jacobi-módszer (vagy Jacobi-féle sajátértékmódszer) néven ismert eljárás olyan iteratív módszer, amely kis méretű (n<10) szimmetrikus valós mátrixok sajátértékeinek és sajátvektorainak a meghatározására használható. Ezen módszer célja a mátrix főátlón kívüli elemeinek iteratív eljárással történő kinullázása. A Jacobi-módszer esetén az iterációs lépéseket addig ismételjük, míg egy általunk meghatározott pontosságig az ismeretleneket meg nem határozzuk. Ez azt fogja jelenteni, hogy akkor állunk meg a lépesekkel, mikor már két egymás utáni lépésben kapott ismeretlen értékek különbsége kisebb egy általunk meghatározott értéknél.

Nevét Carl Gustav Jacob Jacobiról kapta, aki először 1846-ban publikálta, de csak az 1950-es években vált elterjedtté a számítógépek fejlődése miatt.

A Jacobi-módszer esetében az iterációs képlet a következő lesz:

$x_{i}^{(k)}={b_{i} \over a_{ii}}-\sum _{j=1,j\neq i}^{n}a_{ij}x_{j}^{(k-1)}$

Ahhoz, hogy könnyebben megérthessük a módszer elvét, tekintsünk egy példát:

${\begin{pmatrix}a_{11}&a_{12}\\a_{21}&a_{22}\end{pmatrix}}{x_{1} \choose x_{2}}={b_{1} \choose b_{2}}$

Hogy jobban áttekinthető legyen, átírhatjuk egyenletek formájába, amely így nézhet ki:

${\begin{cases}a_{11}x_{1}+a_{12}x_{2}=b_{1}\\a_{21}x_{1}+a_{22}x_{2}=b_{2}\end{cases}}$

Innen kifejezhető az x₁ és x₂ ismeretlen, így a következő egyenleteket kapjuk:

$x_{1}=-{a_{12} \over a_{11}}x_{2}+{b_{1} \over a_{11}}$,

$x_{2}=-{a_{21} \over a_{22}}x_{1}+{b_{2} \over a_{22}}$

Az így kapott egyenletrendszert úgy oldhatjuk meg, hogy kezdetben kiindulunk az $x_{1}^{(0)}$, illetve az $x_{2}^{(0)}$ legjobb becslésünkből, vagy az egyszerűség kedvéért indulhatunk 0-ból is. Ezután felhasználva az

$x_{1}^{(k)}=-{a_{12} \over a_{11}}x_{2}^{(k-1)}+{b_{1} \over a_{11}}$,

$x_{2}^{(k)}=-{a_{21} \over a_{22}}x_{1}^{(k-1)}+{b_{2} \over a_{22}}$

lépéseket, eljuthatunk egy jobb közelítő értékig. Ezt addig alkalmazzuk, amíg az ismeretleneket tetszőleges pontossággal meg nem határozzuk.

Leírás

Az olyan transzformációt, ahol egy mátrixszal jobbról és az inverzével balról szorzunk egy mátrixot, hasonlósági transzformációnak nevezzük. A karakterisztikus egyenletet felírva belátható, hogy a hasonlósági transzformáció nem változtatja meg a sajátértékeket. Valós és szimmetrikus mátrixok esetén $\mathbf {R^{-1}=R^{T}}$, vagyis a hasonlósági transzformáció ortogonális transzformáció is egyben. Ezen az összefüggésen alapul a következőkben ismertetett módszer is. Vagyis a megfelelően megválasztott transzformációval a mátrixot diagonalizáljuk. Mivel a sajátvektorok maguk is valósak és ortogonálisak, az $\mathbf {A}$ szimmetrikus mátrix diagonalizálása megoldható az $\mathbf {R}$ ortogonális hasonlósági transzformáció segítségével, azaz

${\begin{aligned}\mathbf {R^{T}\cdot A\cdot R=\Lambda } .\end{aligned}}$

Vegyük példaként a $2\times 2$ típusú mátrix esetét. Ekkor a transzformációhoz használjuk a

$\mathbf {R} ={\begin{pmatrix}\cos \varphi &-\sin \varphi \\\sin \varphi &\cos \varphi \end{pmatrix}}$

síkforgatást leíró mátrixot, ahol $\varphi$ a forgatás szöge. Ha felírjuk ezzel az $\mathbf {A'=R^{T}\cdot A\cdot R}$ szimmetria transzformációt, a transzformálás után az $\mathbf {A'}$ mátrix elemei

${\begin{aligned}a_{11}^{'}&=a_{11}\cos ^{2}\varphi +2a_{21}\sin \varphi \cos \varphi +a_{22}\sin ^{2}\varphi \\a_{22}^{'}&=a_{11}\sin ^{2}\varphi -2a_{21}\sin \varphi \cos \varphi +a_{22}\cos ^{2}\varphi \\a_{21}^{'}&=a_{21}(\cos ^{2}\varphi -\sin ^{2}\varphi )+(a_{22}-a_{11})\sin \varphi \cos \varphi =a_{12}^{'}\end{aligned}}$

lesznek. Ha a nem átlós $a_{12}^{'}$ és $a_{21}^{'}$ elemeket 0-vá alakítjuk, az elforgatási szögre a következő egyenletet kapjuk:

${\begin{aligned}\cot ^{2}\varphi +{\frac {a_{22}-a_{11}}{a_{21}}}\cot \varphi -1&=0,\end{aligned}}$

melynek alapján

${\begin{aligned}\tan \varphi &={\Biggl [}{\frac {a_{11}-a_{22}}{2a_{12}}}\pm {\sqrt {{\Bigl (}{\frac {a_{11}-a_{22}}{2a_{12}}}{\Bigl )}^{2}+1}}{\Biggl ]}^{-1}.\end{aligned}}$

Innen megkaphatjuk a $\cos \varphi =(1+\tan \varphi )^{-1/2}$ és $\sin \varphi =\tan \varphi \cos \varphi$ függvényeket, melyekkel felépítjük a forgatásmátrixot. Az így kapott $\mathbf {A'}$ mátrix diagonális, tehát az átlóban található együtthatók a sajátértékek, míg az $\mathbf {R}$ forgatásmátrix két oszlopa a sajátértékeknek megfelelő két sajátvektor:

${\begin{aligned}\lambda _{1}&=a_{11}^{'},&\mathbf {x} ^{(1)}&={\begin{bmatrix}\cos \varphi \\\sin \varphi \end{bmatrix}};\\\lambda _{2}&=a_{22}^{'},&\mathbf {x} ^{(2)}&={\begin{bmatrix}-\sin \varphi \\\cos \varphi \end{bmatrix}}.\end{aligned}}$

Általános eset

A következőkben nézzük meg, hogy miként működik ez a módszer általános esetben $n\times n$ méretű mátrixok esetén. A sík-forgatás mátrixunk az egységmátrixtól csak az $r_{ii},r_{ij},r_{ji},r_{jj}$ elemekben tér el, vagyis

$\mathbf {R} _{ij}={\begin{pmatrix}1&\vdots &&\vdots &0\\\cdots &\cos \varphi &\cdots &-\sin \varphi &\cdots \\&\vdots &\ddots &\vdots &&\\\cdots &\sin \varphi &\cdots &\cos \varphi &\cdots \\0&\vdots &&\vdots &1\end{pmatrix}}$

Ezt felhasználva az

$\mathbf {A'=R_{ij}^{T}\cdot A\cdot R_{ij}}$

ortogonális hasonlósági transzformációval nullákat viszünk be az $a_{ij}^{'}$és $a_{ji}^{'}$ elemek helyére. A szorzás elvégzése után az

${\begin{aligned}a_{ik}^{'}&=a_{ki}^{'}=a_{ik}\cos \varphi +a_{jk}\sin \varphi ,k={\overline {1,n}}\\a_{jk}^{'}&=a_{kj}^{'}=a_{ik}\sin \varphi +a_{jk}\cos \varphi ,k\neq i,j\\a_{ii}^{'}&=a_{ii}\cos ^{2}\varphi +2a_{ji}\sin \varphi \cos \varphi +a_{jj}\sin ^{2}\varphi \\a_{jj}^{'}&=a_{ii}\sin ^{2}\varphi -2a_{ji}\sin \varphi \cos \varphi +a_{jj}\cos ^{2}\varphi \\a_{ji}^{'}&=a_{ji}(\cos ^{2}\varphi -\sin ^{2}\varphi )+(a_{jj}-a_{ii})\sin \varphi \cos \varphi =a_{ij}^{'}\end{aligned}}$

mátrixelemeket kapjuk eredményül. Ezek közül megköveteljük, hogy az $a_{ij}^{'}$ , illetve az $a_{ji}^{'}$ elemek 0-ák legyenek. Ekkor a

$\cot ^{2}\varphi +{\frac {a_{jj}-a_{ii}}{a_{ji}}}\cot \varphi -1=0$

egyenlethez jutunk, melyet megoldva a forgatás szöge

$\tan \varphi ={\Bigg [}{\frac {a_{ii}-a_{jj}}{2a_{ij}}}\pm {\sqrt {{\Big (}{\frac {a_{ii}-a_{jj}}{2a_{ji}}}{\Big )}^{2}+1}}{\Bigg ]}^{-1}$

lesz.

Meg kell jegyeznünk, hogy amikor egy másik elemet nullázunk ki a következő lépésben, akkor az előzőekben kinullázott elem elromlik. Viszont belátható, hogy bizonyos feltételek mellett az átlón kívüli elemek négyzetösszege egy lépésben $2\mid a_{ij}\mid ^{2}$-tel csökken, vagyis monoton módon tart 0-hoz.

$A_{l}$-lel jelölve az $l$. transzformáció utáni mátrixot, a transzformáció-sorozatot a következőképpen írhatjuk:

$\mathbf {A_{0}=A,A_{1}=R_{1}^{T}\cdot A_{0}\cdot R_{1},A_{2}=R_{2}^{T}\cdot A_{1}\cdot R_{2},\dots ,A_{l}=R_{l}^{T}\cdot A_{l-1}\cdot R_{l},\dots ,}$

ahol $\mathbf {R_{l}}$-el valamely nem-átlós elemre alkalmazott transzformációt jelöltük. Képezzük a transzformációs mátrixok

${\begin{aligned}\mathbf {X_{l}=R_{0}\cdot R_{1}\cdots R_{l}} ,&&l=0,1,2,\dots ,\end{aligned}}$

szorzatát. Ha végtelen sok transzformációt végzünk, akkor

${\begin{aligned}\lim _{l\to \infty }\mathbf {A_{l}=\Lambda } ,&&\lim _{l\to \infty }\mathbf {X_{l}=X} \end{aligned}}$

lesz. Ez azt jelenti, hogy ha ezeket a transzformációkat egymás után alkalmazzuk, akkor a mátrix diagonalizálódik, és az átlóban a sajátértékeket kapjuk. A sajátvektorok pedig a transzformációk szorzatmátrixának oszlopaiban lesznek.

A módszer konvergenciáját a $\tan \varphi \leq 1$ feltétel tiszteletben tartása biztosítja, ami egy $\varphi \leq \pi /4$ forgatásnak felel meg. Ezt úgy tudjuk biztosítani, hogy a két gyök közül a "+" előjelest választjuk, amennyiben $(a_{ii}-a_{jj})/a_{ji}>0$ és a "−" előjelest az ellenkező esetben. Ezt úgy tudjuk legkönnyebben megvalósítani, hogy a szöget a következőképpen számoljuk:

$\tan \varphi =sign{\Big (}{\frac {a_{ii}-a_{jj}}{2a_{ji}}}{\Big )}{\Bigg [}\left|{\frac {a_{ii}-a_{jj}}{2a_{ji}}}\right|+{\sqrt {{\Big (}{\frac {a_{ii}-a_{jj}}{2a_{ji}}}{\Big )}^{2}+1}}{\Bigg ]}^{-1}.$

Algoritmus

A leírt módszer a következő algoritmus segítségével alkalmazható számítógépre:

from __future__ import division
import math
dim=4

def Jacobi(a,imax,epsilon,x,l):
	for i in range(dim):
		for j in range(dim):
			x[i][j]=0
		x[i][i]=1
		l[i]=a[i][i]
	for it in range (imax):
		amax=0
		for j in range(1,dim,1):
			for i in range (j):
				a[i][i]=l[i]
				a[j][j]=l[j]
				a[j][i]=math.fabs(a[j][i])
				if amax<a[j][i]:
					amax=a[j][i]
				if a[j][i]>epsilon:
					tmp=(a[i][i]-a[j][j])/(2*a[j][i])
					t=1/(math.fabs(tmp)+math.sqrt(1+tmp*tmp))
					if tmp<0:
						t=-t
					c=1.0/(math.sqrt(1+t*t))
					s=c*t
					for k in range(i):
						temp=a[i][k]*c+a[j][k]*s
						a[j][k]=a[j][k]*c-a[i][k]*s
						a[i][k]=temp
					for k in range(i+1,j,1):
						temp=a[k][i]*c+a[j][k]*s
						a[j][k]=a[j][k]*c-a[k][i]*s
						a[k][i]=temp
					for k in range(j+1,dim,1):
						temp=a[k][i]*c+a[k][j]*s
						a[k][j]=a[k][j]*c-a[k][i]*s
						a[k][i]=temp
					for k in range (dim):
						temp=x[k][i]*c+x[k][j]*s
						x[k][j]=x[k][j]*c-x[k][i]*s
						x[k][i]=temp
					tmp=2*s*c*a[j][i]
					l[i]=a[i][i]*c*c+a[j][j]*s*s+tmp
					l[j]=a[i][i]*s*s+a[j][j]*c*c-tmp
					a[j][i]=0
		if amax<=epsilon:
			return 0
	return 666

a=[
	[3,0,2,1],
	[0,1,3,4],
	[2,3,2,1],
	[1,4,1,5]
	]
x=[
	[0,0,0,0],
	[0,0,0,0],
	[0,0,0,0],
	[0,0,0,0]
	]
l=[0,0,0,0]
epsilon=1e-16
imax=1e6
print a
b=Jacobi(a,imax,epsilon,x,l)
print b
print x
print l

Példa

Legyen $A={\begin{pmatrix}3&0&2&1\\0&1&3&4\\2&3&2&1\\1&4&1&5\end{pmatrix}}$

A jacobi a következő sajátértékeket és sajátvektorokat adja:

$\lambda _{1}=3.8614176875601696$

$\eta ^{(1)}={\begin{bmatrix}0.20847934594025172\\-0.9248091113211869\\-0.30940655891140356\\0.07437776036050801\end{bmatrix}}$

$\lambda _{2}=1.0818507981865024$

$\eta ^{(2)}={\begin{bmatrix}0.04702320745376396\\0.334174641833777\\-0.9043336166358197\\0.2613366345126636\end{bmatrix}}$

$\lambda _{3}=-2.741255286528889$

$\eta ^{(3)}={\begin{bmatrix}0.5641570498877014\\0.09332500337954595\\0.2860200054465712\\0.7689016993677129\end{bmatrix}}$

$\lambda _{4}=8.797986800782214$

$\eta ^{(4)}={\begin{bmatrix}-0.7975286849631742\\-0.156031599737972\\0.06812376684627766\\0.5787584029064224\end{bmatrix}}$

Források

• Numerikus módszerek, 1st, Kolozsvári Egyetemi Kiadó (2008) 
• Digitális tankönyvtár/Természettudományok/Matematika/Numerikus módszerek 1./Jacobi-módszer

További információk

• Press, WH; Teukolsky, SA & Vetterling, WT et al. (2007), "Section 11.1. Jacobi Transformations of a Symmetric Matrix", Numerical Recipes: The Art of Scientific Computing (3rd ed.), New York: Cambridge University Press, ISBN 978-0-521-88068-8
• Rutishauser, H. (1966). „Handbook Series Linear Algebra: The Jacobi method for real symmetric matrices.”. Numerische Mathematik 9 (1), 1–10. o. DOI:10.1007/BF02165223.  
• Sameh, A.H. (1971). „On Jacobi and Jacobi-like algorithms for a parallel computer”. Mathematics of Computation 25 (115), 579–590. o. DOI:10.1090/s0025-5718-1971-0297131-6.  
• Shroff, Gautam M. (1991). „A parallel algorithm for the eigenvalues and eigenvectors of a general complex matrix”. Numerische Mathematik 58 (1), 779–805. o. DOI:10.1007/BF01385654.  
• Veselić, K. (1979). „On a class of Jacobi-like procedures for diagonalising arbitrary real matrices”. Numerische Mathematik 33 (2), 157–172. o. DOI:10.1007/BF01399551.  
• Veselić, K. (1979). „A quadratically convergent Jacobi-like method for real matrices with complex eigenvalues”. Numerische Mathematik 33 (4), 425–435. o. DOI:10.1007/BF01399324.