Demo 4: The Spectral Theorem, Diagonalization, and Hermitian Matrices#
Demo by Christian Mikkelstrup, Hans Henrik Hermansen, Magnus Troen, Karl Johan Måstrup Kristensen, and Jakob Lemvig. Revised March 2026 by shsp.
from sympy import *
from dtumathtools import *
init_printing()
Symmetric and Hermitian Matrices#
According to the textbook, a matrix \(M\) is symmetric, respectively Hermitian, if it is equal to its own transpose \(M=M^T\), respectively adjoint \(M=M^*\). We can investigate this in Sympy in different ways. As some examples, let’s investigate for these two matrices:
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
Using Simulated-Manual Calculations#
A.T, A - A.T
B.adjoint(), B - B.adjoint()
Note that an alternative to .adjoint() with Python is .H:
B.H, B - B.H
We could also instead have done a logical comparison:
A == A.T, B == B.H
(True, True)
Using Sympy#
A.is_symmetric(), A.is_hermitian
(True, True)
As expected, a (real) symmetric matrix is also Hermitian.
B.is_symmetric(), B.is_hermitian
(False, True)
Be aware that .is_symmetric() is a Sympy function and must be typed with brackets afterwards, whereas .is_hermitian is a property and should be used without brackets. This is a peculiarity in Sympy. Without going too deep into the details, this has got something to do with what information that is being computed beforehand and stored by Sympy and what that is not. When in doubt, just give it a try - e.g., try uncommenting this code line:
# A.is_symmetric
Here we are being told that we are trying to perform a “bound method”, meaning a Sympy/Python function. Thus brackets must be added.
The Spectral Theorem#
According to the textook, a complex normal matrix \(M\in\mathsf M_n(\Bbb C)\), meaning a matrix that fulfills \(MM^* = M^*M\), can always be unitarily diagonalized by:
where \(U\) is a unitary matrix consisting of orthonormal eigenvectors from \(M\) as columns, and where \(\Lambda\) contains the corresponding eigenvalues in its diagonal. In the real case of \(M\in\mathsf M_n(\Bbb R)\), then \(M\) just has to be symmetric, so fulfilling \(M=M^T\), for us to be able to do orthogonal diagonalization by:
\(Q\) is here a (real) orthogonal matrix. All this is a part of the famous spectral theorem.
Spectral Decomposition of Real Matrix#
As an example, consider the real and symmetric matrix:
A = Matrix([[1,2,0],
[2,3,2],
[0,2,1]])
A
For a spectral decomposition we needs its eigenvalues for \(\Lambda\) and its eigenvectors for \(Q\):
ev = A.eigenvects()
ev
It has three eigenvalues, all with \(\operatorname{am}(\lambda) = \operatorname{gm}(\lambda) = 1\), meaning a one-dimensional eigenspace for each eigenvalue. The (real) orthogonal \(Q\) must have orthonormal eigenvector columns, which we here achieve with a simple normalization of an eigenvector for each eigenvalue. The spectral theorem promises that different eigenspaces are mutually orthogonal, why normalization is the only thing missing:
Q = Matrix.hstack(*[ev[i][2][0].normalized() for i in range(3)])
Q
# Avoid using lambda as variable name since it is a reserved word in Python. We will use lamda instead
lamda = Q.T*A*Q
lamda
We have now achieved orthogonal diagonalization by carrying out a spectral decomposition using a (real) orthogonal matrix.
Spectral Decomposition of Complex Matrix#
Now for an example of a complex matrix:
B = Matrix([[1,2*I,0],
[-2*I,3,1+I],
[0,1-I,1]])
B
Is this a normal matrix? Let’s check:
B*B.adjoint(), B.adjoint()*B
Yes, it is. Note that all Hermitian matrices automatically are normal, since the Hermitian condition \(B = B^*\) implies \(BB^* = B^2 = B^*B\):
B.is_hermitian
True
For the spectral decomposition, let’s first solve the eigenvalue-problem:
ev = B.eigenvects()
ev
Again we see three different eigenspaces, so an eigenvector from each will automatically be orthogonal. We normalize and achieve the orthonormal columns for our unitary matrix \(U\):
U = Matrix.hstack(*[ev[i][2][0].normalized() for i in range(3)])
# U is a unitary matrix
simplify(U*U.H), simplify(U.H*U)
Our spectral decomposition is now complete:
simplify(U.H*B*U)
Diagonalization via Orthogonal Substitution#
Consider the following symmetric and real matrix \(C\):
C = Matrix([[6,2,4],[2,9,-2],[4,-2,6]])
C
We solve the eigenvalue-problem:
C.eigenvects()
and do a quick diagonalization \(\Lambda = V^{-1}CV\) with the three found linearly independent eigenvectors as columns in \(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
When Not All Eigenspaces are One-Dimensional#
Let’s now do an orthogonal diagonalization of \(C\) via the spectral theorem as in the earlier sections. But this time we must be careful! One of the eigenvalues has \(\operatorname{am}(\lambda) = \operatorname{gm}(\lambda) \ge 2\), meaning a two-dimensional eigenspace. .eigenvects() has given us three linearly independent eigenvectors, but two of them come from the same eigenspace. The spectral theorem only ensures that eigenvectors from different eigenspaces are orthogonal, and we cannot assume that two vectors from the same eigenspace are orthogonal.
If we as before simply normalize all three eigenvectors:
q1, q2, q3 = [v.normalized() for v in [v1, v2, v3]]
q1, q2, q3
(or in one go without having to type the \(\boldsymbol{v}\) vectors manually):
V, Lamda = C.diagonalize()
[v1,v2,v3] = [V.col(k) for k in range(3)] # Each column in V is retrieved and stored as a vector
q1, q2, q3 = [v.normalized() for v in [v1, v2, v3]]
display([v1,v2,v3])
display([q1,q2,q3])
and then form the indended \(Q\) matrix from those as columns:
Q = Matrix.hstack(q1,q2,q3)
Q
we see that this matrix indeed is not (real) orthogonal:
Q.T * Q
Fixing the Issue with the Gram-Schmidt Procedure#
The problem can of course be fixed by the Gram-Schmidt-procedure! And it is simpler than we might fear, as Gram-Schmidt only is needed within each eigenspace; that is, only on eigenvectors from the same eigenspace.
In our example, \(\boldsymbol v_1\) comes from an eigenspace with \(\mathrm{am}(\lambda) = 1\) and hence only needs normalization. \(\boldsymbol v_2\) and \(\boldsymbol v_3\) come from an eigenspace with \(\mathrm{am}(\lambda)=\mathrm{gm}(\lambda) = 2\) and we will do the Gram-Schmidt procedure on those:
q1 = v1.normalized()
q2, q3 = GramSchmidt([v2, v3], orthonormal=True)
Q = Matrix.hstack(q1,q2,q3)
Q
Earlier we found a diagonal matrix \(\Lambda\) as \(\Lambda=V^{-1} C V\) using the non-orthogonal matrix \(V\) that contained the at that time non-orthonormalized eigenvectors. The Gram-Schmidt procedure should not have changed \(\Lambda\), so we will expect that \(\Lambda = V^{-1} C V = Q^T C Q\). A quick check:
Lamda == Q.T*C*Q
True
As expected. Note that the order of eigenvalues in \(\Lambda\) will match the order of their corresponding eigenvectors as columns in the matrix \(Q\), respectively \(V\) - so, if you have altered the order, this logical comparison won’t hold.
Positive (Semi-)Definite Matrices#
Positive (semi-)definite matrices have some properties that make them particularly easier to perform calculations with. They are furthermore Hermitian and hence fulfill everything we have worked through in this demo. See the definitions in the textbook. You will mainly notice this during matrix manipulations, and this is thus not something we necessarily notice when using Sympy/CAS. Sympy is, though, constructed with two functions that make it easy to find out whether a matrix is positive definite or positive semi-definite.
Consider the following two Hermitian matrices \(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
We can check whether they are positive definite with this command:
A.is_positive_definite, B.is_positive_definite
(False, True)
We could also try proving that \(B\) is positive definite via simulated-manual calculation. A strategy could be to prove axiom (i) in Definition 2.9.1, but since this would require that we test for all vectors in \(\mathbb{R}^4\), this might not be the best choice in Sympy. Instead, Theorem 2.9.1 (ii) would be easier to work with. Here we want to show that \(B\) has strictly positive eigenvalues. This turns out to be a challenge, though:
B.eigenvals()
If you manage to read this very long output, you’ll find that it still isn’t obvious whether these eigenvalues should be positive. They even seem to be non-real… Let us have a look at them as decimal numbers:
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
We now see a more useful output, but the eigenvalues still appear to be non-real. We might argue in this case that the imaginary parts are so insignificantly tiny that they can be assumed to be zero, which we would be happy with to all practical uses. Indeed those imaginary parts only appear due to some numerical rounding errors - it can be called “floating-point noise”. Even so, this is neither a particularly good nor clear mathematical proof.
Generally speaking, we always hope to be fortunate enough to be dealing with matrices with nicer eigenvalues. Sympy’s functions can be used for checks and control, but you should always carry out proofs for certainty by hand.