Demo 4: Spektralsætningen, diagonalisering og hermitiske matricer#
Demo af Christian Mikkelstrup, Hans Henrik Hermansen, Magnus Troen, Karl Johan Måstrup Kristensen og Jakob Lemvig. Revideret marts 2026 af shsp.
from sympy import *
from dtumathtools import *
init_printing()
Symmetriske og hermitiske matricer#
En matrix kan undersøges for symmetriske og hermitiske egenskaber med Sympy på forskellig vis.
Med simuleret-manuelle beregninger#
A = Matrix([[6,2,4],[2,9,-2],[4,-2,6]])
A, A.T, A - A.T
B = Matrix([[6,2+7*I,4-4*I],[2-7*I,9,-2],[4+4*I,-2,6]])
B, B.adjoint(), B.H, (B - B.adjoint(), B - B.H) # A.H og A.adjoint() returnerer samme output
Bemærk, at da matricen konjugeres under operationen B.H, der producerer \(B^*\), da \(B^* = \bar B ^ T\), så vil en hermitisk matrix altid have reelle tal i diagonalen.
Med Sympy#
# For A
A.is_symmetric(), A.is_hermitian
(True, True)
# For B
B.is_symmetric(), B.is_hermitian
(False, True)
Vær opmærksom på, at .is_symmetric() er en Sympy-funktion, der skal tastes med afsluttende parenteser, hvorimod .is_hermitian er en egenskab, der bør tastes uden parenteser. Det er blot en særegenhed i Sympy; uden at gå for meget i dybden har det noget at gøre med, hvilken information der beregnes på forhånd og lagres af Sympy. Ved tvivl kan man altid bare prøve sig frem - udkommentér og kør fx følgende kodelinje:
# A.is_symmetric
Vi får så at vide, at vi forsøger at udføre en “bound method”, dvs. en Sympy/Python-funktion. Derfor skal parenteserne huskes:
A.is_symmetric()
True
Ved brug af definition#
Det er sjældent nødvendigt at benytte .is_symmetric(), ej heller at udføre manuelle beregninger, da det oftest er muligt direkte at se, om en matrix er symmetrisk. Et simpelt tjek af, om definitionen på en symmetrisk matrix, \(A = A^T\), er opfyldt, er ofte hurtigt og nemt:
A == A.T
True
Spektralsætningen#
Spektralsætningen er et godt eksempel på, hvorfor hermitiske (og reelt symmetriske) matricer er så utroligt brugbare, særligt i forbindelse med diagonalisering.
Spektral dekomposition af reel matrix#
Lad os som et eksempel betragte den reelle og symmetriske matrix:
A = Matrix([[1,2,0],
[2,3,2],
[0,2,1]])
A
Vi finder først dens egenværdier:
ev = A.eigenvects()
ev
Den har tre egenværdier, alle med \(\operatorname{am}(\lambda) = \operatorname{gm}(\lambda) = 1\), dvs. et én-dimensionelt egenrum for hver egenværdi. Spektralsætningen siger, at symmetriske matricer altid kan ortogonalt diagonaliseres ved:
hvor \(Q\) er en (reelt) ortogonal matrix bestående af egenvektorer fra \(A\) som søjler, og \(\Lambda\) er en diagonalmatrix indeholdende egenværdierne i diagonalen.
Da \(Q\) er (reelt) ortogonal, skal dens egenvektor-søjler være ortonormale. En sådan ortonormal mængde af søjler indsamler vi ved blot at normalisere en egenvektor for hver egenværdi. Spektralsætningen lover nemlig også, at forskellige egenrum er indbyrdes ortogonale, så normalisering er det eneste, der mangler:
Q = Matrix.hstack(*[ev[i][2][0].normalized() for i in range(3)])
Q
# Undgå at bruge lambda som variabelnavn, da det er et reserveret ord i Python. Vi bruger lamda i stedet
lamda = Q.T*A*Q
lamda
Vi har nu opnået ortogonal diagonalisering via en spektral dekomposition ved brug af en (reelt) ortogonal matrix.
Spektral dekomposition af kompleks matrix#
For matricer med komplekse elementer kan ovenstående tilgang følges, hvis vi blot arbejder på en normal matrix i stedet for en (reelt) ortogonal matrix. En normal matrix er en matrix, der opfylder:
Den spektrale dekomposition, vi søger, er da, for en unitær matrix \(U\):
Som før skal \(U\) bestå af en ortonormal mængde af egenvektorer fra \(A\) som søjler, mens \(\Lambda\) igen vil indeholde de tilsvarende egenværdier i diagonalen.
Lad os betragte følgende komplekse matrix som et eksempel:
A = Matrix([[1,2*I,0],
[-2*I,3,1+I],
[0,1-I,1]])
A
Er denne en normal matrix? Lad os tjekke efter:
A*A.adjoint(), A.adjoint()*A
Ja, det er den. Bemærk, at enhver hermitisk matrix også er normal, da den hermitiske betingelse \(A = A^*\) medfører \(AA^* = A^2 = A^*A\):
A.is_hermitian
True
For en spektral dekomposition løser vi først og fremmest egenværdiproblemet for \(A\):
ev = A.eigenvects()
ev
Igen har vi tre forskellige egenrum, så en egenvektor fra hvert rum vil give et mængde, der automatisk er indbyrdes ortogonal. Vi normaliserer og opnår de ortonormale søjler, som vores unitære matrix \(U\) skal bestå af:
U = Matrix.hstack(*[ev[i][2][0].normalized() for i in range(3)])
# Bekræfter at U er en unitær matrix
simplify(U*U.H), simplify(U.H*U)
Vores spektrale dekomposition er nu færdig:
simplify(U.H*A*U)
Diagonalisering via ortogonal substitution#
I tilfælde af ikke-én-dimensionelle egenrum#
Løber vi ind i tilfældet, hvor en egenværdi af en symmetrisk matrix har \(\operatorname{am}(\lambda) = \operatorname{gm}(\lambda) \ge 2\), så kan vi ikke være sikre på, at de egenværdier, vi finder, som udgangspunkt er ortogonale. Spektralsætningen sikrer dog fortsat, at en ortonormal basis bestående af egenvektorer eksisterer. Betragt for eksempel den symmetriske matrix \(A\) givet ved:
A = Matrix([[6,2,4],[2,9,-2],[4,-2,6]])
A
Den har en egenværdi med en algebraisk multiplicitet på 2 og to (lineært uafhængige) egenvektorer:
A.eigenvects()
Hvis vi som før blot normaliserer egenvektorerne:
v_1 = Matrix([-1,Rational(1,2),1])
v_2 = Matrix([Rational(1,2),1,0])
v_3 = Matrix([1,0,1])
[v.normalized() for v in [v_1,v_2,v_3]]
(alternativt i ét hug, uden vi behøver at indtaste \(\boldsymbol{v}\)-vektorerne manuelt):
V, Lamda = A.diagonalize()
[v_1,v_2,v_3] = [V.col(k) for k in range(3)] # Hver søjle i V trækkes ud og gemmes som en vektor
q_1, q_2, q_3 = [v.normalized() for v in [v_1, v_2, v_3]]
display([v_1,v_2,v_3])
display([q_1,q_2,q_3])
og derefter samler dem som søjler i en matrix \(Q\):
Q = Matrix.hstack(q_1,q_2,q_3)
Q
ser vi, at denne matrix desværre ikke er (reelt) ortogonal:
Q.T * Q
Problemet er, at de sidste to søjler i \(Q\) ikke er ortogonale.
Fiks problemet med Gram-Schmidt#
Problemet kan naturligvis løses med Gram-Schmidt-proceduren!
Bemærk, at vi kun har behov for Gram-Schmidt-proceduren indenfor de enkelte egenrum; altså, på egenvektorer fra samme egenrum, dvs. på egenvektorer, der deler en egenværdi. Og vi behøver ikke engang at udføre Gram-Schmidt-proceduren på én-dimensionelle egenrum, hvor \(\mathrm{am}(\lambda) = 1\), da egenvektorer herfra allerede er ortogonale på de andre vektorrum - disse behøver kun normalisering.
I vores eksempel kommer \(\boldsymbol v_1\) fra et egenrum med \(\mathrm{am}(\lambda) = 1\) og behøver derfor kun normalisering. \(\boldsymbol v_2\) og \(\boldsymbol v_3\) derimod kommer fra samme egenrum med \(\mathrm{am}(\lambda)=\mathrm{gm}(\lambda) = 2\), og vi udfører derfor Gram-Schmidt-proceduren på dem:
q_1 = v_1.normalized()
q_2, q_3 = GramSchmidt([v_2, v_3], orthonormal=True)
Q = Matrix.hstack(q_1,q_2,q_3)
Q
Tidligere udregnede vi en diagonalmatrix \(\Lambda\) ved udtrykket \(\Lambda=V^{-1} A V\), der gjorde brug af den ikke-ortogonale matrix \(V\), som indeholdt de dengang ikke-ortonormaliserede egenvektorer:
Lamda
Gram-Schmidt-proceduren bør ikke have ændret på \(\Lambda\). Et hurtigt tjek:
Lamda == Q.T*A*Q
True
Vi ser, som forventet, at \(\Lambda = V^{-1} A V = Q^T A Q\). Bemærk, at rækkefølgen af egenværdier i \(\Lambda\) vil stemme med rækkefølgen af deres tilsvarende egenvektorer som kolonner i matricen \(Q\), hhv. \(V\).
Positive (semi)definite matricer#
Positive (semi-)definite matricer har nogle egenskaber, der gør dem endnu nemmere at udføre beregninger med. De er herudover hermitiske og opfylder derfor alt, vi har arbejdet med i denne demo. Se definitionen af begrebet i lærebogen. Særligt i arbejde med matrixmanipulationer er sådanne matricer nemme, hvorfor vi sjældent lægger mærke til dem, hvis vi bruger Sympy/CAS. Sympy har dog to indbyggede funktioner, der gør det nemmere at finde ud af, om en matrix er positiv definit eller positiv semidefinit.
Betragt følgende to hermitiske matricer \(A,B \in \mathsf M_4(\mathbb{R})\):
A = Matrix([[5,8,4,1],
[8,17,8,2],
[4,8,5,0],
[1,2,0,1]])
B = Matrix([[1,-1,2,0],
[-1,4,-1,1],
[2,-1,6,-2],
[0,1,-2,4]])
A,B
Man tjekker, om de er positiv definite, med kommandoen:
A.is_positive_definite, B.is_positive_definite
(False, True)
Vi kan også forsøge at bevise, at \(B\) er positiv definit via simuleret-manuelle beregninger. En strategi kunne være at bevise aksiom (i) i Definition 2.9.1, men da det ville kræve, at vi tester for alle vektorer i \(\mathbb{R}^4\), er det ikke den bedste fremgangsmåde i Sympy. Derimod vil Theorem 2.9.1 (ii) være nemmere at arbejde med. Her skal der vises, at \(B\) har strengt positive egenværdier. Dette kan dog vise sig at være en udfordring for visse matricer:
B.eigenvals()
Selv hvis det lykkes dig at læse dette ganske lange output, vil det stadig ikke være klart, om egenværdierne er positive. De ser endda ud til at være ikke-reelle… Lad os se tallene som decimaltal:
for eigenval in B.eigenvals():
print(eigenval.evalf(15))
3.10533934717727 - 1.39585627777974e-29*I
0.0216077757381638 + 2.41351647874518e-30*I
8.26801827124053 + 3.58918706790694e-31*I
3.60503460584404 + 1.11861275922615e-29*I
Vi har her et mere brugbart udtryk, men egenværdierne ser stadig ikke-reelle ud. Vi kan nu fx argumentere for, at deres imaginærdele er så ubetydeligt små, at de kan ignoreres, som var de nul, hvilket ville være en acceptabel måde at komme videre på i praktisk brug. Og disse imaginærdele opstår da netop også kun pga. nogle numeriske afrundningsfejl - et fænomen, der kaldes “floating-point noise”. Men fra et matematisk perspektiv er dette stadigvæk hverken et særligt godt eller klart bevis.
Generelt set håber vi altid på, at vi er heldige at få med matricer med pæne egenværdier af gøre. Sympys indbyggede funktioner kan benyttes til tjek og kontrol, men ofte bør beviser udføres i hånden, før man kan være helt sikker.