Exercises#

Computing powers of a matrix#

Let’s consider the matrix $A=\left[\begin{array}{cc}7& 2\\ -4& 1\end{array}\right]$. We wish to compute $A^k$ for an arbitrary $k\in \mathbb{N}$, given that $A=PA{P}^{-1}$, where $P=\left[\begin{array}{cc}1& 1\\ -1& -2\end{array}\right]$, and $D=\left[\begin{array}{cc}5& 0\\ 0& 3\end{array}\right]$.

Using this decomposition, we have:

Since $D$ is diagonal, $D^k=\left[\begin{array}{cc}{5}^k& 0\\ 0& 3^k\end{array}\right]$. Additionally, $P^{-1}=\left[\begin{array}{cc}2& 1\\ -1& -1\end{array}\right]$. Therefore,

This decomposition, known as diagonalization, enables us to compute $A^k$ easily for large values of $k$. Without such factorization, computing $A^k$ for large matrices and large $k$ can be a time-consuming task. Therefore, it is beneficial to find such decomposition before calculating the power of $A$.

Unfortunately, this decomposition doesn’t exist for all matrices. The goal of this section is to characterize matrices that admit such a factorization.

Diagonalization#

An $n\times n$ matrix $A$ is called diagonalizable if it is similar to a diagonal matrix. More precisely, there exists an invertible $n\times n$ matrix $P$ and an $n\times n$ diagonal matrix $D$ such that $A=PD{P}^{-1}$.

It turns out an $n\times n$ matrix is diagonalizable if we can form a basis ${\mathbb{R}}^n$ using eigenvectors of $A$:

In fact, suppose ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n$ are eigenvalues of $A$, and $\{\overrightarrow{v_1},{\overrightarrow{v}}_2,\dots ,{\overrightarrow{v}}_n\}$ is the linearly independent set of eigenvectors of $A$. If we form the matrix $P$ by taking the eigenvectors as columns, and $D$ is a diagonal matrix with eigenvalues ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n$ on the diagonal, then we have $A=PD{P}^{-1}$.

Finding eigenvectors in general is not an easy task, and is not clear how to check which $n\times n$ matrices have $n$ linearly independent eigenvectors. The next theorem partially addresses this issue.

Combining Theorem 18 and Theorem 19, we obtain the following corollary:

Solution

To determine if $A$ is diagonalizable we need to find eigenvalues and to check if it admits $3$ linearly independent eigenvectors.

• Step 1

we use linalg.eigvas to find the eigenvalues and eigenvectors of $A$:

A = np.array([[1,3,3],[-3,-5,3],[3,3,1]])
# find evalues and evectors
evalues, evectors = linalg.eig(A)

print("eigenvalues = ",  evalues)

eigenvalues =  [-5. -2.  4.]

$A$ has three distinct eigenvalues and therefore it is diagonalizable by Corollary 1.

• Step 2: Constructing the diagonal matrix $D$

lambda_1 = evalues[0]
lambda_2 = evalues[1]
lambda_3 = evalues[2]

D = np.array([[lambda_1, 0, 0], [0, lambda_2, 0], [0, 0, lambda_3]]) 
D

array([[-5.,  0.,  0.],
       [ 0., -2.,  0.],
       [ 0.,  0.,  4.]])

• Step 3 : Finding the matrix $P$ whose $i$th column is an eigenvector corresponding to ${\lambda }_i$ for $i\in \{1,2,3\}$:

# set up the invirtible matrix P
P = evectors

• Step 4: compute the inverse of matrix P

# compute the inverse of P

Q = linalg.inv(P)

• Step 5 (sanity check): Now that we have $D$, $P$, $Q=P^{-1}$, lets check if $A=PDQ$

P @ D @ Q

array([[ 1.,  3.,  3.],
       [-3., -5.,  3.],
       [ 3.,  3.,  1.]])

Solution:

• Step 1: Finding eigenvalues and eigenvectors of $A$:

A = np.array([[1,3,3],[-3,-5,-3],[3,3,1]])

# find evalues and evectors
evalues, evectors = linalg.eig(A)

print("eigenvalues = ",  evalues)

eigenvalues =  [ 1. -2. -2.]

$A$ has 2 eigenvalues: ${\lambda }_1=1$ with multiplicity 1 and ${\lambda }_2=-2$ with multiplicity 2. Note that we cannot use Corollary 1 anymore because eigenvalues are not distinct. To check if $A$ is diagonalizable, we need to find three linearly independent eigenvectors of $A$. Lets have a look at our eigenvectors:

# A matrix whose columns are eigenvectors:
P = evectors

# computes the rank of P
np.linalg.matrix_rank(P)

3

Since $rank(P)=3$ the columns of $P$ are linearly independents. Thus, $A$ is diagonalizable by Theorem 1.

• Step 2: Finding the diagonal matrix $D$

lambda_1 = evalues[0]
lambda_2 = evalues[1]
lambda_3 = evalues[2]

D = np.array([[lambda_1, 0, 0], [0, lambda_2, 0], [0, 0, lambda_3]]) 
D

array([[ 1.,  0.,  0.],
       [ 0., -2.,  0.],
       [ 0.,  0., -2.]])

• Step 3: Finding the matrix $P$ whose $i$th column is an eigenvector corresponding to ${\lambda }_i$ for $i\in \{1,2,3\}$. We already have this matrix so we go to the next step.

• Step 4: Computing the inverse of $P$

# computing the inverse of P

Q = np.linalg.inv(P)
Q

array([[-1.73205081, -1.73205081, -1.73205081],
       [ 0.89919543,  2.75763333,  1.8584379 ],
       [-0.6541932 ,  0.5254512 ,  1.1796444 ]])

• Step 5 (sanity check): lets check if $A=PDQ$

# check A = PDQ 
P @ D @ Q

array([[ 1.,  3.,  3.],
       [-3., -5., -3.],
       [ 3.,  3.,  1.]])

The next theorem provides another way to characterize the diagonalizability of matrices based on the multiplicity of eigenvalues. Before stating the theorem, we need the following definition:

Suppose $A$ is a square matrix and $\lambda$ is an eigenvalue of $A$. The dimension of the corresponding eigenspace $E_{\lambda }=null(A-\lambda I)$ is called the geometric multiplicity of $\lambda$.

Recall that dim($E_{{\lambda }_i}$) = dim null$(A-{\lambda }_{i}I)$ = the number of non-pivot columns in $(A-{\lambda }_{i}I)$

Solution:

$M$ is an upper triangular matrix, and its eigenvalues are on the diagonal: ${\lambda }_1=1$ with algebraic multiplicity 2 and ${\lambda }_2=2$ with algebraic multiplicity 1.

Since we don’t have distinct eigenvalues, Corollary 1 doesn’t apply here. Let’s take a look at our eigenvectors:

import numpy as np
from numpy import linalg


M = np.array([[1, 1, 0], [0, 1, 1], [0, 0, 2]])

evalues, evectors  = np.linalg.eig(M)
p = evectors

np.linalg.matrix_rank(p)

2

Observe that the rank of matrix $P$ is two, which means we don’t have three linearly independent eigenvectors in $M$. However, using Theorem 3, we can confirm that this is not the case, and consequently, $M$ is not diagonalizable. In fact, it is straightforward to check that the geometric multiplicity of $\lambda =1$ is 1 ($\text{dim}(E_1)=1$), while the algebraic multiplicity is 2. Thus, $M$ is not diagonalizable.

Diagonalization as a change of basis#

Recall that similar matrices represent the same linear map under two different choices of a pair of bases for ${\mathbb{R}}^n$ and ${\mathbb{R}}^m$. More precisely, suppose $A=PC{P}^{-1}$ ($C$ is not necessarily diagonal), then the matrix representation of $T$ with respect to the basis $B$ (formed by the columns of $P$) is $C$. Conversely, changing the basis of ${\mathbb{R}}^n$ leads to a matrix representation of $T$ that is similar to $A$ (section 3.4). This particularly applies to diagonalizable matrices because a matrix is diagonalizable if it is similar to a diagonal matrix.

Solution:

From Example 1, we know that $A=PA{P}^{-1}$ where $P=\left[\begin{array}{cc}1& 1\\ -1& -2\end{array}\right]$ and $D=\left[\begin{array}{cc}5& 0\\ 0& 3\end{array}\right]$.

The columns of $P$ are eigenvectors of $A$, and they form a basis $B$ for ${\mathbb{R}}^2$. By Theorem 2, $D$ is the matrix representation of $T$ with respect to $B$. In fact, the mappings $\overrightarrow{x}\to A\overrightarrow{x}$ and $\overrightarrow{x}\to D\overrightarrow{x}$ describe the same linear transformation, relative to different bases.

Exercises#

1. Suppose

a. Is M diagonalizable? If yes find a diagonalization for it.

b. Use this diagonalization to compute $M^{100}$

2. Suppose

Is $C$ diagonalizable?

3. Let $A$ be a $5\times 5$ matrix, with two eigenvalues. One eigenspace is three-dimensional, and the other is two-dimensional. Is $A$ diagonalizable?

4. Compute $A^{200}$ where

Modeling Dynamical Systems as sequence of linear systems#

In many fields (e.g., ecology, biology, economics, etc.), we wish to mathematically model and study a dynamic system that changes over time. We are usually given several features of the system that are measured at (discrete) different time intervals which can be viewed as a sequence of vectors

We interpret $\overrightarrow{x_i}$ as the state of a system at time $t=i$. We also assume that there is a matrix $A$ such that

Our goal is to develop a formula for $x_n$ and describe what can happen when $n$ approaches to infinity.

Let’s consider the movement of populations or groups of people from one region to another. We want a simple model that considers the changes in the population of a city and its surrounding suburbs over a period of years. Let’s say the initial year is 2020 and denote the population of the city and suburbs by $c_0$ and $s_0$, respectively. Now we can form the initial vector population $\overrightarrow{x_0}=\left[\begin{array}{c}{c}_0\\ r_0\end{array}\right]$. Similarly, for 2021 and subsequent years, denote the populations of the city and suburbs by the vectors

How can be these vectors mathematically related?

Suppose demographic studies show that each year about $5\mathrm{\%}$ of the city’s population moves to the suburbs (and $95\mathrm{\%}$ remains in the city), while $3\mathrm{\%}$ of the suburban population moves to the city (and $97\mathrm{\%}$ remains in the suburbs). With this information, we can find the vector population $\overrightarrow{x_1}$ based on the initial population $\overrightarrow{x_0}$. After the first year, $c_{0}95\mathrm{\%}$ remain in the city and $c_{0}0.05$ move to the suburbs. Additionally, $s_0.97\mathrm{\%}$ remain in the suburbs and $s_{0}0.03$ move to the city. Therefore,

The matrix $M=\left[\begin{array}{cc}0.95& 0.03\\ 0.97& 0.05\end{array}\right]$ is called the migration matrix. Consequently, we can write

and the sequence

describes the population of the city/suburbs over a period of years.

Solution:

import numpy as np

# migration matrix
M = np.array([[0.95, 0.03] ,[0.05, 0.097]])


# initial population vector x_0
x_0 = np.array([600000, 40000])

# population vector x_1
x_1 = M @ x_0

print('the population of city in 2021 was = ', x_1[0],'\n')
print('the population of suburbs in 2021 was = ', x_1[1], '\n\n')

# population vector x_2

x_2 = M @ x_1

print('the population of the city in 2022 was = ', x_2[0],'\n')
print('the population of suburbs in 2022 was = ', x_1[1], '\n\n')

the population of city in 2021 was =  571200.0 

the population of suburbs in 2021 was =  33880.0 


the population of the city in 2022 was =  543656.4 

the population of suburbs in 2022 was =  33880.0 

Long-term behavior of dynamical systems#

Eigenvalues and eigenvectors can be used to understand the long-term behavior or evolution of a dynamical system described by an equation in the form of ${\overrightarrow{x}}_{k+1}=A{\overrightarrow{x}}_k$. We assume that $A$ is diagonalizable with $n$ linearly independent eigenvectors $\{\overrightarrow{v_1},{\overrightarrow{v}}_2,\dots ,{\overrightarrow{v}}_n\}$ and corresponding eigenvalues ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n$. For convenience, let’s assume the eigenvectors are arranged such that $|{\lambda }_1|\ge |{\lambda }_2|\ge \cdots \ge |{\lambda }_n|$. Since $\{\overrightarrow{v_1},{\overrightarrow{v}}_2,\dots ,{\overrightarrow{v}}_n\}$ forms a basis for ${\mathbb{R}}^n$, any initial vector ${\overrightarrow{x}}_0$ can be expressed uniquely as:

By applying $A$ to both sides, we have:

And in general:

The above equation represents the evolution of the system over time. It shows that each component of the vector ${\overrightarrow{x}}_k$ is a linear combination of the eigenvectors $\overrightarrow{v_i}$, scaled by the corresponding eigenvalues ${\lambda }_i$ raised to the power of $k$.

As an example, let’s consider a predator-prey system involving owls and wood rats. The population vectors ${\overrightarrow{x}}_k$ at time $k$ are given by:

where $O_k$ is the population of owls and $R_k$ is the population of rats (measured in thousands). The populations evolve according to the following equations:

Here, the coefficient $0.5$ in the first equation says that, without wood rats for food, only half of the owls will survive each month. The coefficient $1.1$ in the second equation says that, in the absence of owls as predators, the rat population will grow by 10% each month. The parameter $p$ is a positive value to be specified.

To determine the evolution of this system, we can compute the eigenvalues and eigenvectors of the coefficient matrix $A$, and then apply equation $(\ast \ast )$. For $p=0.104$, lets find the evolution of this system:

# coefficient matrix
A = np.array([[0.5, 0.4],[-0.104, 1.1]])

# Computing the eigenvalues and eigenvectors of A
np.linalg.eigvals(A)

array([0.58, 1.02])

Note that $|{\lambda }_1|\le 1$ and $|{\lambda }_2|\ge 1$. A straightforward calculation shows that the corresponding eigenvectors are $\overrightarrow{v_1}=\left[\begin{array}{c}5\\ 1\end{array}\right]$ and $\overrightarrow{v_2}=\left[\begin{array}{c}10\\ 13\end{array}\right]$, respectively. The initial vector $\overrightarrow{x_0}$ can be expressed as a linear combination of the eigenvectors:

Therefore, according to equation $(\ast \ast )$, we have:

As $k\to \mathrm{\infty }$, the term $({\lambda }_1{)}^k\overrightarrow{v_1}$ approaches zero, allowing us to approximate:

Furthermore, observe that:

The above approximation shows that both the number of owls and rats grow by an approximate factor of 1.02 each month.

let’s calculate some population vectors, assuming the initial population vector is $\left[\begin{array}{c}10\\ 20\end{array}\right]$.

import numpy as np

# A matrix whose columns are eigenvectors
evectors = np.array([[5, 10], [1, 13]])

# Initial vector
x_0 = np.array([1, 2])

# Solving the equation for d_1 and d_2
d = np.linalg.solve(evectors, x_0)

# Extracting d_2
d_2 = d[1]

# Eigenvalues corresponding to the eigenvectors
evalues = np.array([10, 20])  # Replace with the actual eigenvalues

# Iterating over the range
for k in range(10):
    # Calculating the population using the eigenvalues, eigenvectors, and d_2
    population = d_2 * pow(evalues[1], k+1) * evectors[1]
    print('At month', k+1, 'the owl population is', population[0],' and ', 'the rat population is', population[1],'\n')

At month 1 the owl population is 3.2727272727272725  and  the rat population is 42.54545454545454 

At month 2 the owl population is 65.45454545454545  and  the rat population is 850.9090909090909 

At month 3 the owl population is 1309.090909090909  and  the rat population is 17018.181818181816 

At month 4 the owl population is 26181.81818181818  and  the rat population is 340363.63636363635 

At month 5 the owl population is 523636.36363636365  and  the rat population is 6807272.7272727275 

At month 6 the owl population is 10472727.272727273  and  the rat population is 136145454.54545456 

At month 7 the owl population is 209454545.45454547  and  the rat population is 2722909090.909091 

At month 8 the owl population is -27786072.436363637  and  the rat population is -361218941.6727273 

At month 9 the owl population is 147091381.52727273  and  the rat population is 1912187959.8545456 

At month 10 the owl population is 130576309.52727273  and  the rat population is 1697492023.8545456 

Note that if the absolute values of both eigenvalues are less than 1, then both populations will tend towards zero for large values of $k$. If ${\lambda }_1=1$ and $\lambda 2<1$, then

Finally, if ${\lambda }_1=1$ and $\lambda 2\ge 1$, then

Exercise#