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

from sympy import *
from dtumathtools import *
init_printing()

Symmetric and Hermitian Matrices

In Sympy, matrices can be checked for symmetry and Hermitian properties as follows.

1. You can use simulated-manual calculations:

A = Matrix([[6,2,4],[2,9,-2],[4,-2,6]])
A, A.T, A - A.T 

\[\begin{split}\displaystyle \left( \left[\begin{matrix}6 & 2 & 4\\2 & 9 & -2\\4 & -2 & 6\end{matrix}\right], \ \left[\begin{matrix}6 & 2 & 4\\2 & 9 & -2\\4 & -2 & 6\end{matrix}\right], \ \left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]\right)\end{split}\]

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 and A.adjoint() return the same output in python

\[\begin{split}\displaystyle \left( \left[\begin{matrix}6 & 2 + 7 i & 4 - 4 i\\2 - 7 i & 9 & -2\\4 + 4 i & -2 & 6\end{matrix}\right], \ \left[\begin{matrix}6 & 2 + 7 i & 4 - 4 i\\2 - 7 i & 9 & -2\\4 + 4 i & -2 & 6\end{matrix}\right], \ \left[\begin{matrix}6 & 2 + 7 i & 4 - 4 i\\2 - 7 i & 9 & -2\\4 + 4 i & -2 & 6\end{matrix}\right], \ \left( \left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right], \ \left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]\right)\right)\end{split}\]

Note that since the matrix is being conjugated in the operation B.H, which produces \(B^*\), since \(B^* = \bar B ^ T\), a Hermitian matrix cannot have non-real numbers in the diagonal.

2. Alternatively, you can use a Sympy command:

# For A
A.is_symmetric(), A.is_hermitian

(True, True)

# For B
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 tiny peculiarity in Sympy. Without going too much into the details, this has got something to do with which information that has been computed beforehand and stored by Sympy and which that have not. When in doubt, one can always just give it a try:

# A.is_symmetric

Here we are being told that we are trying to perform a “bound method”, meaning a function. Thus brackets must be added:

A.is_symmetric()

True

It is rarely necessary to use .is_symmetric(), since it often is possible to clarify whether a matrix is symmetric simply by looking at it. Finally, one can check whether the definition of a symmetric matrix is directly fulfilled, meaning whether \(A = A^T\) is true. In Python:

A == A.T

True

The Spectral Theorem

The spectral theorem provides a good example of why Hermitian (and real symmetric) matrices are so unbelievably practical. In particular in relation to diagonalization.

Consider for example the real and symmetric matrix

A = Matrix([[1,2,0],
            [2,3,2],
            [0,2,1]])
A

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 2 & 0\\2 & 3 & 2\\0 & 2 & 1\end{matrix}\right]\end{split}\]

First we find the eigenvectors:

ev = A.eigenvects()
ev

\[\begin{split}\displaystyle \left[ \left( -1, \ 1, \ \left[ \left[\begin{matrix}1\\-1\\1\end{matrix}\right]\right]\right), \ \left( 1, \ 1, \ \left[ \left[\begin{matrix}-1\\0\\1\end{matrix}\right]\right]\right), \ \left( 5, \ 1, \ \left[ \left[\begin{matrix}1\\2\\1\end{matrix}\right]\right]\right)\right]\end{split}\]

Since we have 3 eigenvalues that all fulfill \(\operatorname{am}(\lambda) = \operatorname{gm}(\lambda) = 1\) we can use Corollary 2.8.6 to easily diagonalize \(A\) by doing

\[\begin{equation*} \Lambda = Q^T A Q. \end{equation*}\]

where \(Q\) is an orthogonal matrix, consisting of eigenvectors of \(A\).

Q = Matrix.hstack(*[ev[i][2][0].normalized() for i in range(3)])
Q

\[\begin{split}\displaystyle \left[\begin{matrix}\frac{\sqrt{3}}{3} & - \frac{\sqrt{2}}{2} & \frac{\sqrt{6}}{6}\\- \frac{\sqrt{3}}{3} & 0 & \frac{\sqrt{6}}{3}\\\frac{\sqrt{3}}{3} & \frac{\sqrt{2}}{2} & \frac{\sqrt{6}}{6}\end{matrix}\right]\end{split}\]

# 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

\[\begin{split}\displaystyle \left[\begin{matrix}-1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 5\end{matrix}\right]\end{split}\]

For matrices with complex values the same approach can be used, as long as we work with normal matrices. Meaning, matrices that fulfill

\[\begin{equation*} AA^* = A^*A. \end{equation*}\]

Note that all Hermitian matrices also are normal since

\[\begin{equation*} A = A^* \end{equation*}\]

implies

\[\begin{equation*} AA^* = A^2 = A^*A. \end{equation*}\]

See for example the matrix:

A = Matrix([[1,2*I,0],
           [-2*I,3,1+I],
           [0,1-I,1]])

display(A, A.is_hermitian)

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 2 i & 0\\- 2 i & 3 & 1 + i\\0 & 1 - i & 1\end{matrix}\right]\end{split}\]

True

# Eigenvectors are found
ev = A.eigenvects()
ev

\[\begin{split}\displaystyle \left[ \left( 1, \ 1, \ \left[ \left[\begin{matrix}\frac{1}{2} - \frac{i}{2}\\0\\1\end{matrix}\right]\right]\right), \ \left( 2 - \sqrt{7}, \ 1, \ \left[ \left[\begin{matrix}-1 + i\\- \frac{1}{2} - \frac{i}{2} + \left(\frac{1}{2} + \frac{i}{2}\right) \left(2 - \sqrt{7}\right)\\1\end{matrix}\right]\right]\right), \ \left( 2 + \sqrt{7}, \ 1, \ \left[ \left[\begin{matrix}-1 + i\\- \frac{1}{2} - \frac{i}{2} + \left(\frac{1}{2} + \frac{i}{2}\right) \left(2 + \sqrt{7}\right)\\1\end{matrix}\right]\right]\right)\right]\end{split}\]

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)

\[\begin{split}\displaystyle \left( \left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right], \ \left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]\right)\end{split}\]

simplify(U.H*A*U)

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 2 - \sqrt{7} & 0\\0 & 0 & 2 + \sqrt{7}\end{matrix}\right]\end{split}\]

Diagonalization of a Symmetric Matrix via Orthogonal Substitution

If we run into the case where an eigenvalue of a symmetric matrix has \(\operatorname{am}(\lambda) = \operatorname{gm}(\lambda) \ge 2\), we can be sure that the first eigenvectors we find are orthogonal. The spectral theorem still ensures that an orthonormal basis consisting of eigenvectors exists. Consider for instance the symmetric matrix \(A\) given by

A = Matrix([[6,2,4],[2,9,-2],[4,-2,6]])
A

\[\begin{split}\displaystyle \left[\begin{matrix}6 & 2 & 4\\2 & 9 & -2\\4 & -2 & 6\end{matrix}\right]\end{split}\]

Here we see that an eigenvalue has an algebraic multiplicity of 2 and two (linearly independent) eigenvectors

A.eigenvects()

\[\begin{split}\displaystyle \left[ \left( 1, \ 1, \ \left[ \left[\begin{matrix}-1\\\frac{1}{2}\\1\end{matrix}\right]\right]\right), \ \left( 10, \ 2, \ \left[ \left[\begin{matrix}\frac{1}{2}\\1\\0\end{matrix}\right], \ \left[\begin{matrix}1\\0\\1\end{matrix}\right]\right]\right)\right]\end{split}\]

So if we, as before, only normalize the eigenvectors:

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]]

\[\begin{split}\displaystyle \left[ \left[\begin{matrix}- \frac{2}{3}\\\frac{1}{3}\\\frac{2}{3}\end{matrix}\right], \ \left[\begin{matrix}\frac{\sqrt{5}}{5}\\\frac{2 \sqrt{5}}{5}\\0\end{matrix}\right], \ \left[\begin{matrix}\frac{\sqrt{2}}{2}\\0\\\frac{\sqrt{2}}{2}\end{matrix}\right]\right]\end{split}\]

(or in one go without having to type the \(\boldsymbol{v}\) vectors manually):

V, Lamda = A.diagonalize()
[v_1,v_2,v_3] = [V.col(k) for k in range(3)]        # Each column in V is retrieved and stored as a vector
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])

\[\begin{split}\displaystyle \left[ \left[\begin{matrix}-2\\1\\2\end{matrix}\right], \ \left[\begin{matrix}1\\2\\0\end{matrix}\right], \ \left[\begin{matrix}1\\0\\1\end{matrix}\right]\right]\end{split}\]

\[\begin{split}\displaystyle \left[ \left[\begin{matrix}- \frac{2}{3}\\\frac{1}{3}\\\frac{2}{3}\end{matrix}\right], \ \left[\begin{matrix}\frac{\sqrt{5}}{5}\\\frac{2 \sqrt{5}}{5}\\0\end{matrix}\right], \ \left[\begin{matrix}\frac{\sqrt{2}}{2}\\0\\\frac{\sqrt{2}}{2}\end{matrix}\right]\right]\end{split}\]

and create the wanted \(Q\) matrix with those as columns

Q = Matrix.hstack(q_1,q_2,q_3)
Q

\[\begin{split}\displaystyle \left[\begin{matrix}- \frac{2}{3} & \frac{\sqrt{5}}{5} & \frac{\sqrt{2}}{2}\\\frac{1}{3} & \frac{2 \sqrt{5}}{5} & 0\\\frac{2}{3} & 0 & \frac{\sqrt{2}}{2}\end{matrix}\right]\end{split}\]

we see that the matrix unfortunately is not real orthogonal,

Q.T * Q

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 1 & \frac{\sqrt{10}}{10}\\0 & \frac{\sqrt{10}}{10} & 1\end{matrix}\right]\end{split}\]

since the last two columns in \(Q\) are not orthogonal to each other!

Fortunately we can use a method we have learned in the course for these cases, that being the Gram-Schmidt procedure. For eigenvectors with equal eigenvalues the Gram-Schmidt procedure is carried out, and eigenvectors for which \(\mathrm{am}(\lambda) = 1\) are just treated with normalization:

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

\[\begin{split}\displaystyle \left[\begin{matrix}- \frac{2}{3} & \frac{\sqrt{5}}{5} & \frac{4 \sqrt{5}}{15}\\\frac{1}{3} & \frac{2 \sqrt{5}}{5} & - \frac{2 \sqrt{5}}{15}\\\frac{2}{3} & 0 & \frac{\sqrt{5}}{3}\end{matrix}\right]\end{split}\]

The order of eigenvalues in the diagonal matrix \(\Lambda\) is determined by the order in which the eigenvectors we chose are merged into a matrix. \(\Lambda\) is equal to \(V^{-1} A V\) and given by:

Lamda

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 10 & 0\\0 & 0 & 10\end{matrix}\right]\end{split}\]

Let us check that the Gram-Schmidt procedure has not changed \(\Lambda\). We write:

Lamda == Q.T*A*Q

True

As expected, we have \(\Lambda = V^{-1} A V = Q^T A Q\). But remember that this Gram-Schmidt method only is necessary when the matrix satisfies the spectral theorem (again for real matrices according to Theorem 2.8.5 and for complex matrices according to Theorem 2.8.9 in the textbook) but has eigenvalues with \(\operatorname{am}(\lambda) = \operatorname{gm}(\lambda) \ge 2\).

Positive (Semi-)Definite Matrices

Positive (semi-)definite matrices have some properties that make them even easier to perform calculations with. They are furthermore Hermitian and hence fulfill everything we have worked through in this demo. 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 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

\[\begin{split}\displaystyle \left( \left[\begin{matrix}5 & 8 & 4 & 1\\8 & 17 & 8 & 2\\4 & 8 & 5 & 0\\1 & 2 & 0 & 1\end{matrix}\right], \ \left[\begin{matrix}1 & -1 & 2 & 0\\-1 & 4 & -1 & 1\\2 & -1 & 6 & -2\\0 & 1 & -2 & 4\end{matrix}\right]\right)\end{split}\]

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. So, we want to show that \(B\) has strictly positive eigenvalues. This turns out to be a challenging case, though:

B.eigenvals()

It is not particularly obvious that these eigenvalues should be positive, and they even seem to be complex. We can instead try looking at the numerical values:

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

Here one could argue that the imaginary parts are so insignificantly tiny that they can be assumed to be zero, which we would be happy with. This is, though, neither a particularly good nor clear proof.

One can be fortunate and be dealing with a matrix with nicer eigenvalues. Sympy’s functions can be used for checks and control, but we recommend that proofs are carried out by hand instead.