import numpy as np
import sympy as sy
sy.init_printing()
import matplotlib.pyplot as plt
plt.style.use('ggplot')
If $A = PBP^{-1}$, we say $A$ is similar to $B$, decomposing $A$ into $PBP^{-1}$ is also called a similarity transformation.
If $n\times n$ matrices $A$ and $B$ are similar, they have the same eigenvalues.
The diagnoalization, which we will explain below, is a process of finding similar matrices.
Let $A$ be an $n\times n$ matrix. If there exists an $n\times n$ invertible matrix $P$ and a diagonal matrix $D$, such that
$$ A=PDP^{-1} $$then matrix $A$ is called a diagonalizable matrix.
And further, the columns of $P$ are linearly independent eigenvectors of $A$, and its corresponding eigenvalues are on the diagonal of $D$. In other words, $A$ is diagonalizable if and only if the dimension of eigenspace basis is $n$.
Let's show why this equation holds.
Let
$$ P = \left[\begin{array}{llll} {v}_{1} & {v}_{2} & \cdots & {v}_{n} \end{array}\right]\\ $$$$ D = \left[\begin{array}{cccc} \lambda_{1} & 0 & \cdots & 0 \\ 0 & \lambda_{2} & \cdots & 0 \\ \vdots & \vdots & & \vdots \\ 0 & 0 & \cdots & \lambda_{n} \end{array}\right] $$where $v_i, i \in (1, 2, ...n)$ is an eigenvector of $A$, $\lambda_i, i \in (1, 2, ...n)$ is an eigenvalue of $A$.
We know that $A{v}_i = \lambda_i{v}_i$, i.e.
$$ AP = PD $$Since $P$ has all independent eigenvectors, then
$$ A = PDP^{-1} $$Strictly speaking, if $A$ is symmetric, i.e. $A=A^T$, what we have just shown is called Spectral decomposition, the similar matrix $D$ holds all the eigenvalues on its diagonal. And $P$ is orthogonal matrix, which means any of of its two columns are perpendicular. Therefore it could be rewritten as $$ A = PDP^{T} $$
What these plots tells are actually the functions of each decomposed matrix, $P$ and $P^T$ are for rotation because they orthogonal, and $D$ are for scaling because it's diagonal.
import numpy as np
import matplotlib.pyplot as plt
# Define matrix A
A = np.array([[1, 3], [3, -5]])
# Calculate eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(A)
# Plot eigenvectors
fig, ax = plt.subplots(figsize=(15, 4), nrows=1, ncols=4)
for i in range(2):
ax[0].quiver(0, 0, eigenvectors[:,i][0], eigenvectors[:,i][1], angles='xy',
scale_units='xy', scale=1, color=['r','b'][i])
ax[0].set_xlim(-2, 2)
ax[0].set_ylim(-2, 2)
for i in range(4):
ax[i].set_xlabel('x')
ax[i].set_ylabel('y')
ax[i].set_aspect('equal')
ax[1].quiver(0, 0, 1, 0, angles='xy',
scale_units='xy', scale=1, color=['r'])
ax[1].quiver(0, 0, 0, 1, angles='xy',
scale_units='xy', scale=1, color=['b'])
ax[1].set_xlim(-2, 2)
ax[1].set_ylim(-2, 2)
ax[1].set_title('$P^T$')
ax[2].quiver(0, 0, eigenvalues[0], 0, angles='xy',
scale_units='xy', scale=1, color=['r'])
ax[2].quiver(0, 0, 0, eigenvalues[1], angles='xy',
scale_units='xy', scale=1, color=['b'])
ax[2].set_xlim(-7, 7)
ax[2].set_ylim(-7, 7)
ax[2].set_title('$DP^T$')
temp = np.array([[eigenvalues[0], 0],
[0, eigenvalues[1]]])
temp1 = temp@eigenvectors
ax[3].quiver(0, 0, temp1[:,0][0], temp1[:,0][1], angles='xy',
scale_units='xy', scale=1, color=['r'])
ax[3].quiver(0, 0, temp1[:,1][0], temp1[:,1][1], angles='xy',
scale_units='xy', scale=1, color=['b'])
ax[3].set_xlim(-7, 7)
ax[3].set_ylim(-7, 7)
ax[3].set_title('$PDP^T$')
plt.show()
However, in reality we don't have many chances of working on symmetric matrices, so spectral decomposition is mostly for theoretical demonstration, not so much application in reality.
Consider a matrix
$$A=\left[\begin{array}{rrr} 1 & 3 & 3 \\ -3 & -5 & -3 \\ 3 & 3 & 1 \end{array}\right]$$We seek to diagonalize the matrix $A$.
Following these steps:
A = sy.Matrix([[1,3,3], [-3, -5, -3], [3,3,1]])
eig = sy.matrices.matrices.MatrixEigen.eigenvects(A)
eig
Reminder the return value takes the form [(eigenval, multiplicity, eigenspace), ...]
.
Construct $P$
P = sy.zeros(3, 3)
P[:, 0] = eig[0][2][0]
P[:, 1] = eig[0][2][1]
P[:, 2] = eig[1][2][0]
P
Construct $D$
D = sy.diag(eig[0][0], eig[0][0], eig[1][0])
D
We can verify if $PDP^{-1}=A$ holds:
P * D * P.inv() == A
True
Of course we don't need to go through this process seperately. There is diagonalize
method in SymPy.
P, D = A.diagonalize()
P
D
We obtain the same results as previous separate steps.
Sometimes we just want to test if a matrix is diagonalizable, then use is_diagonalizable
in SymPy.
A.is_diagonalizable()
True
If $A$ is symmetric, all of its eigenvectors are orthogonal.