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#
Jf. lærebogen er en matrix \(M\) symmetrisk, hhv. hermitisk, hvis den er lig med sin egen transponerede \(M=M^T\), hhv. adjungerede \(M=M^*\). Dette kan undersøges med Sympy på forskellig vis. Lad os fx undersøge følgende to matricer:
A = Matrix([[6,2,4],[2,9,-2],[4,-2,6]])
B = Matrix([[6,2+7*I,4-4*I],[2-7*I,9,-2],[4+4*I,-2,6]])
A,B
Simuleret-manuelt#
A.T, A - A.T
B.adjoint(), B - B.adjoint()
Bemærk, at et alternativ til .adjoint(), når man skal bestemme den adjungerede med Python, er .H:
B.H, B - B.H
Man kunne også i stedet have udført en logisk sammenligning:
A == A.T, B == B.H
(True, True)
Med Sympy#
A.is_symmetric(), A.is_hermitian
(True, True)
Som forventet er en (reel) symmetrisk matrix også hermitisk.
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 parenteser huskes.
Spektralsætningen#
Jf. lærebogen kan en komplex normal matrix \(M\in\mathsf M_n(\Bbb C)\), altså en matrix, der opfylder \(MM^* = M^*M\), altid unitært diagonaliseres ved:
hvor \(U\) er en unitær matrix bestående af ortonormale egenvektorer fra \(M\) som søjler, og hvor \(\Lambda\) indeholder de tilhørende egenværdier i diagonalen. I det reelle tilfælde \(M\in\mathsf M_n(\Bbb R)\), da skal \(M\) blot være symmetrisk, altså opfylde \(M=M^T\), før vi kan udføre ortogonal diagonalisering ved:
\(Q\) er her en (reelt) ortogonal matrix. Alt dette kommer af den berømte spektralsætning.
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
Til en spektral dekomposition skal vi bruge dens egenværdier til \(\Lambda\) og dens egenvektorer til \(Q\):
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. Den (reelt) ortogonale \(Q\) skal have ortonormale egenvektor-søjler, hvilket vi her opnår ved simpelthen at normalisere en egenvektor for hver egenværdi. Spektralsætningen lover nemlig, 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#
Lad os betragte følgende komplekse matrix:
B = Matrix([[1,2*I,0],
[-2*I,3,1+I],
[0,1-I,1]])
B
Er den normal? Lad os undersøge:
B*B.adjoint(), B.adjoint()*B
Ja, det er den. Bemærk, at enhver hermitisk matrix også er normal, da den hermitiske betingelse \(B = B^*\) medfører \(BB^* = B^2 = B^*B\):
B.is_hermitian
True
For en spektral dekomposition løser vi først og fremmest egenværdiproblemet:
ev = B.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*B*U)
Diagonalisering via ortogonal substitution#
Betragt følgende symmetriske og reelle matrix \(C\):
C = Matrix([[6,2,4],[2,9,-2],[4,-2,6]])
C
Vi løser egenværdiproblemet:
C.eigenvects()
og udfører en hurtig diagonalisering \(\Lambda = V^{-1}CV\) med de tre fundne lineært uafhængige egenvektorer som søjler i \(V\):
v1 = Matrix([-1,Rational(1,2),1])
v2 = Matrix([Rational(1,2),1,0])
v3 = Matrix([1,0,1])
V = Matrix.hstack(v1,v2,v3)
Lamda = diag(1, 10, 10)
V, Lamda
Problemet med flerdimensionelle egenrum#
Lad os udføre en ortogonal diagonalisering af \(C\) via spektralsætningen som i de tidligere afsnit. Men denne gang skal vi passe på! For én af egenværdierne har \(\operatorname{am}(\lambda) = \operatorname{gm}(\lambda) \ge 2\), altså et egenrum er to-dimensionelt. .eigenvects() har givet os tre lineært uafhængige egenvektorer, men to af dem kommer fra samme egenrum. Spektralsætningen sikrer kun, at egenvektorer fra forskellige egenrum er ortogonale, og vi kan ikke regne med, at to vektorer fra samme egenrum er ortogonale.
Hvis vi som før blot normaliserer alle tre egenvektorer:
q1, q2, q3 = [v.normalized() for v in [v1, v2, v3]]
(alternativt i ét hug, uden vi behøver at indtaste \(\boldsymbol{v}\)-vektorerne manuelt):
V, Lamda = C.diagonalize()
[v1,v2,v3] = [V.col(k) for k in range(3)] # Hver søjle i V trækkes ud og gemmes som en vektor
q1, q2, q3 = [v.normalized() for v in [v1, v2, v3]]
display([v1,v2,v3])
display([q1,q2,q3])
og derefter samler dem som søjler i en matrix \(Q\):
Q = Matrix.hstack(q1,q2,q3)
Q
ser vi, at denne matrix rigtig nok ikke er (reelt) ortogonal:
Q.T * Q
Løs problemet med Gram-Schmidt#
Problemet kan naturligvis løses med Gram-Schmidt-proceduren! Og det bliver simplere, end man kunne frygte, da der kun er behov for Gram-Schmidt indenfor de enkelte egenrum; altså, på egenvektorer fra samme egenrum.
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:
q1 = v1.normalized()
q2, q3 = GramSchmidt([v2, v3], orthonormal=True)
Q = Matrix.hstack(q1,q2,q3)
Q
Tidligere fandt vi en diagonalmatrix \(\Lambda\) ved udtrykket \(\Lambda=V^{-1} C V\), der gjorde brug af den ikke-ortogonale matrix \(V\), som indeholdt de dengang ikke-ortonormaliserede egenvektorer. Gram-Schmidt-proceduren bør ikke have ændret på \(\Lambda\), så vi vil forvente \(\Lambda = V^{-1} C V = Q^T C Q\). Et hurtigt tjek:
Lamda == Q.T*C*Q
True
Som forventet. 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\) - så, hvis du har ændret i rækkefølgen, vil denne logiske sammenligning ikke holde.
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.