Chapter 7

Here we create an adjacency matrix for a graph on four nodes.

A = [0 1 0 0; 1 0 0 0; 1 1 0 1; 0 1 1 0]

The graphplot function makes a visual representation of this graph.

using Plots, GraphRecipes
graphplot(A, names=1:4, markersize=0.2, arrow=6)

Since this adjacency matrix is not symmetric, the edges are all directed, as indicated by the arrows. Here are the counts of all walks of length 3 in the graph:

A^3

If the adjacency matrix is symmetric, the result is an undirected graph: all edges connect in both directions.

A = [0 1 1 0; 1 0 0 1; 1 0 0 0; 0 1 0 0]
graphplot(A, names=1:4, markersize=0.2)

The Images package has many functions for image manipulation, and TestImages has some standard images to play with.

using Images, TestImages
img = testimage("mandrill")

The variable img is a matrix.

size(img)

However, its entries are colors, not numbers.

img[100, 10]

You can use eltype to find out the type of the elements of any array.

eltype(img)

It’s possible to extract matrices of red, green, and blue intensities, scaled from 0 to 1.

R, G, B = red.(img), green.(img), blue.(img);
@show minB, maxB = extrema(B);

(minB, maxB) = extrema(B) = (0.0N0f8, 1.0N0f8)

Or we can convert the pixels to gray, each pixel again scaled from 0 to 1.

Gray.(img)

In order to do our usual operations, we need to tell Julia that we want to interpret the elements of the image matrix as floating-point values.

A = Float64.(Gray.(img))
A[1:4, 1:5]

We can use Gray to reinterpret a matrix of floating-point values as grayscale pixels.

Gray.(reverse(A, dims=1))

The eigvals function returns a vector of the eigenvalues of a matrix.

A = π * ones(2, 2)

λ = eigvals(A)

If you want the eigenvectors as well, use eigen.

λ, V = eigen(A)

norm(A * V[:, 2] - λ[2] * V[:, 2])

Both functions allow you to sort the eigenvalues by specified criteria.

A = diagm(-2.3:1.7)
@show eigvals(A, sortby=real);
@show eigvals(A, sortby=abs);

eigvals(A, sortby = real) = [-2.3, -1.3, -0.3, 0.7, 1.7]

eigvals(A, sortby = abs) = 

[-0.3, 0.7, -1.3, 1.7, -2.3]

If the matrix is not diagonalizable, no message is given, but V will be singular. The robust way to detect that circumstance is via $\kappa(\mathbf{V})$.

A = [-1 1; 0 -1]
λ, V = eigen(A)

cond(V)

Even in the nondiagonalizable case, $\mathbf{A}\mathbf{V} = \mathbf{V}\mathbf{D}$ holds.

opnorm(A * V - V * diagm(λ))

We first define a hermitian matrix. Note that the ' operation is the adjoint and includes complex conjugation.

n = 7
A = randn(n, n) + 1im * randn(n, n)
A = (A + A') / 2

We confirm that the matrix $\mathbf{A}$ is normal by checking that $\kappa(\mathbf{V}) = 1$ (to within roundoff).

λ, V = eigen(A)
@show cond(V);

cond(V) = 1.0000000000000007

Now we perturb $\mathbf{A}$ and measure the effect on the eigenvalues. The Bauer–Fike theorem uses absolute differences, not relative ones.

ΔA = 1e-8 * normalize(randn(n, n) + 1im * randn(n, n))
λ̃ = eigvals(A + ΔA)
dist = minimum([abs(x - y) for x in λ̃, y in λ], dims=2)

As promised, the perturbations in the eigenvalues do not exceed the normwise perturbation to the original matrix.

Now we see what happens for a triangular matrix.

n = 20
x = 1:n
A = triu(x * ones(n)')
A[1:5, 1:5]

This matrix is not especially close to normal.

λ, V = eigen(A)
@show cond(V);

cond(V) = 6.149906664929389e9

As a result, the eigenvalues can change by a good deal more.

ΔA = 1e-8 * normalize(randn(n, n) + 1im * randn(n, n))
λ̃ = eigvals(A + ΔA)
dist = minimum([abs(x - y) for x in λ̃, y in λ], dims=2)
BF_bound = cond(V) * norm(ΔA)
@show maximum(dist), BF_bound;

(maximum(dist), BF_bound) = (0.494851347444697, 61.499066649293894)

If we plot the eigenvalues of many perturbations, we get a cloud of points that roughly represents all the possible eigenvalues when representing this matrix with single-precision accuracy.

using Plots
plt = scatter(λ, zeros(n), aspect_ratio=1)
for _ in 1:200
    ΔA = eps(Float32) * normalize(randn(n, n) + 1im * randn(n, n))
    λ̃ = eigvals(A + ΔA)
    scatter!(real(λ̃), imag(λ̃), m=1, color=:black)
end
plt

The plot shows that some eigenvalues are much more affected than others. This situation is not unusual, but it is not explained by the Bauer–Fike theorem.

Let’s start with a known set of eigenvalues and an orthogonal eigenvector basis.

D = diagm([-6, -1, 2, 4, 5])
V, R = qr(randn(5, 5))    # V is unitary
A = V * D * V'

eigvals(A)

Now we will take the QR factorization and just reverse the factors.

Q, R = qr(A)
A = R * Q;

It turns out that this is a similarity transformation, so the eigenvalues are unchanged.

eigvals(A)

What’s remarkable, and not elementary, is that if we repeat this transformation many times, the resulting matrix converges to $\mathbf{D}$.

for k in 1:40
    Q, R = qr(A)
    A = R * Q
end
A

We verify some of the fundamental SVD properties using standard Julia functions from LinearAlgebra.

A = [i^j for i = 1:5, j = 0:3]

To get only the singular values, use svdvals.

σ = svdvals(A)

Here is verification of the connections between the singular values, norm, and condition number.

@show opnorm(A, 2);
@show σ[1];

opnorm(A, 2) = 146.69715365883005

σ[1] = 146.69715365883005

@show cond(A, 2);
@show σ[1] / σ[end];

cond(A, 2) = 1229.846887633767

σ[1] / σ[end] = 1229.846887633767

To get singular vectors as well, use svd. The thin form of the factorization is the default.

U, σ, V = svd(A);
@show size(U);
@show size(V);

size(U) = (5, 4)

size(V) = (4, 4)

We verify the orthogonality of the singular vectors as follows:

@show opnorm(U' * U - I);
@show opnorm(V' * V - I);

opnorm(U' * U - I) = 4.947764483928556e-16

opnorm(V' * V - I) = 7.570241024277813e-16

We will use a symmetric matrix with a known EVD and eigenvalues equal to the integers from 1 to 20.

n = 20;
λ = 1:n
D = diagm(λ)
V, _ = qr(randn(n, n))   # get a random orthogonal V
A = V * D * V';

The Rayleigh quotient is a scalar-valued function of a vector.

R = x -> (x' * A * x) / (x' * x);

The Rayleigh quotient evaluated at an eigenvector gives the corresponding eigenvalue.

R(V[:, 7])

If the input to he Rayleigh quotient is within a small δ of an eigenvector, its output is within $O(\delta^2)$ of the corresponding eigenvalue. In this experiment, we observe that each additional digit of accuracy in an approximate eigenvector gives two more digits to the eigenvalue estimate coming from the Rayleigh quotient.

δ = @. 1 ./ 10^(1:5)
eval_diff = zeros(size(δ))
for (k, delta) in enumerate(δ)
    e = randn(n)
    e = delta * e / norm(e)
    x = V[:, 7] + e
    eval_diff[k] = R(x) - 7
end
labels = ["perturbation δ", "δ²", "R(x) - λ"]
@pt :header=labels [δ δ .^ 2 eval_diff]

We make an image from some text, then reload it as a matrix.

using Plots, Images
plot(annotations=(0.5, 0.5, text("Hello world", 44, :center, :center)),
    grid=:none, frame=:none, size=(400, 150))
savefig("hello.png")
img = load("hello.png")
A = @. Float64(Gray(img))
Gray.(A)

Next we show that the singular values decrease until they reach zero (more precisely, until they are about $\epsilon_\text{mach}$ times the norm of the matrix) at around $k=45$.

U, σ, V = svd(A)
scatter(σ;
    xaxis=(L"i"),  yaxis=(:log10, L"\sigma_i"),
    title="Singular values")

The rapid decrease suggests that we can get fairly good low-rank approximations.

plt = plot(layout=(2, 2), frame=:none, aspect_ratio=1, titlefontsize=10)
for i in 1:4
    k = 3i
    Ak = U[:, 1:k] * diagm(σ[1:k]) * V[:, 1:k]'
    plot!(Gray.(Ak), subplot=i, title="rank = $k")
end
plt

Consider how little data is needed to reconstruct these images. For rank-9, for instance, we have 9 left and right singular vectors plus 9 singular values, for a compression ratio of better than 12:1.

m, n = size(A)
compression = m * n / (9 * (m + n + 1))

This matrix describes the votes on bills in the 111th session of the United States Senate. (The data set was obtained from [https://voteview.com].) Each row is one senator, and each column is a vote item.

using JLD2
@load "voting.jld2" A;

If we visualize the votes (yellow is “yea,” blue is “nay”), we can see great similarity between many rows, reflecting party unity.

heatmap(A;
    color=:viridis,  xlabel="bill",  ylabel="senator",
    title="Votes in 111th U.S. Senate")

We use (7.5.4) to quantify the decay rate of the values.

U, σ, V = svd(A)
τ = cumsum(σ .^ 2) / sum(σ .^ 2)
scatter(τ[1:16];
    xaxis=("k"),  yaxis=(L"\tau_k"),
    title="Fraction of singular value energy")

The first and second singular triples contain about 58% and 17%, respectively, of the energy of the matrix. All others have far less effect, suggesting that the information is primarily two-dimensional. The first left and right singular vectors also contain interesting structure.

scatter(U[:, 1], label="", layout=(1, 2),
    xlabel="senator",  title="left singular vector")
scatter!(V[:, 1], label="", subplot=2,
    xlabel="bill",  title="right singular vector")

Both vectors have values greatly clustered near $\pm C$ for a constant $C$. These can be roughly interpreted as how partisan a particular senator or bill was, and for which political party. Projecting the senators’ vectors into the first two $\mathbf{V}$-coordinates gives a particularly nice way to reduce them to two dimensions. Political scientists label these dimensions partisanship and bipartisanship. Here we color them by actual party affiliation (also given in the data file): red for Republican, blue for Democrat, and yellow for independent.

x1 = A * V[:, 1];
x2 = A * V[:, 2];

@load "voting.jld2" Rep Dem Ind
Rep = vec(Rep);
Dem = vec(Dem);
Ind = vec(Ind);
scatter(x1[Dem], x2[Dem];
    color=:blue,  label="D",
    xaxis=("partisanship"),  yaxis=("bipartisanship"), 
    title="111th US Senate by voting record")
scatter!(x1[Rep], x2[Rep], color=:red, label="R")
scatter!(x1[Ind], x2[Ind], color=:yellow, label="I")