17. Eigenvalues and Eigenvectors#

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 \times m$ matrix $A$ with an $m \times 1$ column vector $x$ to obtain an $n \times 1$ column vector $y$:

If we fix $A$ and consider different choices of $x$, we can understand $A$ as a map transforming $x$ to $Ax$.

Because $A$ is $n \times m$, it transforms $m$-vectors to $n$-vectors.

We can write this formally as $A \colon \mathbb{R}^m \rightarrow \mathbb{R}^n$.

You might argue that if $A$ is a function then we should write $A(x) = y$ rather than $Ax = 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 $\mathbb R^n$ to itself.

This means $A$ is an $n \times n$ matrix that maps (or “transforms”) a vector $x$ in $\mathbb{R}^n$ to a new vector $y=Ax$ also in $\mathbb{R}^n$.

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 $\theta$ and

• then scales it by some scalar $\gamma$ to obtain the image $y$ of $x$.

17.3.1. Scaling#

A matrix of the form

scales vectors across the x-axis by a factor $\alpha$ and along the y-axis by a factor $\beta$.

Here we illustrate a simple example where $\alpha = \beta = 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

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

is called a rotation matrix.

This matrix rotates vectors clockwise by an angle $\theta$.

θ = 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

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.1. Linear compositions#

Consider the two matrices

What will the output be when we try to obtain $ABx$ for some $2 \times 1$ vector $x$?

We can observe that applying the transformation $AB$ on the vector $x$ is the same as first applying $B$ on $x$ and then applying $A$ on the vector $Bx$.

Thus the matrix product $AB$ 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 \times m$ matrix $A$ with an $m \times k$ matrix $B$ the obtained matrix product is an $n \times k$ matrix $AB$.

Thus, if $A$ and $B$ are transformations such that $A \colon \mathbb{R}^m \to \mathbb{R}^n$ and $B \colon \mathbb{R}^k \to \mathbb{R}^m$, then $AB$ transforms $\mathbb{R}^k$ to $\mathbb{R}^n$.

Viewing matrix multiplication as composition of maps helps us understand why, under matrix multiplication, $AB$ is generally not equal to $BA$.

(After all, when we compose functions, the order usually matters.)

17.4.2. Examples#

Let $A$ be the $90^{\circ}$ clockwise rotation matrix given by $\begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix}$ and let $B$ be a shear matrix along the x-axis given by $\begin{bmatrix} 1 & 2 \\ 0 & 1 \end{bmatrix}$.

We will visualize how a grid of points changes when we apply the transformation $AB$ and then compare it with the transformation $BA$.

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 $AB$ is not the same as the transformation $BA$.

17.6.1. Definitions#

Let $A$ be an $n \times n$ square matrix.

If $\lambda$ is scalar and $v$ is a non-zero $n$-vector such that

Then we say that $\lambda$ 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 $Av$ 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()

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 $Av = \lambda v$,

• $\lambda$ 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 - \lambda I) v = 0$.

This equation has a nonzero solution $v$ only when the columns of $A - \lambda I$ are linearly dependent.

This in turn is equivalent to stating the determinant is zero.

Hence, to find all eigenvalues, we can look for $\lambda$ such that the determinant of $A - \lambda I$ is zero.

This problem can be expressed as one of solving for the roots of a polynomial in $\lambda$ 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:

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 $\lambda_1, \ldots, \lambda_n$ are its eigenvalues, then the eigenvalues of $A^{-1}$ are $1/\lambda_1, \ldots, 1/\lambda_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.1. Scalar series#

Here’s a fundamental result about series:

If $a$ is a number and $|a| < 1$, then

For a one-dimensional linear equation $x = ax + b$ where x is unknown we can thus conclude that the solution $x^{*}$ is given by:

17.7.2. Matrix series#

A generalization of this idea exists in the matrix setting.

Consider the system of equations $x = Ax + b$ where $A$ is an $n \times n$ square matrix and $x$ and $b$ are both column vectors in $\mathbb{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.

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.