ℹ️ Skipped - page is already crawled
| Filter | Status | Condition | Details |
|---|---|---|---|
| HTTP status | PASS | download_http_code = 200 | HTTP 200 |
| Age cutoff | PASS | download_stamp > now() - 6 MONTH | 0 months ago |
| History drop | PASS | isNull(history_drop_reason) | No drop reason |
| Spam/ban | PASS | fh_dont_index != 1 AND ml_spam_score = 0 | ml_spam_score=0 |
| Canonical | PASS | meta_canonical IS NULL OR = '' OR = src_unparsed | Not set |
| Property | Value |
|---|---|
| URL | https://intro.quantecon.org/eigen_I.html |
| Last Crawled | 2026-04-14 17:13:53 (6 hours ago) |
| First Indexed | 2023-06-15 13:03:43 (2 years ago) |
| HTTP Status Code | 200 |
| Meta Title | 17. Eigenvalues and Eigenvectors — A First Course in Quantitative Economics with Python |
| Meta Description | This |
| Meta Canonical | null |
| Boilerpipe Text | 17.1.
Overview
#
Eigenvalues and eigenvectors are a relatively advanced topic in linear algebra.
At the same time, these concepts are extremely useful for
economic modeling (especially dynamics!)
statistics
some parts of applied mathematics
machine learning
and many other fields of science.
In this lecture we explain the basics of eigenvalues and eigenvectors and introduce the Neumann Series Lemma.
We assume in this lecture that students are familiar with matrices
and understand
the basics of matrix algebra
.
We will use the following imports:
import
matplotlib.pyplot
as
plt
import
numpy
as
np
from
numpy.linalg
import
matrix_power
from
matplotlib.lines
import
Line2D
from
matplotlib.patches
import
FancyArrowPatch
from
mpl_toolkits.mplot3d
import
proj3d
17.2.
Matrices as transformations
#
Let’s start by discussing an important concept concerning matrices.
17.2.1.
Mapping vectors to vectors
#
One way to think about a matrix is as a rectangular collection of
numbers.
Another way to think about a matrix is as a
map
(i.e., as a function) that
transforms vectors to new vectors.
To understand the second point of view, suppose we multiply an
n
×
m
matrix
A
with an
m
×
1
column vector
x
to obtain an
n
×
1
column vector
y
:
A
x
=
y
If we fix
A
and consider different choices of
x
, we can understand
A
as
a map transforming
x
to
A
x
.
Because
A
is
n
×
m
, it transforms
m
-vectors to
n
-vectors.
We can write this formally as
A
:
R
m
→
R
n
.
You might argue that if
A
is a function then we should write
A
(
x
)
=
y
rather than
A
x
=
y
but the second notation is more conventional.
17.2.2.
Square matrices
#
Let’s restrict our discussion to square matrices.
In the above discussion, this means that
m
=
n
and
A
maps
R
n
to
itself.
This means
A
is an
n
×
n
matrix that maps (or “transforms”) a vector
x
in
R
n
to a new vector
y
=
A
x
also in
R
n
.
Example 17.1
[
2
1
−
1
1
]
[
1
3
]
=
[
5
2
]
Here, the matrix
A
=
[
2
1
−
1
1
]
transforms the vector
x
=
[
1
3
]
to the vector
y
=
[
5
2
]
.
Let’s visualize this using Python:
A
=
np
.
array
([[
2
,
1
],
[
-
1
,
1
]])
from
math
import
sqrt
fig
,
ax
=
plt
.
subplots
()
# Set the axes through the origin
for
spine
in
[
'left'
,
'bottom'
]:
ax
.
spines
[
spine
]
.
set_position
(
'zero'
)
for
spine
in
[
'right'
,
'top'
]:
ax
.
spines
[
spine
]
.
set_color
(
'none'
)
ax
.
set
(
xlim
=
(
-
2
,
6
),
ylim
=
(
-
2
,
4
),
aspect
=
1
)
vecs
=
((
1
,
3
),
(
5
,
2
))
c
=
[
'r'
,
'black'
]
for
i
,
v
in
enumerate
(
vecs
):
ax
.
annotate
(
''
,
xy
=
v
,
xytext
=
(
0
,
0
),
arrowprops
=
dict
(
color
=
c
[
i
],
shrink
=
0
,
alpha
=
0.7
,
width
=
0.5
))
ax
.
text
(
0.2
+
1
,
0.2
+
3
,
'x=$(1,3)$'
)
ax
.
text
(
0.2
+
5
,
0.2
+
2
,
'Ax=$(5,2)$'
)
ax
.
annotate
(
''
,
xy
=
(
sqrt
(
10
/
29
)
*
5
,
sqrt
(
10
/
29
)
*
2
),
xytext
=
(
0
,
0
),
arrowprops
=
dict
(
color
=
'purple'
,
shrink
=
0
,
alpha
=
0.7
,
width
=
0.5
))
ax
.
annotate
(
''
,
xy
=
(
1
,
2
/
5
),
xytext
=
(
1
/
3
,
1
),
arrowprops
=
{
'arrowstyle'
:
'->'
,
'connectionstyle'
:
'arc3,rad=-0.3'
},
horizontalalignment
=
'center'
)
ax
.
text
(
0.8
,
0.8
,
f
'θ'
,
fontsize
=
14
)
plt
.
show
()
One way to understand this transformation is that
A
first rotates
x
by some angle
θ
and
then scales it by some scalar
γ
to obtain the image
y
of
x
.
17.3.
Types of transformations
#
Let’s examine some standard transformations we can perform with matrices.
Below we visualize transformations by thinking of vectors as points
instead of arrows.
We consider how a given matrix transforms
a grid of points and
a set of points located on the unit circle in
R
2
.
To build the transformations we will use two functions, called
grid_transform
and
circle_transform
.
Each of these functions visualizes the actions of a given
2
×
2
matrix
A
.
Show source
Hide code cell source
def
colorizer
(
x
,
y
):
r
=
min
(
1
,
1
-
y
/
3
)
g
=
min
(
1
,
1
+
y
/
3
)
b
=
1
/
4
+
x
/
16
return
(
r
,
g
,
b
)
def
grid_transform
(
A
=
np
.
array
([[
1
,
-
1
],
[
1
,
1
]])):
xvals
=
np
.
linspace
(
-
4
,
4
,
9
)
yvals
=
np
.
linspace
(
-
3
,
3
,
7
)
xygrid
=
np
.
column_stack
([[
x
,
y
]
for
x
in
xvals
for
y
in
yvals
])
uvgrid
=
A
@
xygrid
colors
=
list
(
map
(
colorizer
,
xygrid
[
0
],
xygrid
[
1
]))
fig
,
ax
=
plt
.
subplots
(
1
,
2
,
figsize
=
(
10
,
5
))
for
axes
in
ax
:
axes
.
set
(
xlim
=
(
-
11
,
11
),
ylim
=
(
-
11
,
11
))
axes
.
set_xticks
([])
axes
.
set_yticks
([])
for
spine
in
[
'left'
,
'bottom'
]:
axes
.
spines
[
spine
]
.
set_position
(
'zero'
)
for
spine
in
[
'right'
,
'top'
]:
axes
.
spines
[
spine
]
.
set_color
(
'none'
)
# Plot x-y grid points
ax
[
0
]
.
scatter
(
xygrid
[
0
],
xygrid
[
1
],
s
=
36
,
c
=
colors
,
edgecolor
=
"none"
)
# ax[0].grid(True)
# ax[0].axis("equal")
ax
[
0
]
.
set_title
(
"points $x_1, x_2, \cdots, x_k$"
)
# Plot transformed grid points
ax
[
1
]
.
scatter
(
uvgrid
[
0
],
uvgrid
[
1
],
s
=
36
,
c
=
colors
,
edgecolor
=
"none"
)
# ax[1].grid(True)
# ax[1].axis("equal")
ax
[
1
]
.
set_title
(
"points $Ax_1, Ax_2, \cdots, Ax_k$"
)
plt
.
show
()
def
circle_transform
(
A
=
np
.
array
([[
-
1
,
2
],
[
0
,
1
]])):
fig
,
ax
=
plt
.
subplots
(
1
,
2
,
figsize
=
(
10
,
5
))
for
axes
in
ax
:
axes
.
set
(
xlim
=
(
-
4
,
4
),
ylim
=
(
-
4
,
4
))
axes
.
set_xticks
([])
axes
.
set_yticks
([])
for
spine
in
[
'left'
,
'bottom'
]:
axes
.
spines
[
spine
]
.
set_position
(
'zero'
)
for
spine
in
[
'right'
,
'top'
]:
axes
.
spines
[
spine
]
.
set_color
(
'none'
)
θ
=
np
.
linspace
(
0
,
2
*
np
.
pi
,
150
)
r
=
1
θ_1
=
np
.
empty
(
12
)
for
i
in
range
(
12
):
θ_1
[
i
]
=
2
*
np
.
pi
*
(
i
/
12
)
x
=
r
*
np
.
cos
(
θ
)
y
=
r
*
np
.
sin
(
θ
)
a
=
r
*
np
.
cos
(
θ_1
)
b
=
r
*
np
.
sin
(
θ_1
)
a_1
=
a
.
reshape
(
1
,
-
1
)
b_1
=
b
.
reshape
(
1
,
-
1
)
colors
=
list
(
map
(
colorizer
,
a
,
b
))
ax
[
0
]
.
plot
(
x
,
y
,
color
=
'black'
,
zorder
=
1
)
ax
[
0
]
.
scatter
(
a_1
,
b_1
,
c
=
colors
,
alpha
=
1
,
s
=
60
,
edgecolors
=
'black'
,
zorder
=
2
)
ax
[
0
]
.
set_title
(
r
"unit circle in $\mathbb
{R}
^2$"
)
x1
=
x
.
reshape
(
1
,
-
1
)
y1
=
y
.
reshape
(
1
,
-
1
)
ab
=
np
.
concatenate
((
a_1
,
b_1
),
axis
=
0
)
transformed_ab
=
A
@
ab
transformed_circle_input
=
np
.
concatenate
((
x1
,
y1
),
axis
=
0
)
transformed_circle
=
A
@
transformed_circle_input
ax
[
1
]
.
plot
(
transformed_circle
[
0
,
:],
transformed_circle
[
1
,
:],
color
=
'black'
,
zorder
=
1
)
ax
[
1
]
.
scatter
(
transformed_ab
[
0
,
:],
transformed_ab
[
1
:,],
color
=
colors
,
alpha
=
1
,
s
=
60
,
edgecolors
=
'black'
,
zorder
=
2
)
ax
[
1
]
.
set_title
(
"transformed circle"
)
plt
.
show
()
17.3.1.
Scaling
#
A matrix of the form
[
α
0
0
β
]
scales vectors across the x-axis by a factor
α
and along the y-axis by
a factor
β
.
Here we illustrate a simple example where
α
=
β
=
3
.
A
=
np
.
array
([[
3
,
0
],
# scaling by 3 in both directions
[
0
,
3
]])
grid_transform
(
A
)
circle_transform
(
A
)
17.3.2.
Shearing
#
A “shear” matrix of the form
[
1
λ
0
1
]
stretches vectors along the x-axis by an amount proportional to the
y-coordinate of a point.
A
=
np
.
array
([[
1
,
2
],
# shear along x-axis
[
0
,
1
]])
grid_transform
(
A
)
circle_transform
(
A
)
17.3.3.
Rotation
#
A matrix of the form
[
cos
θ
sin
θ
−
sin
θ
cos
θ
]
is called a
rotation matrix
.
This matrix rotates vectors clockwise by an angle
θ
.
θ
=
np
.
pi
/
4
# 45 degree clockwise rotation
A
=
np
.
array
([[
np
.
cos
(
θ
),
np
.
sin
(
θ
)],
[
-
np
.
sin
(
θ
),
np
.
cos
(
θ
)]])
grid_transform
(
A
)
17.3.4.
Permutation
#
The permutation matrix
[
0
1
1
0
]
interchanges the coordinates of a vector.
A
=
np
.
column_stack
([[
0
,
1
],
[
1
,
0
]])
grid_transform
(
A
)
More examples of common transition matrices can be found
here
.
17.4.
Matrix multiplication as composition
#
Since matrices act as functions that transform one vector to another, we can
apply the concept of function composition to matrices as well.
17.4.1.
Linear compositions
#
Consider the two matrices
A
=
[
0
1
−
1
0
]
and
B
=
[
1
2
0
1
]
What will the output be when we try to obtain
A
B
x
for some
2
×
1
vector
x
?
[
0
1
−
1
0
]
⏟
A
[
1
2
0
1
]
⏟
B
[
1
3
]
⏞
x
→
[
0
1
−
1
−
2
]
⏟
A
B
[
1
3
]
⏞
x
→
[
3
−
7
]
⏞
y
[
0
1
−
1
0
]
⏟
A
[
1
2
0
1
]
⏟
B
[
1
3
]
⏞
x
→
[
0
1
−
1
0
]
⏟
A
[
7
3
]
⏞
B
x
→
[
3
−
7
]
⏞
y
We can observe that applying the transformation
A
B
on the vector
x
is the
same as first applying
B
on
x
and then applying
A
on the vector
B
x
.
Thus the matrix product
A
B
is the
composition
of the
matrix transformations
A
and
B
This means first apply transformation
B
and then
transformation
A
.
When we matrix multiply an
n
×
m
matrix
A
with an
m
×
k
matrix
B
the obtained matrix product is an
n
×
k
matrix
A
B
.
Thus, if
A
and
B
are transformations such that
A
:
R
m
→
R
n
and
B
:
R
k
→
R
m
, then
A
B
transforms
R
k
to
R
n
.
Viewing matrix multiplication as composition of maps helps us
understand why, under matrix multiplication,
A
B
is generally not equal to
B
A
.
(After all, when we compose functions, the order usually matters.)
17.4.2.
Examples
#
Let
A
be the
90
∘
clockwise rotation matrix given by
[
0
1
−
1
0
]
and let
B
be a shear matrix
along the x-axis given by
[
1
2
0
1
]
.
We will visualize how a grid of points changes when we apply the
transformation
A
B
and then compare it with the transformation
B
A
.
Show source
Hide code cell source
def
grid_composition_transform
(
A
=
np
.
array
([[
1
,
-
1
],
[
1
,
1
]]),
B
=
np
.
array
([[
1
,
-
1
],
[
1
,
1
]])):
xvals
=
np
.
linspace
(
-
4
,
4
,
9
)
yvals
=
np
.
linspace
(
-
3
,
3
,
7
)
xygrid
=
np
.
column_stack
([[
x
,
y
]
for
x
in
xvals
for
y
in
yvals
])
uvgrid
=
B
@
xygrid
abgrid
=
A
@
uvgrid
colors
=
list
(
map
(
colorizer
,
xygrid
[
0
],
xygrid
[
1
]))
fig
,
ax
=
plt
.
subplots
(
1
,
3
,
figsize
=
(
15
,
5
))
for
axes
in
ax
:
axes
.
set
(
xlim
=
(
-
12
,
12
),
ylim
=
(
-
12
,
12
))
axes
.
set_xticks
([])
axes
.
set_yticks
([])
for
spine
in
[
'left'
,
'bottom'
]:
axes
.
spines
[
spine
]
.
set_position
(
'zero'
)
for
spine
in
[
'right'
,
'top'
]:
axes
.
spines
[
spine
]
.
set_color
(
'none'
)
# Plot grid points
ax
[
0
]
.
scatter
(
xygrid
[
0
],
xygrid
[
1
],
s
=
36
,
c
=
colors
,
edgecolor
=
"none"
)
ax
[
0
]
.
set_title
(
r
"points $x_1, x_2, \cdots, x_k$"
)
# Plot intermediate grid points
ax
[
1
]
.
scatter
(
uvgrid
[
0
],
uvgrid
[
1
],
s
=
36
,
c
=
colors
,
edgecolor
=
"none"
)
ax
[
1
]
.
set_title
(
r
"points $Bx_1, Bx_2, \cdots, Bx_k$"
)
# Plot transformed grid points
ax
[
2
]
.
scatter
(
abgrid
[
0
],
abgrid
[
1
],
s
=
36
,
c
=
colors
,
edgecolor
=
"none"
)
ax
[
2
]
.
set_title
(
r
"points $ABx_1, ABx_2, \cdots, ABx_k$"
)
plt
.
show
()
A
=
np
.
array
([[
0
,
1
],
# 90 degree clockwise rotation
[
-
1
,
0
]])
B
=
np
.
array
([[
1
,
2
],
# shear along x-axis
[
0
,
1
]])
17.4.2.1.
Shear then rotate
#
grid_composition_transform
(
A
,
B
)
# transformation AB
17.4.2.2.
Rotate then shear
#
grid_composition_transform
(
B
,
A
)
# transformation BA
It is evident that the transformation
A
B
is not the same as the transformation
B
A
.
17.5.
Iterating on a fixed map
#
In economics (and especially in dynamic modeling), we are often interested in
analyzing behavior where we repeatedly apply a fixed matrix.
For example, given a vector
v
and a matrix
A
, we are interested in
studying the sequence
v
,
A
v
,
A
A
v
=
A
2
v
,
…
Let’s first see examples of a sequence of iterates
(
A
k
v
)
k
≥
0
under
different maps
A
.
def
plot_series
(
A
,
v
,
n
):
B
=
np
.
array
([[
1
,
-
1
],
[
1
,
0
]])
fig
,
ax
=
plt
.
subplots
()
ax
.
set
(
xlim
=
(
-
4
,
4
),
ylim
=
(
-
4
,
4
))
ax
.
set_xticks
([])
ax
.
set_yticks
([])
for
spine
in
[
'left'
,
'bottom'
]:
ax
.
spines
[
spine
]
.
set_position
(
'zero'
)
for
spine
in
[
'right'
,
'top'
]:
ax
.
spines
[
spine
]
.
set_color
(
'none'
)
θ
=
np
.
linspace
(
0
,
2
*
np
.
pi
,
150
)
r
=
2.5
x
=
r
*
np
.
cos
(
θ
)
y
=
r
*
np
.
sin
(
θ
)
x1
=
x
.
reshape
(
1
,
-
1
)
y1
=
y
.
reshape
(
1
,
-
1
)
xy
=
np
.
concatenate
((
x1
,
y1
),
axis
=
0
)
ellipse
=
B
@
xy
ax
.
plot
(
ellipse
[
0
,
:],
ellipse
[
1
,
:],
color
=
'black'
,
linestyle
=
(
0
,
(
5
,
10
)),
linewidth
=
0.5
)
# Initialize holder for trajectories
colors
=
plt
.
cm
.
rainbow
(
np
.
linspace
(
0
,
1
,
20
))
for
i
in
range
(
n
):
iteration
=
matrix_power
(
A
,
i
)
@
v
v1
=
iteration
[
0
]
v2
=
iteration
[
1
]
ax
.
scatter
(
v1
,
v2
,
color
=
colors
[
i
])
if
i
==
0
:
ax
.
text
(
v1
+
0.25
,
v2
,
f
'$v$'
)
elif
i
==
1
:
ax
.
text
(
v1
+
0.25
,
v2
,
f
'$Av$'
)
elif
1
<
i
<
4
:
ax
.
text
(
v1
+
0.25
,
v2
,
f
'$A^
{
i
}
v$'
)
plt
.
show
()
A
=
np
.
array
([[
sqrt
(
3
)
+
1
,
-
2
],
[
1
,
sqrt
(
3
)
-
1
]])
A
=
(
1
/
(
2
*
sqrt
(
2
)))
*
A
v
=
(
-
3
,
-
3
)
n
=
12
plot_series
(
A
,
v
,
n
)
Here with each iteration the vectors get shorter, i.e., move closer to the origin.
In this case, repeatedly multiplying a vector by
A
makes the vector “spiral in”.
B
=
np
.
array
([[
sqrt
(
3
)
+
1
,
-
2
],
[
1
,
sqrt
(
3
)
-
1
]])
B
=
(
1
/
2
)
*
B
v
=
(
2.5
,
0
)
n
=
12
plot_series
(
B
,
v
,
n
)
Here with each iteration vectors do not tend to get longer or shorter.
In this case, repeatedly multiplying a vector by
A
simply “rotates it around
an ellipse”.
B
=
np
.
array
([[
sqrt
(
3
)
+
1
,
-
2
],
[
1
,
sqrt
(
3
)
-
1
]])
B
=
(
1
/
sqrt
(
2
))
*
B
v
=
(
-
1
,
-
0.25
)
n
=
6
plot_series
(
B
,
v
,
n
)
Here with each iteration vectors tend to get longer, i.e., farther from the
origin.
In this case, repeatedly multiplying a vector by
A
makes the vector “spiral out”.
We thus observe that the sequence
(
A
k
v
)
k
≥
0
behaves differently depending on the map
A
itself.
We now discuss the property of A that determines this behavior.
17.6.
Eigenvalues
#
In this section we introduce the notions of eigenvalues and eigenvectors.
17.6.1.
Definitions
#
Let
A
be an
n
×
n
square matrix.
If
λ
is scalar and
v
is a non-zero
n
-vector such that
A
v
=
λ
v
.
Then we say that
λ
is an
eigenvalue
of
A
, and
v
is the corresponding
eigenvector
.
Thus, an eigenvector of
A
is a nonzero vector
v
such that when the map
A
is
applied,
v
is merely scaled.
The next figure shows two eigenvectors (blue arrows) and their images under
A
(red arrows).
As expected, the image
A
v
of each
v
is just a scaled version of the original
17.6.2.
Complex values
#
So far our definition of eigenvalues and eigenvectors seems straightforward.
There is one complication we haven’t mentioned yet:
When solving
A
v
=
λ
v
,
λ
is allowed to be a complex number and
v
is allowed to be an
n
-vector of complex numbers.
We will see some examples below.
17.6.3.
Some mathematical details
#
We note some mathematical details for more advanced readers.
(Other readers can skip to the next section.)
The eigenvalue equation is equivalent to
(
A
−
λ
I
)
v
=
0
.
This equation has a nonzero solution
v
only when the columns of
A
−
λ
I
are linearly dependent.
This in turn is equivalent to stating the determinant is zero.
Hence, to find all eigenvalues, we can look for
λ
such that the
determinant of
A
−
λ
I
is zero.
This problem can be expressed as one of solving for the roots of a polynomial
in
λ
of degree
n
.
This in turn implies the existence of
n
solutions in the complex
plane, although some might be repeated.
17.6.4.
Facts
#
Some nice facts about the eigenvalues of a square matrix
A
are as follows:
the determinant of
A
equals the product of the eigenvalues
the trace of
A
(the sum of the elements on the principal diagonal) equals the sum of the eigenvalues
if
A
is symmetric, then all of its eigenvalues are real
if
A
is invertible and
λ
1
,
…
,
λ
n
are its eigenvalues, then the eigenvalues of
A
−
1
are
1
/
λ
1
,
…
,
1
/
λ
n
.
A corollary of the last statement is that a matrix is invertible if and only if all its eigenvalues are nonzero.
17.6.5.
Computation
#
Using NumPy, we can solve for the eigenvalues and eigenvectors of a matrix as follows
from
numpy.linalg
import
eig
A
=
((
1
,
2
),
(
2
,
1
))
A
=
np
.
array
(
A
)
evals
,
evecs
=
eig
(
A
)
evals
# eigenvalues
array([ 3., -1.])
evecs
# eigenvectors
array([[ 0.70710678, -0.70710678],
[ 0.70710678, 0.70710678]])
Note that the
columns
of
evecs
are the eigenvectors.
Since any scalar multiple of an eigenvector is an eigenvector with the same
eigenvalue (which can be verified), the
eig
routine normalizes the length of each eigenvector
to one.
The eigenvectors and eigenvalues of a map
A
determine how a vector
v
is transformed when we repeatedly multiply by
A
.
This is discussed further later.
17.7.
The Neumann Series Lemma
#
In this section we present a famous result about series of matrices that has
many applications in economics.
17.7.1.
Scalar series
#
Here’s a fundamental result about series:
If
a
is a number and
|
a
|
<
1
, then
(17.1)
#
∑
k
=
0
∞
a
k
=
1
1
−
a
=
(
1
−
a
)
−
1
For a one-dimensional linear equation
x
=
a
x
+
b
where x is unknown we can thus conclude that the solution
x
∗
is given by:
x
∗
=
b
1
−
a
=
∑
k
=
0
∞
a
k
b
17.7.2.
Matrix series
#
A generalization of this idea exists in the matrix setting.
Consider the system of equations
x
=
A
x
+
b
where
A
is an
n
×
n
square matrix and
x
and
b
are both column vectors in
R
n
.
Using matrix algebra we can conclude that the solution to this system of equations will be given by:
What guarantees the existence of a unique vector
x
∗
that satisfies
(17.2)
?
The following is a fundamental result in functional analysis that generalizes
(17.1)
to a multivariate case.
Theorem 17.1
(Neumann Series Lemma)
Let
A
be a square matrix and let
A
k
be the
k
-th power of
A
.
Let
r
(
A
)
be the
spectral radius
of
A
, defined as
max
i
|
λ
i
|
, where
{
λ
i
}
i
is the set of eigenvalues of
A
and
|
λ
i
|
is the modulus of the complex number
λ
i
Neumann’s Theorem states the following: If
r
(
A
)
<
1
, then
I
−
A
is invertible, and
(
I
−
A
)
−
1
=
∑
k
=
0
∞
A
k
We can see the Neumann Series Lemma in action in the following example.
A
=
np
.
array
([[
0.4
,
0.1
],
[
0.7
,
0.2
]])
evals
,
evecs
=
eig
(
A
)
# finding eigenvalues and eigenvectors
r
=
max
(
abs
(
λ
)
for
λ
in
evals
)
# compute spectral radius
print
(
r
)
0.5828427124746189
The spectral radius
r
(
A
)
obtained is less than 1.
Thus, we can apply the Neumann Series Lemma to find
(
I
−
A
)
−
1
.
I
=
np
.
identity
(
2
)
# 2 x 2 identity matrix
B
=
I
-
A
B_inverse
=
np
.
linalg
.
inv
(
B
)
# direct inverse method
A_sum
=
np
.
zeros
((
2
,
2
))
# power series sum of A
A_power
=
I
for
i
in
range
(
50
):
A_sum
+=
A_power
A_power
=
A_power
@
A
Let’s check equality between the sum and the inverse methods.
np
.
allclose
(
A_sum
,
B_inverse
)
True
Although we truncate the infinite sum at
k
=
50
, both methods give us the same
result which illustrates the result of the Neumann Series Lemma.
17.8.
Exercises
#
Exercise 17.1
Power iteration is a method for finding the greatest absolute eigenvalue of a diagonalizable matrix.
The method starts with a random vector
b
0
and repeatedly applies the matrix
A
to it
b
k
+
1
=
A
b
k
‖
A
b
k
‖
A thorough discussion of the method can be found
here
.
In this exercise, first implement the power iteration method and use it to find the greatest absolute eigenvalue and its corresponding eigenvector.
Then visualize the convergence.
Exercise 17.2
We have discussed the trajectory of the vector
v
after being transformed by
A
.
Consider the matrix
A
=
[
1
2
1
1
]
and the vector
v
=
[
2
−
2
]
.
Try to compute the trajectory of
v
after being transformed by
A
for
n
=
4
iterations and plot the result.
Exercise 17.3
Previously
, we demonstrated the trajectory of the vector
v
after being transformed by
A
for three different matrices.
Use the visualization in the previous exercise to explain the trajectory of the vector
v
after being transformed by
A
for the three different matrices. |
| Markdown | On this page
- [17\.1. Overview](https://intro.quantecon.org/eigen_I.html#overview)
- [17\.2. Matrices as transformations](https://intro.quantecon.org/eigen_I.html#matrices-as-transformations)
- [17\.2.1. Mapping vectors to vectors](https://intro.quantecon.org/eigen_I.html#mapping-vectors-to-vectors)
- [17\.2.2. Square matrices](https://intro.quantecon.org/eigen_I.html#square-matrices)
- [17\.3. Types of transformations](https://intro.quantecon.org/eigen_I.html#types-of-transformations)
- [17\.3.1. Scaling](https://intro.quantecon.org/eigen_I.html#scaling)
- [17\.3.2. Shearing](https://intro.quantecon.org/eigen_I.html#shearing)
- [17\.3.3. Rotation](https://intro.quantecon.org/eigen_I.html#rotation)
- [17\.3.4. Permutation](https://intro.quantecon.org/eigen_I.html#permutation)
- [17\.4. Matrix multiplication as composition](https://intro.quantecon.org/eigen_I.html#matrix-multiplication-as-composition)
- [17\.4.1. Linear compositions](https://intro.quantecon.org/eigen_I.html#linear-compositions)
- [17\.4.2. Examples](https://intro.quantecon.org/eigen_I.html#examples)
- [17\.4.2.1. Shear then rotate](https://intro.quantecon.org/eigen_I.html#shear-then-rotate)
- [17\.4.2.2. Rotate then shear](https://intro.quantecon.org/eigen_I.html#rotate-then-shear)
- [17\.5. Iterating on a fixed map](https://intro.quantecon.org/eigen_I.html#iterating-on-a-fixed-map)
- [17\.6. Eigenvalues](https://intro.quantecon.org/eigen_I.html#eigenvalues)
- [17\.6.1. Definitions](https://intro.quantecon.org/eigen_I.html#definitions)
- [17\.6.2. Complex values](https://intro.quantecon.org/eigen_I.html#complex-values)
- [17\.6.3. Some mathematical details](https://intro.quantecon.org/eigen_I.html#some-mathematical-details)
- [17\.6.4. Facts](https://intro.quantecon.org/eigen_I.html#facts)
- [17\.6.5. Computation](https://intro.quantecon.org/eigen_I.html#computation)
- [17\.7. The Neumann Series Lemma](https://intro.quantecon.org/eigen_I.html#the-neumann-series-lemma)
- [17\.7.1. Scalar series](https://intro.quantecon.org/eigen_I.html#scalar-series)
- [17\.7.2. Matrix series](https://intro.quantecon.org/eigen_I.html#matrix-series)
- [17\.8. Exercises](https://intro.quantecon.org/eigen_I.html#exercises)
[](https://quantecon.org/) [](https://quantecon.org/)
Powered by [Jupyter Book](https://jupyterbook.org/v1/)
[**Back to top**](https://intro.quantecon.org/eigen_I.html#top)
[A First Course in Quantitative Economics with Python](https://intro.quantecon.org/intro.html)
[Thomas J. Sargent](http://www.tomsargent.com/) and [John Stachurski](https://johnstachurski.net/)
Last changed: Aug 01, 2025 ▼
#### Changelog ([full history](https://github.com/QuantEcon/lecture-python-intro/commits/main/lectures/eigen_I.md))
- [1ed3d4f5](https://github.com/QuantEcon/lecture-python-intro/commit/1ed3d4f5) Matt McKay 6 months ago MAINT: update anaconda=2025.06 and python=3.13 (\#637)
- [6cc6c4bc](https://github.com/QuantEcon/lecture-python-intro/commit/6cc6c4bc) Humphrey Yang 1 year ago \[FIX\] Fix string warnings (\#566)
- [2ecd0ff2](https://github.com/QuantEcon/lecture-python-intro/commit/2ecd0ff2) Longye Tian 1 year ago \[ar1\_processes\] to \[eigen\_II\] \[heavy\_tails\]\[inequality\]Example and spelling (\#538)
- [0636d75e](https://github.com/QuantEcon/lecture-python-intro/commit/0636d75e) hengchengzhang 2 years ago Fix building error
- [1dfeb593](https://github.com/QuantEcon/lecture-python-intro/commit/1dfeb593) hengchengzhang 2 years ago Minor updates
- [98a08b87](https://github.com/QuantEcon/lecture-python-intro/commit/98a08b87) hengchengzhang 2 years ago Adjust to PEP8
- [9b19b600](https://github.com/QuantEcon/lecture-python-intro/commit/9b19b600) hengchengzhang 2 years ago Minor updates
- [1d35cb53](https://github.com/QuantEcon/lecture-python-intro/commit/1d35cb53) hengchengzhang 2 years ago Fix typo
- [e49728b6](https://github.com/QuantEcon/lecture-python-intro/commit/e49728b6) mmcky 2 years ago ENH: enable linkchecker (\#174)
- [5d351b53](https://github.com/QuantEcon/lecture-python-intro/commit/5d351b53) hengchengzhang 2 years ago Updates
# 17\. Eigenvalues and Eigenvectors[\#](https://intro.quantecon.org/eigen_I.html#eigenvalues-and-eigenvectors "Link to this heading")
## 17\.1. Overview[\#](https://intro.quantecon.org/eigen_I.html#overview "Link to this heading")
Eigenvalues and eigenvectors are a relatively advanced topic in linear algebra.
At the same time, these concepts are extremely useful for
- economic modeling (especially dynamics!)
- statistics
- some parts of applied mathematics
- machine learning
- and many other fields of science.
In this lecture we explain the basics of eigenvalues and eigenvectors and introduce the Neumann Series Lemma.
We assume in this lecture that students are familiar with matrices and understand [the basics of matrix algebra](https://intro.quantecon.org/linear_equations.html).
We will use the following imports:
```
import matplotlib.pyplot as plt
import numpy as np
from numpy.linalg import matrix_power
from matplotlib.lines import Line2D
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import proj3d
```
## 17\.2. Matrices as transformations[\#](https://intro.quantecon.org/eigen_I.html#matrices-as-transformations "Link to this heading")
Let’s start by discussing an important concept concerning matrices.
### 17\.2.1. Mapping vectors to vectors[\#](https://intro.quantecon.org/eigen_I.html#mapping-vectors-to-vectors "Link to this heading")
One way to think about a matrix is as a rectangular collection of numbers.
Another way to think about a matrix is as a *map* (i.e., as a function) that transforms vectors to new vectors.
To understand the second point of view, suppose we multiply an n × m matrix A with an m × 1 column vector x to obtain an n × 1 column vector y:
A
x
\=
y
If we fix A and consider different choices of x, we can understand A as a map transforming x to A x.
Because A is n × m, it transforms m\-vectors to n\-vectors.
We can write this formally as A : R m → R n.
You might argue that if A is a function then we should write A ( x ) \= y rather than A x \= y but the second notation is more conventional.
### 17\.2.2. Square matrices[\#](https://intro.quantecon.org/eigen_I.html#square-matrices "Link to this heading")
Let’s restrict our discussion to square matrices.
In the above discussion, this means that m \= n and A maps R n to itself.
This means A is an n × n matrix that maps (or “transforms”) a vector x in R n to a new vector y \= A x also in R n.
Example 17.1
\[
2
1
−
1
1
\]
\[
1
3
\]
\=
\[
5
2
\]
Here, the matrix
A
\=
\[
2
1
−
1
1
\]
transforms the vector x \= \[ 1 3 \] to the vector y \= \[ 5 2 \].
Let’s visualize this using Python:
```
A = np.array([[2, 1],
[-1, 1]])
```
```
from math import sqrt
fig, ax = plt.subplots()
# Set the axes through the origin
for spine in ['left', 'bottom']:
ax.spines[spine].set_position('zero')
for spine in ['right', 'top']:
ax.spines[spine].set_color('none')
ax.set(xlim=(-2, 6), ylim=(-2, 4), aspect=1)
vecs = ((1, 3), (5, 2))
c = ['r', 'black']
for i, v in enumerate(vecs):
ax.annotate('', xy=v, xytext=(0, 0),
arrowprops=dict(color=c[i],
shrink=0,
alpha=0.7,
width=0.5))
ax.text(0.2 + 1, 0.2 + 3, 'x=$(1,3)$')
ax.text(0.2 + 5, 0.2 + 2, 'Ax=$(5,2)$')
ax.annotate('', xy=(sqrt(10/29) * 5, sqrt(10/29) * 2), xytext=(0, 0),
arrowprops=dict(color='purple',
shrink=0,
alpha=0.7,
width=0.5))
ax.annotate('', xy=(1, 2/5), xytext=(1/3, 1),
arrowprops={'arrowstyle': '->',
'connectionstyle': 'arc3,rad=-0.3'},
horizontalalignment='center')
ax.text(0.8, 0.8, f'θ', fontsize=14)
plt.show()
```
[](https://intro.quantecon.org/_images/7d23b5ff1fa1c168900029b65833d23823e73164c9bb8ca00e4c7958d775e168.png)
One way to understand this transformation is that A
- first rotates x by some angle θ and
- then scales it by some scalar γ to obtain the image y of x.
## 17\.3. Types of transformations[\#](https://intro.quantecon.org/eigen_I.html#types-of-transformations "Link to this heading")
Let’s examine some standard transformations we can perform with matrices.
Below we visualize transformations by thinking of vectors as points instead of arrows.
We consider how a given matrix transforms
- a grid of points and
- a set of points located on the unit circle in R 2.
To build the transformations we will use two functions, called `grid_transform` and `circle_transform`.
Each of these functions visualizes the actions of a given 2 × 2 matrix A.
Show source
Hide code cell source
```
def colorizer(x, y):
r = min(1, 1-y/3)
g = min(1, 1+y/3)
b = 1/4 + x/16
return (r, g, b)
def grid_transform(A=np.array([[1, -1], [1, 1]])):
xvals = np.linspace(-4, 4, 9)
yvals = np.linspace(-3, 3, 7)
xygrid = np.column_stack([[x, y] for x in xvals for y in yvals])
uvgrid = A @ xygrid
colors = list(map(colorizer, xygrid[0], xygrid[1]))
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
for axes in ax:
axes.set(xlim=(-11, 11), ylim=(-11, 11))
axes.set_xticks([])
axes.set_yticks([])
for spine in ['left', 'bottom']:
axes.spines[spine].set_position('zero')
for spine in ['right', 'top']:
axes.spines[spine].set_color('none')
# Plot x-y grid points
ax[0].scatter(xygrid[0], xygrid[1], s=36, c=colors, edgecolor="none")
# ax[0].grid(True)
# ax[0].axis("equal")
ax[0].set_title("points $x_1, x_2, \cdots, x_k$")
# Plot transformed grid points
ax[1].scatter(uvgrid[0], uvgrid[1], s=36, c=colors, edgecolor="none")
# ax[1].grid(True)
# ax[1].axis("equal")
ax[1].set_title("points $Ax_1, Ax_2, \cdots, Ax_k$")
plt.show()
def circle_transform(A=np.array([[-1, 2], [0, 1]])):
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
for axes in ax:
axes.set(xlim=(-4, 4), ylim=(-4, 4))
axes.set_xticks([])
axes.set_yticks([])
for spine in ['left', 'bottom']:
axes.spines[spine].set_position('zero')
for spine in ['right', 'top']:
axes.spines[spine].set_color('none')
θ = np.linspace(0, 2 * np.pi, 150)
r = 1
θ_1 = np.empty(12)
for i in range(12):
θ_1[i] = 2 * np.pi * (i/12)
x = r * np.cos(θ)
y = r * np.sin(θ)
a = r * np.cos(θ_1)
b = r * np.sin(θ_1)
a_1 = a.reshape(1, -1)
b_1 = b.reshape(1, -1)
colors = list(map(colorizer, a, b))
ax[0].plot(x, y, color='black', zorder=1)
ax[0].scatter(a_1, b_1, c=colors, alpha=1, s=60,
edgecolors='black', zorder=2)
ax[0].set_title(r"unit circle in $\mathbb{R}^2$")
x1 = x.reshape(1, -1)
y1 = y.reshape(1, -1)
ab = np.concatenate((a_1, b_1), axis=0)
transformed_ab = A @ ab
transformed_circle_input = np.concatenate((x1, y1), axis=0)
transformed_circle = A @ transformed_circle_input
ax[1].plot(transformed_circle[0, :],
transformed_circle[1, :], color='black', zorder=1)
ax[1].scatter(transformed_ab[0, :], transformed_ab[1:,],
color=colors, alpha=1, s=60, edgecolors='black', zorder=2)
ax[1].set_title("transformed circle")
plt.show()
```
### 17\.3.1. Scaling[\#](https://intro.quantecon.org/eigen_I.html#scaling "Link to this heading")
A matrix of the form
\[
α
0
0
β
\]
scales vectors across the x-axis by a factor α and along the y-axis by a factor β.
Here we illustrate a simple example where α \= β \= 3.
```
A = np.array([[3, 0], # scaling by 3 in both directions
[0, 3]])
grid_transform(A)
circle_transform(A)
```
[](https://intro.quantecon.org/_images/ec0670fe128a1127c8a21ade7dc684ffca7906f98f85553754ba44792f7a676b.png) [](https://intro.quantecon.org/_images/ba052802122543f4266caf4abc251f2b2d03910aea58d04bd117eddac946856e.png)
### 17\.3.2. Shearing[\#](https://intro.quantecon.org/eigen_I.html#shearing "Link to this heading")
A “shear” matrix of the form
\[
1
λ
0
1
\]
stretches vectors along the x-axis by an amount proportional to the y-coordinate of a point.
```
A = np.array([[1, 2], # shear along x-axis
[0, 1]])
grid_transform(A)
circle_transform(A)
```
[](https://intro.quantecon.org/_images/33b8a5c284c05b4b78d4cdb1aa3a2b30734bb6a773eac0ecb126ea02fff561a1.png) [](https://intro.quantecon.org/_images/d9f347261bdae764c86b586f84422867681f926b3ba795bbe4b8dfb46b13661b.png)
### 17\.3.3. Rotation[\#](https://intro.quantecon.org/eigen_I.html#rotation "Link to this heading")
A matrix of the form
\[
cos
θ
sin
θ
−
sin
θ
cos
θ
\]
is called a *rotation matrix*.
This matrix rotates vectors clockwise by an angle θ.
```
θ = np.pi/4 # 45 degree clockwise rotation
A = np.array([[np.cos(θ), np.sin(θ)],
[-np.sin(θ), np.cos(θ)]])
grid_transform(A)
```
[](https://intro.quantecon.org/_images/7224517b2b0b8251812eb69305b5422d9b07e7f2a92ec795ab13d3652d60dfa6.png)
### 17\.3.4. Permutation[\#](https://intro.quantecon.org/eigen_I.html#permutation "Link to this heading")
The permutation matrix
\[
0
1
1
0
\]
interchanges the coordinates of a vector.
```
A = np.column_stack([[0, 1], [1, 0]])
grid_transform(A)
```
[](https://intro.quantecon.org/_images/b30e756cf362d88e19e6a97bde7f588f3c86656dc5f64d4ef29c96966f0a227d.png)
More examples of common transition matrices can be found [here](https://en.wikipedia.org/wiki/Transformation_matrix#Examples_in_2_dimensions).
## 17\.4. Matrix multiplication as composition[\#](https://intro.quantecon.org/eigen_I.html#matrix-multiplication-as-composition "Link to this heading")
Since matrices act as functions that transform one vector to another, we can apply the concept of function composition to matrices as well.
### 17\.4.1. Linear compositions[\#](https://intro.quantecon.org/eigen_I.html#linear-compositions "Link to this heading")
Consider the two matrices
A
\=
\[
0
1
−
1
0
\]
and
B
\=
\[
1
2
0
1
\]
What will the output be when we try to obtain A B x for some 2 × 1 vector x?
\[
0
1
−
1
0
\]
⏟
A
\[
1
2
0
1
\]
⏟
B
\[
1
3
\]
⏞
x
→
\[
0
1
−
1
−
2
\]
⏟
A
B
\[
1
3
\]
⏞
x
→
\[
3
−
7
\]
⏞
y
\[
0
1
−
1
0
\]
⏟
A
\[
1
2
0
1
\]
⏟
B
\[
1
3
\]
⏞
x
→
\[
0
1
−
1
0
\]
⏟
A
\[
7
3
\]
⏞
B
x
→
\[
3
−
7
\]
⏞
y
We can observe that applying the transformation A B on the vector x is the same as first applying B on x and then applying A on the vector B x.
Thus the matrix product A B is the [composition](https://en.wikipedia.org/wiki/Function_composition) of the matrix transformations A and B
This means first apply transformation B and then transformation A.
When we matrix multiply an n × m matrix A with an m × k matrix B the obtained matrix product is an n × k matrix A B.
Thus, if A and B are transformations such that A : R m → R n and B : R k → R m, then A B transforms R k to R n.
Viewing matrix multiplication as composition of maps helps us understand why, under matrix multiplication, A B is generally not equal to B A.
(After all, when we compose functions, the order usually matters.)
### 17\.4.2. Examples[\#](https://intro.quantecon.org/eigen_I.html#examples "Link to this heading")
Let A be the 90 ∘ clockwise rotation matrix given by \[ 0 1 − 1 0 \] and let B be a shear matrix along the x-axis given by \[ 1 2 0 1 \].
We will visualize how a grid of points changes when we apply the transformation A B and then compare it with the transformation B A.
Show source
Hide code cell source
```
def grid_composition_transform(A=np.array([[1, -1], [1, 1]]),
B=np.array([[1, -1], [1, 1]])):
xvals = np.linspace(-4, 4, 9)
yvals = np.linspace(-3, 3, 7)
xygrid = np.column_stack([[x, y] for x in xvals for y in yvals])
uvgrid = B @ xygrid
abgrid = A @ uvgrid
colors = list(map(colorizer, xygrid[0], xygrid[1]))
fig, ax = plt.subplots(1, 3, figsize=(15, 5))
for axes in ax:
axes.set(xlim=(-12, 12), ylim=(-12, 12))
axes.set_xticks([])
axes.set_yticks([])
for spine in ['left', 'bottom']:
axes.spines[spine].set_position('zero')
for spine in ['right', 'top']:
axes.spines[spine].set_color('none')
# Plot grid points
ax[0].scatter(xygrid[0], xygrid[1], s=36, c=colors, edgecolor="none")
ax[0].set_title(r"points $x_1, x_2, \cdots, x_k$")
# Plot intermediate grid points
ax[1].scatter(uvgrid[0], uvgrid[1], s=36, c=colors, edgecolor="none")
ax[1].set_title(r"points $Bx_1, Bx_2, \cdots, Bx_k$")
# Plot transformed grid points
ax[2].scatter(abgrid[0], abgrid[1], s=36, c=colors, edgecolor="none")
ax[2].set_title(r"points $ABx_1, ABx_2, \cdots, ABx_k$")
plt.show()
```
```
A = np.array([[0, 1], # 90 degree clockwise rotation
[-1, 0]])
B = np.array([[1, 2], # shear along x-axis
[0, 1]])
```
#### 17\.4.2.1. Shear then rotate[\#](https://intro.quantecon.org/eigen_I.html#shear-then-rotate "Link to this heading")
```
grid_composition_transform(A, B) # transformation AB
```
[](https://intro.quantecon.org/_images/e9fbe980cdead9393a947c429afab06dff2833e0fcd8c6e916ad45aff09aa61b.png)
#### 17\.4.2.2. Rotate then shear[\#](https://intro.quantecon.org/eigen_I.html#rotate-then-shear "Link to this heading")
```
grid_composition_transform(B,A) # transformation BA
```
[](https://intro.quantecon.org/_images/c925d9e1e85a6a322bbec778a56c0f49b068ad58e929a63569c291d1b0e7714d.png)
It is evident that the transformation A B is not the same as the transformation B A.
## 17\.5. Iterating on a fixed map[\#](https://intro.quantecon.org/eigen_I.html#iterating-on-a-fixed-map "Link to this heading")
In economics (and especially in dynamic modeling), we are often interested in analyzing behavior where we repeatedly apply a fixed matrix.
For example, given a vector v and a matrix A, we are interested in studying the sequence
v
,
A
v
,
A
A
v
\=
A
2
v
,
…
Let’s first see examples of a sequence of iterates ( A k v ) k ≥ 0 under different maps A.
```
def plot_series(A, v, n):
B = np.array([[1, -1],
[1, 0]])
fig, ax = plt.subplots()
ax.set(xlim=(-4, 4), ylim=(-4, 4))
ax.set_xticks([])
ax.set_yticks([])
for spine in ['left', 'bottom']:
ax.spines[spine].set_position('zero')
for spine in ['right', 'top']:
ax.spines[spine].set_color('none')
θ = np.linspace(0, 2 * np.pi, 150)
r = 2.5
x = r * np.cos(θ)
y = r * np.sin(θ)
x1 = x.reshape(1, -1)
y1 = y.reshape(1, -1)
xy = np.concatenate((x1, y1), axis=0)
ellipse = B @ xy
ax.plot(ellipse[0, :], ellipse[1, :], color='black',
linestyle=(0, (5, 10)), linewidth=0.5)
# Initialize holder for trajectories
colors = plt.cm.rainbow(np.linspace(0, 1, 20))
for i in range(n):
iteration = matrix_power(A, i) @ v
v1 = iteration[0]
v2 = iteration[1]
ax.scatter(v1, v2, color=colors[i])
if i == 0:
ax.text(v1+0.25, v2, f'$v$')
elif i == 1:
ax.text(v1+0.25, v2, f'$Av$')
elif 1 < i < 4:
ax.text(v1+0.25, v2, f'$A^{i}v$')
plt.show()
```
```
A = np.array([[sqrt(3) + 1, -2],
[1, sqrt(3) - 1]])
A = (1/(2*sqrt(2))) * A
v = (-3, -3)
n = 12
plot_series(A, v, n)
```
[](https://intro.quantecon.org/_images/d1ccd8bebb39aa9e536e913e0debb86cf5c0792a51abb85be4640953a15680c8.png)
Here with each iteration the vectors get shorter, i.e., move closer to the origin.
In this case, repeatedly multiplying a vector by A makes the vector “spiral in”.
```
B = np.array([[sqrt(3) + 1, -2],
[1, sqrt(3) - 1]])
B = (1/2) * B
v = (2.5, 0)
n = 12
plot_series(B, v, n)
```
[](https://intro.quantecon.org/_images/d7017d7b767cbbd55360a5dc5c4b2c9cb517f41b142884606b2e6d09abe070cb.png)
Here with each iteration vectors do not tend to get longer or shorter.
In this case, repeatedly multiplying a vector by A simply “rotates it around an ellipse”.
```
B = np.array([[sqrt(3) + 1, -2],
[1, sqrt(3) - 1]])
B = (1/sqrt(2)) * B
v = (-1, -0.25)
n = 6
plot_series(B, v, n)
```
[](https://intro.quantecon.org/_images/98508a529d57df518cbf11b072dfda23b1643217eb299fdb6448349b6d304a76.png)
Here with each iteration vectors tend to get longer, i.e., farther from the origin.
In this case, repeatedly multiplying a vector by A makes the vector “spiral out”.
We thus observe that the sequence ( A k v ) k ≥ 0 behaves differently depending on the map A itself.
We now discuss the property of A that determines this behavior.
## 17\.6. Eigenvalues[\#](https://intro.quantecon.org/eigen_I.html#eigenvalues "Link to this heading")
In this section we introduce the notions of eigenvalues and eigenvectors.
### 17\.6.1. Definitions[\#](https://intro.quantecon.org/eigen_I.html#definitions "Link to this heading")
Let A be an n × n square matrix.
If λ is scalar and v is a non-zero n\-vector such that
A
v
\=
λ
v
.
Then we say that λ is an *eigenvalue* of A, and v is the corresponding *eigenvector*.
Thus, an eigenvector of A is a nonzero vector v such that when the map A is applied, v is merely scaled.
The next figure shows two eigenvectors (blue arrows) and their images under A (red arrows).
As expected, the image A v of each v is just a scaled version of the original
```
from numpy.linalg import eig
A = [[1, 2],
[2, 1]]
A = np.array(A)
evals, evecs = eig(A)
evecs = evecs[:, 0], evecs[:, 1]
fig, ax = plt.subplots(figsize=(10, 8))
# Set the axes through the origin
for spine in ['left', 'bottom']:
ax.spines[spine].set_position('zero')
for spine in ['right', 'top']:
ax.spines[spine].set_color('none')
# ax.grid(alpha=0.4)
xmin, xmax = -3, 3
ymin, ymax = -3, 3
ax.set(xlim=(xmin, xmax), ylim=(ymin, ymax))
# Plot each eigenvector
for v in evecs:
ax.annotate('', xy=v, xytext=(0, 0),
arrowprops=dict(facecolor='blue',
shrink=0,
alpha=0.6,
width=0.5))
# Plot the image of each eigenvector
for v in evecs:
v = A @ v
ax.annotate('', xy=v, xytext=(0, 0),
arrowprops=dict(facecolor='red',
shrink=0,
alpha=0.6,
width=0.5))
# Plot the lines they run through
x = np.linspace(xmin, xmax, 3)
for v in evecs:
a = v[1] / v[0]
ax.plot(x, a * x, 'b-', lw=0.4)
plt.show()
```
[](https://intro.quantecon.org/_images/6d26bc96ac56471100124f54422cdb6e6e7b1a1fa712882895d0268a91a113b9.png)
### 17\.6.2. Complex values[\#](https://intro.quantecon.org/eigen_I.html#complex-values "Link to this heading")
So far our definition of eigenvalues and eigenvectors seems straightforward.
There is one complication we haven’t mentioned yet:
When solving A v \= λ v,
- λ is allowed to be a complex number and
- v is allowed to be an n\-vector of complex numbers.
We will see some examples below.
### 17\.6.3. Some mathematical details[\#](https://intro.quantecon.org/eigen_I.html#some-mathematical-details "Link to this heading")
We note some mathematical details for more advanced readers.
(Other readers can skip to the next section.)
The eigenvalue equation is equivalent to ( A − λ I ) v \= 0.
This equation has a nonzero solution v only when the columns of A − λ I are linearly dependent.
This in turn is equivalent to stating the determinant is zero.
Hence, to find all eigenvalues, we can look for λ such that the determinant of A − λ I is zero.
This problem can be expressed as one of solving for the roots of a polynomial in λ of degree n.
This in turn implies the existence of n solutions in the complex plane, although some might be repeated.
### 17\.6.4. Facts[\#](https://intro.quantecon.org/eigen_I.html#facts "Link to this heading")
Some nice facts about the eigenvalues of a square matrix A are as follows:
1. the determinant of A equals the product of the eigenvalues
2. the trace of A (the sum of the elements on the principal diagonal) equals the sum of the eigenvalues
3. if A is symmetric, then all of its eigenvalues are real
4. if A is invertible and λ 1 , … , λ n are its eigenvalues, then the eigenvalues of A − 1 are 1 / λ 1 , … , 1 / λ n.
A corollary of the last statement is that a matrix is invertible if and only if all its eigenvalues are nonzero.
### 17\.6.5. Computation[\#](https://intro.quantecon.org/eigen_I.html#computation "Link to this heading")
Using NumPy, we can solve for the eigenvalues and eigenvectors of a matrix as follows
```
from numpy.linalg import eig
A = ((1, 2),
(2, 1))
A = np.array(A)
evals, evecs = eig(A)
evals # eigenvalues
```
```
array([ 3., -1.])
```
```
evecs # eigenvectors
```
```
array([[ 0.70710678, -0.70710678],
[ 0.70710678, 0.70710678]])
```
Note that the *columns* of `evecs` are the eigenvectors.
Since any scalar multiple of an eigenvector is an eigenvector with the same eigenvalue (which can be verified), the `eig` routine normalizes the length of each eigenvector to one.
The eigenvectors and eigenvalues of a map A determine how a vector v is transformed when we repeatedly multiply by A.
This is discussed further later.
## 17\.7. The Neumann Series Lemma[\#](https://intro.quantecon.org/eigen_I.html#the-neumann-series-lemma "Link to this heading")
In this section we present a famous result about series of matrices that has many applications in economics.
### 17\.7.1. Scalar series[\#](https://intro.quantecon.org/eigen_I.html#scalar-series "Link to this heading")
Here’s a fundamental result about series:
If a is a number and \| a \| \< 1, then
(17.1)[\#](https://intro.quantecon.org/eigen_I.html#equation-gp-sum "Link to this equation")
∑
k
\=
0
∞
a
k
\=
1
1
−
a
\=
(
1
−
a
)
−
1
For a one-dimensional linear equation x \= a x \+ b where x is unknown we can thus conclude that the solution x ∗ is given by:
x
∗
\=
b
1
−
a
\=
∑
k
\=
0
∞
a
k
b
### 17\.7.2. Matrix series[\#](https://intro.quantecon.org/eigen_I.html#matrix-series "Link to this heading")
A generalization of this idea exists in the matrix setting.
Consider the system of equations x \= A x \+ b where A is an n × n square matrix and x and b are both column vectors in R n.
Using matrix algebra we can conclude that the solution to this system of equations will be given by:
(17.2)[\#](https://intro.quantecon.org/eigen_I.html#equation-neumann-eqn "Link to this equation")
x
∗
\=
(
I
−
A
)
−
1
b
What guarantees the existence of a unique vector x ∗ that satisfies [(17.2)](https://intro.quantecon.org/eigen_I.html#equation-neumann-eqn)?
The following is a fundamental result in functional analysis that generalizes [(17.1)](https://intro.quantecon.org/eigen_I.html#equation-gp-sum) to a multivariate case.
Theorem 17.1 (Neumann Series Lemma)
Let A be a square matrix and let A k be the k\-th power of A.
Let r ( A ) be the **spectral radius** of A, defined as max i \| λ i \|, where
- { λ i } i is the set of eigenvalues of A and
- \| λ i \| is the modulus of the complex number λ i
Neumann’s Theorem states the following: If r ( A ) \< 1, then I − A is invertible, and
(
I
−
A
)
−
1
\=
∑
k
\=
0
∞
A
k
We can see the Neumann Series Lemma in action in the following example.
```
A = np.array([[0.4, 0.1],
[0.7, 0.2]])
evals, evecs = eig(A) # finding eigenvalues and eigenvectors
r = max(abs(λ) for λ in evals) # compute spectral radius
print(r)
```
```
0.5828427124746189
```
The spectral radius r ( A ) obtained is less than 1.
Thus, we can apply the Neumann Series Lemma to find ( I − A ) − 1.
```
I = np.identity(2) # 2 x 2 identity matrix
B = I - A
```
```
B_inverse = np.linalg.inv(B) # direct inverse method
```
```
A_sum = np.zeros((2, 2)) # power series sum of A
A_power = I
for i in range(50):
A_sum += A_power
A_power = A_power @ A
```
Let’s check equality between the sum and the inverse methods.
```
np.allclose(A_sum, B_inverse)
```
```
True
```
Although we truncate the infinite sum at k \= 50, both methods give us the same result which illustrates the result of the Neumann Series Lemma.
## 17\.8. Exercises[\#](https://intro.quantecon.org/eigen_I.html#exercises "Link to this heading")
Exercise 17.1
Power iteration is a method for finding the greatest absolute eigenvalue of a diagonalizable matrix.
The method starts with a random vector b 0 and repeatedly applies the matrix A to it
b
k
\+
1
\=
A
b
k
‖
A
b
k
‖
A thorough discussion of the method can be found [here](https://pythonnumericalmethods.berkeley.edu/notebooks/chapter15.02-The-Power-Method.html).
In this exercise, first implement the power iteration method and use it to find the greatest absolute eigenvalue and its corresponding eigenvector.
Then visualize the convergence.
Solution to [Exercise 17.1](https://intro.quantecon.org/eigen_I.html#eig1_ex1)
Here is one solution.
We start by looking into the distance between the eigenvector approximation and the true eigenvector.
```
# Define a matrix A
A = np.array([[1, 0, 3],
[0, 2, 0],
[3, 0, 1]])
num_iters = 20
# Define a random starting vector b
b = np.random.rand(A.shape[1])
# Get the leading eigenvector of matrix A
eigenvector = np.linalg.eig(A)[1][:, 0]
errors = []
res = []
# Power iteration loop
for i in range(num_iters):
# Multiply b by A
b = A @ b
# Normalize b
b = b / np.linalg.norm(b)
# Append b to the list of eigenvector approximations
res.append(b)
err = np.linalg.norm(np.array(b)
- eigenvector)
errors.append(err)
greatest_eigenvalue = np.dot(A @ b, b) / np.dot(b, b)
print(f'The approximated greatest absolute eigenvalue is \
{greatest_eigenvalue:.2f}')
print('The real eigenvalue is', np.linalg.eig(A)[0])
# Plot the eigenvector approximations for each iteration
plt.figure(figsize=(10, 6))
plt.xlabel('iterations')
plt.ylabel('error')
plt.title('Power iteration')
_ = plt.plot(errors)
```
```
The approximated greatest absolute eigenvalue is 4.00
The real eigenvalue is [ 4. -2. 2.]
```
[](https://intro.quantecon.org/_images/51fef2c19a430458673ba258a12f0e6d86ccebadb77e0708ae8df5128adec36a.png)
Then we can look at the trajectory of the eigenvector approximation.
```
# Set up the figure and axis for 3D plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# Plot the eigenvectors
ax.scatter(eigenvector[0],
eigenvector[1],
eigenvector[2],
color='r', s=80)
for i, vec in enumerate(res):
ax.scatter(vec[0], vec[1], vec[2],
color='b',
alpha=(i+1)/(num_iters+1),
s=80)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.tick_params(axis='both', which='major', labelsize=7)
points = [plt.Line2D([0], [0], linestyle='none',
c=i, marker='o') for i in ['r', 'b']]
ax.legend(points, ['actual eigenvector',
r'approximated eigenvector ($b_k$)'])
ax.set_box_aspect(aspect=None, zoom=0.8)
ax.set_title('Power iteration trajectory')
plt.show()
```
[](https://intro.quantecon.org/_images/3d07750e23b7689140205ef32dd772d0491f8faeda1a389c10836101a1ca5a0a.png)
Exercise 17.2
We have discussed the trajectory of the vector v after being transformed by A.
Consider the matrix A \= \[ 1 2 1 1 \] and the vector v \= \[ 2 − 2 \].
Try to compute the trajectory of v after being transformed by A for n \= 4 iterations and plot the result.
Solution to [Exercise 17.2](https://intro.quantecon.org/eigen_I.html#eig1_ex2)
```
A = np.array([[1, 2],
[1, 1]])
v = (0.4, -0.4)
n = 11
# Compute eigenvectors and eigenvalues
eigenvalues, eigenvectors = np.linalg.eig(A)
print(f'eigenvalues:\n {eigenvalues}')
print(f'eigenvectors:\n {eigenvectors}')
plot_series(A, v, n)
```
```
eigenvalues:
[ 2.41421356 -0.41421356]
eigenvectors:
[[ 0.81649658 -0.81649658]
[ 0.57735027 0.57735027]]
```
[](https://intro.quantecon.org/_images/8024ccd681013d0450608b1b83b3f16dd85603792b2ef5b3046748b275a47b67.png)
The result seems to converge to the eigenvector of A with the largest eigenvalue.
Let’s use a [vector field](https://en.wikipedia.org/wiki/Vector_field) to visualize the transformation brought by A.
(This is a more advanced topic in linear algebra, please step ahead if you are comfortable with the math.)
```
# Create a grid of points
x, y = np.meshgrid(np.linspace(-5, 5, 15),
np.linspace(-5, 5, 20))
# Apply the matrix A to each point in the vector field
vec_field = np.stack([x, y])
u, v = np.tensordot(A, vec_field, axes=1)
# Plot the transformed vector field
c = plt.streamplot(x, y, u - x, v - y,
density=1, linewidth=None, color='#A23BEC')
c.lines.set_alpha(0.5)
c.arrows.set_alpha(0.5)
# Draw eigenvectors
origin = np.zeros((2, len(eigenvectors)))
parameters = {'color': ['b', 'g'], 'angles': 'xy',
'scale_units': 'xy', 'scale': 0.1, 'width': 0.01}
plt.quiver(*origin, eigenvectors[0],
eigenvectors[1], **parameters)
plt.quiver(*origin, - eigenvectors[0],
- eigenvectors[1], **parameters)
colors = ['b', 'g']
lines = [Line2D([0], [0], color=c, linewidth=3) for c in colors]
labels = ["2.4 eigenspace", "0.4 eigenspace"]
plt.legend(lines, labels, loc='center left',
bbox_to_anchor=(1, 0.5))
plt.xlabel("x")
plt.ylabel("y")
plt.title("Convergence towards eigenvectors")
plt.grid()
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
```
[](https://intro.quantecon.org/_images/ffafeb31eda11db6138a7d19b1f3b83d3164cb11eb56a6e5920f5710c0a071ff.png)
Note that the vector field converges to the eigenvector of A with the largest eigenvalue and diverges from the eigenvector of A with the smallest eigenvalue.
In fact, the eigenvectors are also the directions in which the matrix A stretches or shrinks the space.
Specifically, the eigenvector with the largest eigenvalue is the direction in which the matrix A stretches the space the most.
We will see more intriguing examples in the following exercise.
Exercise 17.3
[Previously](https://intro.quantecon.org/eigen_I.html#plot-series), we demonstrated the trajectory of the vector v after being transformed by A for three different matrices.
Use the visualization in the previous exercise to explain the trajectory of the vector v after being transformed by A for the three different matrices.
Solution to [Exercise 17.3](https://intro.quantecon.org/eigen_I.html#eig1_ex3)
Here is one solution
```
fig, ax = plt.subplots(1, 3, figsize=(15, 5))
A = np.array([[sqrt(3) + 1, -2],
[1, sqrt(3) - 1]])
A = (1/(2*sqrt(2))) * A
B = np.array([[sqrt(3) + 1, -2],
[1, sqrt(3) - 1]])
B = (1/2) * B
C = np.array([[sqrt(3) + 1, -2],
[1, sqrt(3) - 1]])
C = (1/sqrt(2)) * C
examples = [A, B, C]
for i, example in enumerate(examples):
M = example
# Compute right eigenvectors and eigenvalues
eigenvalues, eigenvectors = np.linalg.eig(M)
print(f'Example {i+1}:\n')
print(f'eigenvalues:\n {eigenvalues}')
print(f'eigenvectors:\n {eigenvectors}\n')
eigenvalues_real = eigenvalues.real
eigenvectors_real = eigenvectors.real
# Create a grid of points
x, y = np.meshgrid(np.linspace(-20, 20, 15),
np.linspace(-20, 20, 20))
# Apply the matrix A to each point in the vector field
vec_field = np.stack([x, y])
u, v = np.tensordot(M, vec_field, axes=1)
# Plot the transformed vector field
c = ax[i].streamplot(x, y, u - x, v - y, density=1,
linewidth=None, color='#A23BEC')
c.lines.set_alpha(0.5)
c.arrows.set_alpha(0.5)
# Draw eigenvectors
parameters = {'color': ['b', 'g'], 'angles': 'xy',
'scale_units': 'xy', 'scale': 1,
'width': 0.01, 'alpha': 0.5}
origin = np.zeros((2, len(eigenvectors)))
ax[i].quiver(*origin, eigenvectors_real[0],
eigenvectors_real[1], **parameters)
ax[i].quiver(*origin,
- eigenvectors_real[0],
- eigenvectors_real[1],
**parameters)
ax[i].set_xlabel("x-axis")
ax[i].set_ylabel("y-axis")
ax[i].grid()
ax[i].set_aspect('equal', adjustable='box')
fig.suptitle("Vector fields of the three matrices")
plt.show()
```
```
Example 1:
eigenvalues:
[0.61237244+0.35355339j 0.61237244-0.35355339j]
eigenvectors:
[[0.81649658+0.j 0.81649658-0.j ]
[0.40824829-0.40824829j 0.40824829+0.40824829j]]
Example 2:
eigenvalues:
[0.8660254+0.5j 0.8660254-0.5j]
eigenvectors:
[[0.81649658+0.j 0.81649658-0.j ]
[0.40824829-0.40824829j 0.40824829+0.40824829j]]
Example 3:
eigenvalues:
[1.22474487+0.70710678j 1.22474487-0.70710678j]
eigenvectors:
[[0.81649658+0.j 0.81649658-0.j ]
[0.40824829-0.40824829j 0.40824829+0.40824829j]]
```
[](https://intro.quantecon.org/_images/f5b355bc139671ffc6677677600d5a2df12e79a964d52b6dddedee42141fd7c7.png)
The vector fields explain why we observed the trajectories of the vector v multiplied by A iteratively before.
The pattern demonstrated here is because we have complex eigenvalues and eigenvectors.
We can plot the complex plane for one of the matrices using `Arrow3D` class retrieved from [stackoverflow](https://stackoverflow.com/questions/22867620/putting-arrowheads-on-vectors-in-a-3d-plot).
```
class Arrow3D(FancyArrowPatch):
def __init__(self, xs, ys, zs, *args, **kwargs):
super().__init__((0, 0), (0, 0), *args, **kwargs)
self._verts3d = xs, ys, zs
def do_3d_projection(self):
xs3d, ys3d, zs3d = self._verts3d
xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d,
self.axes.M)
self.set_positions((0.1*xs[0], 0.1*ys[0]),
(0.1*xs[1], 0.1*ys[1]))
return np.min(zs)
eigenvalues, eigenvectors = np.linalg.eig(A)
# Create meshgrid for vector field
x, y = np.meshgrid(np.linspace(-2, 2, 15),
np.linspace(-2, 2, 15))
# Calculate vector field (real and imaginary parts)
u_real = A[0][0] * x + A[0][1] * y
v_real = A[1][0] * x + A[1][1] * y
u_imag = np.zeros_like(x)
v_imag = np.zeros_like(y)
# Create 3D figure
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
vlength = np.linalg.norm(eigenvectors)
ax.quiver(x, y, u_imag, u_real-x, v_real-y, v_imag-u_imag,
colors='b', alpha=0.3, length=.2,
arrow_length_ratio=0.01)
arrow_prop_dict = dict(mutation_scale=5,
arrowstyle='-|>', shrinkA=0, shrinkB=0)
# Plot 3D eigenvectors
for c, i in zip(['b', 'g'], [0, 1]):
a = Arrow3D([0, eigenvectors[0][i].real],
[0, eigenvectors[1][i].real],
[0, eigenvectors[1][i].imag],
color=c, **arrow_prop_dict)
ax.add_artist(a)
# Set axis labels and title
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('Im')
ax.set_box_aspect(aspect=None, zoom=0.8)
plt.title("3D plot of the vector field")
plt.draw()
plt.show()
```
[](https://intro.quantecon.org/_images/b75bf6c2b8cd93365476b2c7e8c16b483dc78a7a0cfb84588aba9c318a9c72ae.png)
[](https://creativecommons.org/licenses/by-sa/4.0/)
Creative Commons License – This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International.
A theme by [QuantEcon](https://quantecon.org/)
Contents
Introduction
- [1\. About These Lectures](https://intro.quantecon.org/about.html)
Economic Data
- [2\. Long-Run Growth](https://intro.quantecon.org/long_run_growth.html)
- [3\. Business Cycles](https://intro.quantecon.org/business_cycle.html)
- [4\. Price Level Histories](https://intro.quantecon.org/inflation_history.html)
- [5\. Inflation During French Revolution](https://intro.quantecon.org/french_rev.html)
- [6\. Income and Wealth Inequality](https://intro.quantecon.org/inequality.html)
Foundations
- [7\. Introduction to Supply and Demand](https://intro.quantecon.org/intro_supply_demand.html)
- [8\. Linear Equations and Matrix Algebra](https://intro.quantecon.org/linear_equations.html)
- [9\. Complex Numbers and Trigonometry](https://intro.quantecon.org/complex_and_trig.html)
- [10\. Geometric Series for Elementary Economics](https://intro.quantecon.org/geom_series.html)
Linear Dynamics: Finite Horizons
- [11\. Present Values](https://intro.quantecon.org/pv.html)
- [12\. Consumption Smoothing](https://intro.quantecon.org/cons_smooth.html)
- [13\. Tax Smoothing](https://intro.quantecon.org/tax_smooth.html)
- [14\. Equalizing Difference Model](https://intro.quantecon.org/equalizing_difference.html)
- [15\. A Monetarist Theory of Price Levels](https://intro.quantecon.org/cagan_ree.html)
- [16\. Monetarist Theory of Price Levels with Adaptive Expectations](https://intro.quantecon.org/cagan_adaptive.html)
Linear Dynamics: Infinite Horizons
- [17\. Eigenvalues and Eigenvectors](https://intro.quantecon.org/eigen_I.html)
- [18\. Computing Square Roots](https://intro.quantecon.org/greek_square.html)
Probability and Distributions
- [19\. Distributions and Probabilities](https://intro.quantecon.org/prob_dist.html)
- [20\. LLN and CLT](https://intro.quantecon.org/lln_clt.html)
- [21\. Monte Carlo and Option Pricing](https://intro.quantecon.org/monte_carlo.html)
- [22\. Heavy-Tailed Distributions](https://intro.quantecon.org/heavy_tails.html)
- [23\. Racial Segregation](https://intro.quantecon.org/schelling.html)
Nonlinear Dynamics
- [24\. Dynamics in One Dimension](https://intro.quantecon.org/scalar_dynam.html)
- [25\. The Solow-Swan Growth Model](https://intro.quantecon.org/solow.html)
- [26\. The Cobweb Model](https://intro.quantecon.org/cobweb.html)
- [27\. The Overlapping Generations Model](https://intro.quantecon.org/olg.html)
- [28\. Commodity Prices](https://intro.quantecon.org/commod_price.html)
Monetary-Fiscal Policy Interactions
- [29\. Money Financed Government Deficits and Price Levels](https://intro.quantecon.org/money_inflation.html)
- [30\. Some Unpleasant Monetarist Arithmetic](https://intro.quantecon.org/unpleasant.html)
- [31\. Inflation Rate Laffer Curves](https://intro.quantecon.org/money_inflation_nonlinear.html)
- [32\. Laffer Curves with Adaptive Expectations](https://intro.quantecon.org/laffer_adaptive.html)
Stochastic Dynamics
- [33\. AR(1) Processes](https://intro.quantecon.org/ar1_processes.html)
- [34\. Markov Chains: Basic Concepts](https://intro.quantecon.org/markov_chains_I.html)
- [35\. Markov Chains: Irreducibility and Ergodicity](https://intro.quantecon.org/markov_chains_II.html)
- [36\. Univariate Time Series with Matrix Algebra](https://intro.quantecon.org/time_series_with_matrices.html)
Optimization
- [37\. Linear Programming](https://intro.quantecon.org/lp_intro.html)
- [38\. Shortest Paths](https://intro.quantecon.org/short_path.html)
Modeling in Higher Dimensions
- [39\. The Perron-Frobenius Theorem](https://intro.quantecon.org/eigen_II.html)
- [40\. Input-Output Models](https://intro.quantecon.org/input_output.html)
- [41\. A Lake Model of Employment](https://intro.quantecon.org/lake_model.html)
- [42\. Networks](https://intro.quantecon.org/networks.html)
Markets and Competitive Equilibrium
- [43\. Supply and Demand with Many Goods](https://intro.quantecon.org/supply_demand_multiple_goods.html)
- [44\. Market Equilibrium with Heterogeneity](https://intro.quantecon.org/supply_demand_heterogeneity.html)
Estimation
- [45\. Simple Linear Regression Model](https://intro.quantecon.org/simple_linear_regression.html)
- [46\. Maximum Likelihood Estimation](https://intro.quantecon.org/mle.html)
Other
- [47\. Troubleshooting](https://intro.quantecon.org/troubleshooting.html)
- [48\. References](https://intro.quantecon.org/zreferences.html)
- [49\. Execution Statistics](https://intro.quantecon.org/status.html)
- [QuantEcon](https://quantecon.org/) |
| Readable Markdown | ## 17\.1. Overview[\#](https://intro.quantecon.org/eigen_I.html#overview "Link to this heading")
Eigenvalues and eigenvectors are a relatively advanced topic in linear algebra.
At the same time, these concepts are extremely useful for
- economic modeling (especially dynamics!)
- statistics
- some parts of applied mathematics
- machine learning
- and many other fields of science.
In this lecture we explain the basics of eigenvalues and eigenvectors and introduce the Neumann Series Lemma.
We assume in this lecture that students are familiar with matrices and understand [the basics of matrix algebra](https://intro.quantecon.org/linear_equations.html).
We will use the following imports:
```
import matplotlib.pyplot as plt
import numpy as np
from numpy.linalg import matrix_power
from matplotlib.lines import Line2D
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import proj3d
```
## 17\.2. Matrices as transformations[\#](https://intro.quantecon.org/eigen_I.html#matrices-as-transformations "Link to this heading")
Let’s start by discussing an important concept concerning matrices.
### 17\.2.1. Mapping vectors to vectors[\#](https://intro.quantecon.org/eigen_I.html#mapping-vectors-to-vectors "Link to this heading")
One way to think about a matrix is as a rectangular collection of numbers.
Another way to think about a matrix is as a *map* (i.e., as a function) that transforms vectors to new vectors.
To understand the second point of view, suppose we multiply an n × m matrix A with an m × 1 column vector x to obtain an n × 1 column vector y:
A x \= y
If we fix A and consider different choices of x, we can understand A as a map transforming x to A x.
Because A is n × m, it transforms m\-vectors to n\-vectors.
We can write this formally as A : R m → R n.
You might argue that if A is a function then we should write A ( x ) \= y rather than A x \= y but the second notation is more conventional.
### 17\.2.2. Square matrices[\#](https://intro.quantecon.org/eigen_I.html#square-matrices "Link to this heading")
Let’s restrict our discussion to square matrices.
In the above discussion, this means that m \= n and A maps R n to itself.
This means A is an n × n matrix that maps (or “transforms”) a vector x in R n to a new vector y \= A x also in R n.
Example 17.1
\[ 2 1 − 1 1 \] \[ 1 3 \] \= \[ 5 2 \]
Here, the matrix
A \= \[ 2 1 − 1 1 \]
transforms the vector x \= \[ 1 3 \] to the vector y \= \[ 5 2 \].
Let’s visualize this using Python:
```
A = np.array([[2, 1],
[-1, 1]])
```
```
from math import sqrt
fig, ax = plt.subplots()
# Set the axes through the origin
for spine in ['left', 'bottom']:
ax.spines[spine].set_position('zero')
for spine in ['right', 'top']:
ax.spines[spine].set_color('none')
ax.set(xlim=(-2, 6), ylim=(-2, 4), aspect=1)
vecs = ((1, 3), (5, 2))
c = ['r', 'black']
for i, v in enumerate(vecs):
ax.annotate('', xy=v, xytext=(0, 0),
arrowprops=dict(color=c[i],
shrink=0,
alpha=0.7,
width=0.5))
ax.text(0.2 + 1, 0.2 + 3, 'x=$(1,3)$')
ax.text(0.2 + 5, 0.2 + 2, 'Ax=$(5,2)$')
ax.annotate('', xy=(sqrt(10/29) * 5, sqrt(10/29) * 2), xytext=(0, 0),
arrowprops=dict(color='purple',
shrink=0,
alpha=0.7,
width=0.5))
ax.annotate('', xy=(1, 2/5), xytext=(1/3, 1),
arrowprops={'arrowstyle': '->',
'connectionstyle': 'arc3,rad=-0.3'},
horizontalalignment='center')
ax.text(0.8, 0.8, f'θ', fontsize=14)
plt.show()
```
[](https://intro.quantecon.org/_images/7d23b5ff1fa1c168900029b65833d23823e73164c9bb8ca00e4c7958d775e168.png)
One way to understand this transformation is that A
- first rotates x by some angle θ and
- then scales it by some scalar γ to obtain the image y of x.
## 17\.3. Types of transformations[\#](https://intro.quantecon.org/eigen_I.html#types-of-transformations "Link to this heading")
Let’s examine some standard transformations we can perform with matrices.
Below we visualize transformations by thinking of vectors as points instead of arrows.
We consider how a given matrix transforms
- a grid of points and
- a set of points located on the unit circle in R 2.
To build the transformations we will use two functions, called `grid_transform` and `circle_transform`.
Each of these functions visualizes the actions of a given 2 × 2 matrix A.
Show source
Hide code cell source
```
def colorizer(x, y):
r = min(1, 1-y/3)
g = min(1, 1+y/3)
b = 1/4 + x/16
return (r, g, b)
def grid_transform(A=np.array([[1, -1], [1, 1]])):
xvals = np.linspace(-4, 4, 9)
yvals = np.linspace(-3, 3, 7)
xygrid = np.column_stack([[x, y] for x in xvals for y in yvals])
uvgrid = A @ xygrid
colors = list(map(colorizer, xygrid[0], xygrid[1]))
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
for axes in ax:
axes.set(xlim=(-11, 11), ylim=(-11, 11))
axes.set_xticks([])
axes.set_yticks([])
for spine in ['left', 'bottom']:
axes.spines[spine].set_position('zero')
for spine in ['right', 'top']:
axes.spines[spine].set_color('none')
# Plot x-y grid points
ax[0].scatter(xygrid[0], xygrid[1], s=36, c=colors, edgecolor="none")
# ax[0].grid(True)
# ax[0].axis("equal")
ax[0].set_title("points $x_1, x_2, \cdots, x_k$")
# Plot transformed grid points
ax[1].scatter(uvgrid[0], uvgrid[1], s=36, c=colors, edgecolor="none")
# ax[1].grid(True)
# ax[1].axis("equal")
ax[1].set_title("points $Ax_1, Ax_2, \cdots, Ax_k$")
plt.show()
def circle_transform(A=np.array([[-1, 2], [0, 1]])):
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
for axes in ax:
axes.set(xlim=(-4, 4), ylim=(-4, 4))
axes.set_xticks([])
axes.set_yticks([])
for spine in ['left', 'bottom']:
axes.spines[spine].set_position('zero')
for spine in ['right', 'top']:
axes.spines[spine].set_color('none')
θ = np.linspace(0, 2 * np.pi, 150)
r = 1
θ_1 = np.empty(12)
for i in range(12):
θ_1[i] = 2 * np.pi * (i/12)
x = r * np.cos(θ)
y = r * np.sin(θ)
a = r * np.cos(θ_1)
b = r * np.sin(θ_1)
a_1 = a.reshape(1, -1)
b_1 = b.reshape(1, -1)
colors = list(map(colorizer, a, b))
ax[0].plot(x, y, color='black', zorder=1)
ax[0].scatter(a_1, b_1, c=colors, alpha=1, s=60,
edgecolors='black', zorder=2)
ax[0].set_title(r"unit circle in $\mathbb{R}^2$")
x1 = x.reshape(1, -1)
y1 = y.reshape(1, -1)
ab = np.concatenate((a_1, b_1), axis=0)
transformed_ab = A @ ab
transformed_circle_input = np.concatenate((x1, y1), axis=0)
transformed_circle = A @ transformed_circle_input
ax[1].plot(transformed_circle[0, :],
transformed_circle[1, :], color='black', zorder=1)
ax[1].scatter(transformed_ab[0, :], transformed_ab[1:,],
color=colors, alpha=1, s=60, edgecolors='black', zorder=2)
ax[1].set_title("transformed circle")
plt.show()
```
### 17\.3.1. Scaling[\#](https://intro.quantecon.org/eigen_I.html#scaling "Link to this heading")
A matrix of the form
\[ α 0 0 β \]
scales vectors across the x-axis by a factor α and along the y-axis by a factor β.
Here we illustrate a simple example where α \= β \= 3.
```
A = np.array([[3, 0], # scaling by 3 in both directions
[0, 3]])
grid_transform(A)
circle_transform(A)
```
[](https://intro.quantecon.org/_images/ec0670fe128a1127c8a21ade7dc684ffca7906f98f85553754ba44792f7a676b.png) [](https://intro.quantecon.org/_images/ba052802122543f4266caf4abc251f2b2d03910aea58d04bd117eddac946856e.png)
### 17\.3.2. Shearing[\#](https://intro.quantecon.org/eigen_I.html#shearing "Link to this heading")
A “shear” matrix of the form
\[ 1 λ 0 1 \]
stretches vectors along the x-axis by an amount proportional to the y-coordinate of a point.
```
A = np.array([[1, 2], # shear along x-axis
[0, 1]])
grid_transform(A)
circle_transform(A)
```
[](https://intro.quantecon.org/_images/33b8a5c284c05b4b78d4cdb1aa3a2b30734bb6a773eac0ecb126ea02fff561a1.png) [](https://intro.quantecon.org/_images/d9f347261bdae764c86b586f84422867681f926b3ba795bbe4b8dfb46b13661b.png)
### 17\.3.3. Rotation[\#](https://intro.quantecon.org/eigen_I.html#rotation "Link to this heading")
A matrix of the form
\[ cos θ sin θ − sin θ cos θ \]
is called a *rotation matrix*.
This matrix rotates vectors clockwise by an angle θ.
```
θ = np.pi/4 # 45 degree clockwise rotation
A = np.array([[np.cos(θ), np.sin(θ)],
[-np.sin(θ), np.cos(θ)]])
grid_transform(A)
```
[](https://intro.quantecon.org/_images/7224517b2b0b8251812eb69305b5422d9b07e7f2a92ec795ab13d3652d60dfa6.png)
### 17\.3.4. Permutation[\#](https://intro.quantecon.org/eigen_I.html#permutation "Link to this heading")
The permutation matrix
\[ 0 1 1 0 \]
interchanges the coordinates of a vector.
```
A = np.column_stack([[0, 1], [1, 0]])
grid_transform(A)
```
[](https://intro.quantecon.org/_images/b30e756cf362d88e19e6a97bde7f588f3c86656dc5f64d4ef29c96966f0a227d.png)
More examples of common transition matrices can be found [here](https://en.wikipedia.org/wiki/Transformation_matrix#Examples_in_2_dimensions).
## 17\.4. Matrix multiplication as composition[\#](https://intro.quantecon.org/eigen_I.html#matrix-multiplication-as-composition "Link to this heading")
Since matrices act as functions that transform one vector to another, we can apply the concept of function composition to matrices as well.
### 17\.4.1. Linear compositions[\#](https://intro.quantecon.org/eigen_I.html#linear-compositions "Link to this heading")
Consider the two matrices
A \= \[ 0 1 − 1 0 \] and B \= \[ 1 2 0 1 \]
What will the output be when we try to obtain A B x for some 2 × 1 vector x?
\[ 0 1 − 1 0 \] ⏟ A \[ 1 2 0 1 \] ⏟ B \[ 1 3 \] ⏞ x → \[ 0 1 − 1 − 2 \] ⏟ A B \[ 1 3 \] ⏞ x → \[ 3 − 7 \] ⏞ y
\[ 0 1 − 1 0 \] ⏟ A \[ 1 2 0 1 \] ⏟ B \[ 1 3 \] ⏞ x → \[ 0 1 − 1 0 \] ⏟ A \[ 7 3 \] ⏞ B x → \[ 3 − 7 \] ⏞ y
We can observe that applying the transformation A B on the vector x is the same as first applying B on x and then applying A on the vector B x.
Thus the matrix product A B is the [composition](https://en.wikipedia.org/wiki/Function_composition) of the matrix transformations A and B
This means first apply transformation B and then transformation A.
When we matrix multiply an n × m matrix A with an m × k matrix B the obtained matrix product is an n × k matrix A B.
Thus, if A and B are transformations such that A : R m → R n and B : R k → R m, then A B transforms R k to R n.
Viewing matrix multiplication as composition of maps helps us understand why, under matrix multiplication, A B is generally not equal to B A.
(After all, when we compose functions, the order usually matters.)
### 17\.4.2. Examples[\#](https://intro.quantecon.org/eigen_I.html#examples "Link to this heading")
Let A be the 90 ∘ clockwise rotation matrix given by \[ 0 1 − 1 0 \] and let B be a shear matrix along the x-axis given by \[ 1 2 0 1 \].
We will visualize how a grid of points changes when we apply the transformation A B and then compare it with the transformation B A.
Show source
Hide code cell source
```
def grid_composition_transform(A=np.array([[1, -1], [1, 1]]),
B=np.array([[1, -1], [1, 1]])):
xvals = np.linspace(-4, 4, 9)
yvals = np.linspace(-3, 3, 7)
xygrid = np.column_stack([[x, y] for x in xvals for y in yvals])
uvgrid = B @ xygrid
abgrid = A @ uvgrid
colors = list(map(colorizer, xygrid[0], xygrid[1]))
fig, ax = plt.subplots(1, 3, figsize=(15, 5))
for axes in ax:
axes.set(xlim=(-12, 12), ylim=(-12, 12))
axes.set_xticks([])
axes.set_yticks([])
for spine in ['left', 'bottom']:
axes.spines[spine].set_position('zero')
for spine in ['right', 'top']:
axes.spines[spine].set_color('none')
# Plot grid points
ax[0].scatter(xygrid[0], xygrid[1], s=36, c=colors, edgecolor="none")
ax[0].set_title(r"points $x_1, x_2, \cdots, x_k$")
# Plot intermediate grid points
ax[1].scatter(uvgrid[0], uvgrid[1], s=36, c=colors, edgecolor="none")
ax[1].set_title(r"points $Bx_1, Bx_2, \cdots, Bx_k$")
# Plot transformed grid points
ax[2].scatter(abgrid[0], abgrid[1], s=36, c=colors, edgecolor="none")
ax[2].set_title(r"points $ABx_1, ABx_2, \cdots, ABx_k$")
plt.show()
```
```
A = np.array([[0, 1], # 90 degree clockwise rotation
[-1, 0]])
B = np.array([[1, 2], # shear along x-axis
[0, 1]])
```
#### 17\.4.2.1. Shear then rotate[\#](https://intro.quantecon.org/eigen_I.html#shear-then-rotate "Link to this heading")
```
grid_composition_transform(A, B) # transformation AB
```
[](https://intro.quantecon.org/_images/e9fbe980cdead9393a947c429afab06dff2833e0fcd8c6e916ad45aff09aa61b.png)
#### 17\.4.2.2. Rotate then shear[\#](https://intro.quantecon.org/eigen_I.html#rotate-then-shear "Link to this heading")
```
grid_composition_transform(B,A) # transformation BA
```
[](https://intro.quantecon.org/_images/c925d9e1e85a6a322bbec778a56c0f49b068ad58e929a63569c291d1b0e7714d.png)
It is evident that the transformation A B is not the same as the transformation B A.
## 17\.5. Iterating on a fixed map[\#](https://intro.quantecon.org/eigen_I.html#iterating-on-a-fixed-map "Link to this heading")
In economics (and especially in dynamic modeling), we are often interested in analyzing behavior where we repeatedly apply a fixed matrix.
For example, given a vector v and a matrix A, we are interested in studying the sequence
v , A v , A A v \= A 2 v , …
Let’s first see examples of a sequence of iterates ( A k v ) k ≥ 0 under different maps A.
```
def plot_series(A, v, n):
B = np.array([[1, -1],
[1, 0]])
fig, ax = plt.subplots()
ax.set(xlim=(-4, 4), ylim=(-4, 4))
ax.set_xticks([])
ax.set_yticks([])
for spine in ['left', 'bottom']:
ax.spines[spine].set_position('zero')
for spine in ['right', 'top']:
ax.spines[spine].set_color('none')
θ = np.linspace(0, 2 * np.pi, 150)
r = 2.5
x = r * np.cos(θ)
y = r * np.sin(θ)
x1 = x.reshape(1, -1)
y1 = y.reshape(1, -1)
xy = np.concatenate((x1, y1), axis=0)
ellipse = B @ xy
ax.plot(ellipse[0, :], ellipse[1, :], color='black',
linestyle=(0, (5, 10)), linewidth=0.5)
# Initialize holder for trajectories
colors = plt.cm.rainbow(np.linspace(0, 1, 20))
for i in range(n):
iteration = matrix_power(A, i) @ v
v1 = iteration[0]
v2 = iteration[1]
ax.scatter(v1, v2, color=colors[i])
if i == 0:
ax.text(v1+0.25, v2, f'$v$')
elif i == 1:
ax.text(v1+0.25, v2, f'$Av$')
elif 1 < i < 4:
ax.text(v1+0.25, v2, f'$A^{i}v$')
plt.show()
```
```
A = np.array([[sqrt(3) + 1, -2],
[1, sqrt(3) - 1]])
A = (1/(2*sqrt(2))) * A
v = (-3, -3)
n = 12
plot_series(A, v, n)
```
[](https://intro.quantecon.org/_images/d1ccd8bebb39aa9e536e913e0debb86cf5c0792a51abb85be4640953a15680c8.png)
Here with each iteration the vectors get shorter, i.e., move closer to the origin.
In this case, repeatedly multiplying a vector by A makes the vector “spiral in”.
```
B = np.array([[sqrt(3) + 1, -2],
[1, sqrt(3) - 1]])
B = (1/2) * B
v = (2.5, 0)
n = 12
plot_series(B, v, n)
```
[](https://intro.quantecon.org/_images/d7017d7b767cbbd55360a5dc5c4b2c9cb517f41b142884606b2e6d09abe070cb.png)
Here with each iteration vectors do not tend to get longer or shorter.
In this case, repeatedly multiplying a vector by A simply “rotates it around an ellipse”.
```
B = np.array([[sqrt(3) + 1, -2],
[1, sqrt(3) - 1]])
B = (1/sqrt(2)) * B
v = (-1, -0.25)
n = 6
plot_series(B, v, n)
```
[](https://intro.quantecon.org/_images/98508a529d57df518cbf11b072dfda23b1643217eb299fdb6448349b6d304a76.png)
Here with each iteration vectors tend to get longer, i.e., farther from the origin.
In this case, repeatedly multiplying a vector by A makes the vector “spiral out”.
We thus observe that the sequence ( A k v ) k ≥ 0 behaves differently depending on the map A itself.
We now discuss the property of A that determines this behavior.
## 17\.6. Eigenvalues[\#](https://intro.quantecon.org/eigen_I.html#eigenvalues "Link to this heading")
In this section we introduce the notions of eigenvalues and eigenvectors.
### 17\.6.1. Definitions[\#](https://intro.quantecon.org/eigen_I.html#definitions "Link to this heading")
Let A be an n × n square matrix.
If λ is scalar and v is a non-zero n\-vector such that
A v \= λ v .
Then we say that λ is an *eigenvalue* of A, and v is the corresponding *eigenvector*.
Thus, an eigenvector of A is a nonzero vector v such that when the map A is applied, v is merely scaled.
The next figure shows two eigenvectors (blue arrows) and their images under A (red arrows).
As expected, the image A v of each v is just a scaled version of the original
### 17\.6.2. Complex values[\#](https://intro.quantecon.org/eigen_I.html#complex-values "Link to this heading")
So far our definition of eigenvalues and eigenvectors seems straightforward.
There is one complication we haven’t mentioned yet:
When solving A v \= λ v,
- λ is allowed to be a complex number and
- v is allowed to be an n\-vector of complex numbers.
We will see some examples below.
### 17\.6.3. Some mathematical details[\#](https://intro.quantecon.org/eigen_I.html#some-mathematical-details "Link to this heading")
We note some mathematical details for more advanced readers.
(Other readers can skip to the next section.)
The eigenvalue equation is equivalent to ( A − λ I ) v \= 0.
This equation has a nonzero solution v only when the columns of A − λ I are linearly dependent.
This in turn is equivalent to stating the determinant is zero.
Hence, to find all eigenvalues, we can look for λ such that the determinant of A − λ I is zero.
This problem can be expressed as one of solving for the roots of a polynomial in λ of degree n.
This in turn implies the existence of n solutions in the complex plane, although some might be repeated.
### 17\.6.4. Facts[\#](https://intro.quantecon.org/eigen_I.html#facts "Link to this heading")
Some nice facts about the eigenvalues of a square matrix A are as follows:
1. the determinant of A equals the product of the eigenvalues
2. the trace of A (the sum of the elements on the principal diagonal) equals the sum of the eigenvalues
3. if A is symmetric, then all of its eigenvalues are real
4. if A is invertible and λ 1 , … , λ n are its eigenvalues, then the eigenvalues of A − 1 are 1 / λ 1 , … , 1 / λ n.
A corollary of the last statement is that a matrix is invertible if and only if all its eigenvalues are nonzero.
### 17\.6.5. Computation[\#](https://intro.quantecon.org/eigen_I.html#computation "Link to this heading")
Using NumPy, we can solve for the eigenvalues and eigenvectors of a matrix as follows
```
from numpy.linalg import eig
A = ((1, 2),
(2, 1))
A = np.array(A)
evals, evecs = eig(A)
evals # eigenvalues
```
```
array([ 3., -1.])
```
```
evecs # eigenvectors
```
```
array([[ 0.70710678, -0.70710678],
[ 0.70710678, 0.70710678]])
```
Note that the *columns* of `evecs` are the eigenvectors.
Since any scalar multiple of an eigenvector is an eigenvector with the same eigenvalue (which can be verified), the `eig` routine normalizes the length of each eigenvector to one.
The eigenvectors and eigenvalues of a map A determine how a vector v is transformed when we repeatedly multiply by A.
This is discussed further later.
## 17\.7. The Neumann Series Lemma[\#](https://intro.quantecon.org/eigen_I.html#the-neumann-series-lemma "Link to this heading")
In this section we present a famous result about series of matrices that has many applications in economics.
### 17\.7.1. Scalar series[\#](https://intro.quantecon.org/eigen_I.html#scalar-series "Link to this heading")
Here’s a fundamental result about series:
If a is a number and \| a \| \< 1, then
(17.1)[\#](https://intro.quantecon.org/eigen_I.html#equation-gp-sum "Link to this equation")
∑
k
\=
0
∞
a
k
\=
1
1
−
a
\=
(
1
−
a
)
−
1
For a one-dimensional linear equation x \= a x \+ b where x is unknown we can thus conclude that the solution x ∗ is given by:
x ∗ \= b 1 − a \= ∑ k \= 0 ∞ a k b
### 17\.7.2. Matrix series[\#](https://intro.quantecon.org/eigen_I.html#matrix-series "Link to this heading")
A generalization of this idea exists in the matrix setting.
Consider the system of equations x \= A x \+ b where A is an n × n square matrix and x and b are both column vectors in R n.
Using matrix algebra we can conclude that the solution to this system of equations will be given by:
What guarantees the existence of a unique vector x ∗ that satisfies [(17.2)](https://intro.quantecon.org/eigen_I.html#equation-neumann-eqn)?
The following is a fundamental result in functional analysis that generalizes [(17.1)](https://intro.quantecon.org/eigen_I.html#equation-gp-sum) to a multivariate case.
Theorem 17.1 (Neumann Series Lemma)
Let A be a square matrix and let A k be the k\-th power of A.
Let r ( A ) be the **spectral radius** of A, defined as max i \| λ i \|, where
- { λ i } i is the set of eigenvalues of A and
- \| λ i \| is the modulus of the complex number λ i
Neumann’s Theorem states the following: If r ( A ) \< 1, then I − A is invertible, and
( I − A ) − 1 \= ∑ k \= 0 ∞ A k
We can see the Neumann Series Lemma in action in the following example.
```
A = np.array([[0.4, 0.1],
[0.7, 0.2]])
evals, evecs = eig(A) # finding eigenvalues and eigenvectors
r = max(abs(λ) for λ in evals) # compute spectral radius
print(r)
```
```
0.5828427124746189
```
The spectral radius r ( A ) obtained is less than 1.
Thus, we can apply the Neumann Series Lemma to find ( I − A ) − 1.
```
I = np.identity(2) # 2 x 2 identity matrix
B = I - A
```
```
B_inverse = np.linalg.inv(B) # direct inverse method
```
```
A_sum = np.zeros((2, 2)) # power series sum of A
A_power = I
for i in range(50):
A_sum += A_power
A_power = A_power @ A
```
Let’s check equality between the sum and the inverse methods.
```
np.allclose(A_sum, B_inverse)
```
```
True
```
Although we truncate the infinite sum at k \= 50, both methods give us the same result which illustrates the result of the Neumann Series Lemma.
## 17\.8. Exercises[\#](https://intro.quantecon.org/eigen_I.html#exercises "Link to this heading")
Exercise 17.1
Power iteration is a method for finding the greatest absolute eigenvalue of a diagonalizable matrix.
The method starts with a random vector b 0 and repeatedly applies the matrix A to it
b k \+ 1 \= A b k ‖ A b k ‖
A thorough discussion of the method can be found [here](https://pythonnumericalmethods.berkeley.edu/notebooks/chapter15.02-The-Power-Method.html).
In this exercise, first implement the power iteration method and use it to find the greatest absolute eigenvalue and its corresponding eigenvector.
Then visualize the convergence.
Exercise 17.2
We have discussed the trajectory of the vector v after being transformed by A.
Consider the matrix A \= \[ 1 2 1 1 \] and the vector v \= \[ 2 − 2 \].
Try to compute the trajectory of v after being transformed by A for n \= 4 iterations and plot the result.
Exercise 17.3
[Previously](https://intro.quantecon.org/eigen_I.html#plot-series), we demonstrated the trajectory of the vector v after being transformed by A for three different matrices.
Use the visualization in the previous exercise to explain the trajectory of the vector v after being transformed by A for the three different matrices. |
| Shard | 82 (laksa) |
| Root Hash | 3160520800126211882 |
| Unparsed URL | org,quantecon!intro,/eigen_I.html s443 |