# run Pkgs on first use
# import Pkg; Pkg.add("LightGraphs")
# import Pkg; Pkg.add("GraphPlot")
# import Pkg; Pkg.add("Images")
# import Pkg; Pkg.add("TestImages")
# import Pkg; Pkg.add("JLD")
# import Pkg; Pkg.add("Graphs") # replaces LightGraphs
using Plots,LaTeXStrings
default(markersize=3,linewidth=1.5)
#using LightGraphs,GraphPlot
using Graphs,GraphPlot
using Images,TestImages
using DataFrames,JLD
using LinearAlgebra
#include("FNC.jl");
Here is the adjacency matrix of a "small-world" network on 200 nodes. Each node is connected to 4 neighbors, and then some edges are randomly changed to distant connections.
g = watts_strogatz(200,4,0.06);
gplot(g)
The adjacency matrix for this graph reveals the connections as mostly local (i.e., the nonzeros are near the diagonal).
A = adjacency_matrix(g,Float64);
spy(A,m=2,color=:blues,title="Adjacency matrix",leg=:none)
We will use the Images
package for working with images. We also load here the TestImages
package for a large library of well-known standard images.
img = testimage("peppers")
The details vary by image type, but for the most part an image is an array of color values.
@show size(img),eltype(img)
(size(img), eltype(img)) = ((512, 512), RGB{N0f8})
((512, 512), RGB{N0f8})
The elements here have four values, for red, green, blue, and alpha (opacity). We can convert each of those "planes" into an ordinary matrix.
R = red.(img)
R[1:5,1:5]
5×5 Array{N0f8,2} with eltype N0f8: 0.396 0.549 0.592 0.627 0.631 0.482 0.749 0.745 0.729 0.714 0.494 0.725 0.725 0.722 0.702 0.482 0.733 0.729 0.718 0.729 0.498 0.733 0.714 0.757 0.722
The values above go from zero (no red) to one (full red). It may also be convenient to convert the image to grayscale, which has just one "layer" from zero (black) to one (white).
G = Gray.(img)
A = @. gray(Gray(img))
A[1:5,1:5]
5×5 Array{N0f8,2} with eltype N0f8: 0.118 0.227 0.243 0.251 0.255 0.145 0.467 0.447 0.451 0.451 0.149 0.475 0.451 0.427 0.439 0.145 0.424 0.475 0.427 0.447 0.149 0.451 0.486 0.463 0.431
Finally, we can save an image locally for reloading later.
save("peppers.png",Gray.(img))
load("peppers.png")
The eigvals
command will return just the eigenvalues, as a vector.
A = pi*ones(2,2)
lambda = eigvals(A)
2-element Vector{Float64}: 0.0 6.283185307179586
If we also want the eigenvectors (returned as the matrix $V$), we use eigen
.
lambda,V = eigen(A)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}} values: 2-element Vector{Float64}: 0.0 6.283185307179586 vectors: 2×2 Matrix{Float64}: -0.707107 0.707107 0.707107 0.707107
We can check the fact that this is an EVD.
D = diagm(0=>lambda)
opnorm( A - V*D/V ) # "/V" is like "*inv(V)""
8.881784197001252e-16
Even if the matrix is not diagonalizable, eigen
will run successfully, but the matrix ${V}$ will not be invertible.
lambda,V = eigen([1 1;0 1])
@show rank(V);
rank(V) = 1
We will confirm the Bauer-Fike theorem on a triangular matrix. These tend to be far from normal.
n = 15;
lambda = 1:n
A = triu( ones(n)*lambda' );
A[1:5,1:5]
5×5 Matrix{Float64}: 1.0 2.0 3.0 4.0 5.0 0.0 2.0 3.0 4.0 5.0 0.0 0.0 3.0 4.0 5.0 0.0 0.0 0.0 4.0 5.0 0.0 0.0 0.0 0.0 5.0
The Bauer-Fike theorem provides an upper bound on the condition number of these eigenvalues.
lambda,V = eigen(A)
@show cond(V);
cond(V) = 7.197767264538147e7
The theorem suggests that eigenvalue changes may be up to 7 orders of magnitude larger than a perturbation to the matrix. A few random experiments show that effects of nearly that size are not hard to observe.
for k = 1:3
E = randn(n,n); E = 1e-7*E/opnorm(E);
mu = eigvals(A+E)
@show max_change = norm( sort(mu)-lambda, Inf )
end
max_change = norm(sort(mu) - lambda, Inf) = 0.35474142721455415 max_change = norm(sort(mu) - lambda, Inf) = 0.1399620421565455 max_change = norm(sort(mu) - lambda, Inf) = 0.15959698059187488
Let's start with a known set of eigenvalues and an orthogonal eigenvector basis.
D = diagm(0=>[-6,-1,2,4,5])
V,R = qr(randn(5,5))
A = V*D*V'; # note that V' = inv(V)
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.
sort( eigvals(A) )
5-element Vector{Float64}: -6.0000000000000036 -0.9999999999999999 1.9999999999999991 3.9999999999999996 5.000000000000002
What's remarkable is that if we repeat the transformation many times, the process converges to $D$.
for k = 1:40
Q,R = qr(A)
A = R*Q
end
A
5×5 Matrix{Float64}: -5.99998 0.0129051 -1.21121e-6 3.54782e-16 -1.86016e-15 0.0129051 4.99998 -0.000131068 -2.16625e-15 6.69248e-16 -1.21121e-6 -0.000131068 4.0 -8.8957e-13 -5.33649e-17 -6.2024e-19 -7.22661e-18 -8.90899e-13 2.0 1.2894e-12 1.24676e-31 6.16686e-29 -6.5909e-25 1.28981e-12 -1.0
We verify some of the fundamental SVD properties using the built-in Julia command svd
.
A = [i^j for i=1:5, j=0:3]
5×4 Matrix{Int64}: 1 1 1 1 1 2 4 8 1 3 9 27 1 4 16 64 1 5 25 125
U,sigma,V = svd(A);
Note that while the "full" SVD has a square $U$, the "thin" form is the default. Here the columns are orthonormal even though $U$ is not square.
@show size(U),opnorm(U'*U - I);
(size(U), opnorm(U' * U - I)) = ((5, 4), 1.3853853390917798e-15)
@show size(V),opnorm(V'*V - I);
(size(V), opnorm(V' * V - I)) = ((4, 4), 7.154803690039982e-16)
sigma
4-element Vector{Float64}: 146.69715365883005 5.738569780953701 0.9998486640841024 0.1192808268524193
@show opnorm(A),sigma[1];
(opnorm(A), sigma[1]) = (146.69715365883005, 146.69715365883005)
@show cond(A), sigma[1]/sigma[end];
(cond(A), sigma[1] / sigma[end]) = (1229.8468876337665, 1229.8468876337663)
The following matrix is not hermitian.
A = [0 2; -2 0]
2×2 Matrix{Int64}: 0 2 -2 0
It has an eigenvalue decomposition with a unitary matrix of eigenvectors, though, so it is normal.
lambda,V = eigen(A)
opnorm( V'*V - I )
2.220446049250313e-16
The eigenvalues are pure imaginary.
lambda
2-element Vector{ComplexF64}: 0.0 - 2.0000000000000004im 0.0 + 2.0000000000000004im
The singular values are the complex magnitudes of the eigenvalues.
svdvals(A)
2-element Vector{Float64}: 2.0 2.0
We construct a real symmetric matrix with known eigenvalues by using the QR factorization to produce a random orthogonal set of eigenvectors.
n = 30;
lambda = 1:n
D = diagm(0=>lambda)
V,R = qr(randn(n,n)) # get a random orthogonal V
A = V*D*V';
The condition number of these eigenvalues is one. Thus the effect on them is bounded by the norm of the perturbation to $A$.
for k = 1:3
E = randn(n,n); E = 1e-4*E/opnorm(E);
mu = sort(eigvals(A+E))
@show max_change = norm(mu-lambda,Inf)
end
max_change = norm(mu - lambda, Inf) = 2.575407583194078e-5 max_change = norm(mu - lambda, Inf) = 1.8812706223503284e-5 max_change = norm(mu - lambda, Inf) = 1.9617580392150558e-5
We construct a symmetric matrix with a known EVD.
n = 20;
lambda = 1:n
D = diagm(0=>lambda)
V,R = qr(randn(n,n)) # get a random orthogonal V
A = V*D*V';
The Rayleigh quotient of an eigenvector is its eigenvalue.
R = x -> (x'*A*x)/(x'*x);
R(V[:,7])
7.0
The Rayleigh quotient's value is much closer to an eigenvalue than its input is to an eigenvector. In this experiment, each additional digit of accuracy in the eigenvector estimate gives two more digits to the eigenvalue estimate.
delta = @. 1 ./10^(1:4)
quotient = zeros(size(delta))
for (k,delta) = enumerate(delta)
e = randn(n); e = delta*e/norm(e);
x = V[:,7] + e
quotient[k] = R(x)
end
DataFrame(perturbation=delta,RQminus7=quotient.-7)
4 rows × 2 columns
perturbation | RQminus7 | |
---|---|---|
Float64 | Float64 | |
1 | 0.1 | 0.0242551 |
2 | 0.01 | 0.00032191 |
3 | 0.001 | 4.02396e-6 |
4 | 0.0001 | 1.21167e-8 |
We make an image from some text, then reload it as a matrix.
plot([],[],leg=:none,annotations=(0.5,0.5,text("Hello world",44,:center,:middle)),
grid=:none,frame=:none)
┌ Warning: Keyword argument letter not supported with Plots.GRBackend(). Choose from: Set([:top_margin, :group, :inset_subplots, :background_color, :ytickfontsize, :yforeground_color_text, :yguidefontcolor, :tickfontfamily, :show_empty_bins, :seriesalpha, :seriescolor, :ztick_direction, :xgrid, :ygridalpha, :zlims, :xtick_direction, :colorbar, :legend_font_family, :zflip, :ticks, :linealpha, :overwrite_figure, :arrow, :xguidefonthalign, :normalize, :linestyle, :xtickfontvalign, :xflip, :zgrid, :fillcolor, :ygrid, :bar_width, :colorbar_scale, :background_color_inside, :zguidefonthalign, :bins, :zguide, :zforeground_color_text, :legend_font_valign, :yscale, :legend_font_color, :weights, :xgridalpha, :ygridstyle, :clims, :xtickfontcolor, :fill_z, :xguide, :markershape, :background_color_subplot, :ztickfontfamily, :fillalpha, :markerstrokewidth, :tick_direction, :xguidefontvalign, :xguidefontfamily, :gridlinewidth, :foreground_color_subplot, :xgridlinewidth, :yguidefontsize, :foreground_color, :foreground_color_text, :titlefonthalign, :yerror, :x, :xtickfonthalign, :zgridlinewidth, :ytickfontrotation, :discrete_values, :ytick_direction, :grid, :xguidefontrotation, :ribbon, :xguidefontsize, :tickfontrotation, :xforeground_color_axis, :xdiscrete_values, :background_color_outside, :titlefontcolor, :xgridstyle, :line_z, :size, :orientation, :gridstyle, :projection, :markersize, :legend_foreground_color, :camera, :zguidefontrotation, :ydiscrete_values, :xforeground_color_grid, :seriestype, :yflip, :quiver, :zticks, :markerstrokecolor, :ztickfontrotation, :ztickfonthalign, :fillrange, :ztickfontvalign, :xlims, :xforeground_color_border, :markercolor, :xtickfontsize, :ylink, :levels, :color_palette, :connections, :yforeground_color_grid, :lims, :zgridstyle, :foreground_color_border, :zguidefontvalign, :xscale, :marker_z, :markerstrokealpha, :left_margin, :markeralpha, :legend_font_halign, :annotations, :window_title, :tickfontvalign, :foreground_color_axis, :zguidefontcolor, :ygridlinewidth, :zlink, :zscale, :smooth, :yguidefontrotation, :xticks, :guidefontsize, :zguidefontsize, :y, :margin, :ytickfontcolor, :zdiscrete_values, :tickfonthalign, :bottom_margin, :yforeground_color_border, :zguidefontfamily, :framestyle, :yguidefontvalign, :yguidefonthalign, :zerror, :zgridalpha, :ztickfontcolor, :scale, :legend_position, :linecolor, :html_output_format, :legend_title, :zforeground_color_border, :legend_font_pointsize, :title, :tickfontcolor, :subplot_index, :flip, :titlefontrotation, :legend_background_color, :tickfontsize, :titlefontvalign, :z, :yforeground_color_axis, :foreground_color_grid, :xtickfontrotation, :linewidth, :ztickfontsize, :gridalpha, :xerror, :guidefontfamily, :ylims, :contour_labels, :xguidefontcolor, :primary, :xtickfontfamily, :ytickfontvalign, :guidefonthalign, :ytickfontfamily, :aspect_ratio, :xforeground_color_text, :show, :link, :colorbar_title, :guidefontrotation, :subplot, :label, :ytickfonthalign, :guide, :guidefontcolor, :yguide, :titlefontsize, :titlefontfamily, :guidefontvalign, :zforeground_color_axis, :zforeground_color_grid, :layout, :legend_font_rotation, :colorbar_entry, :yguidefontfamily, :polar, :right_margin, :xlink, :series_annotations, :yticks]) └ @ Plots /home/hermanr/.julia/packages/Plots/FI0vT/src/args.jl:1607
savefig("hello.png")
img = load("hello.png")
A = @. Float64(Gray(img));
@show m,n = size(A);
(m, n) = size(A) = (400, 600)
Next we show that the singular values decrease exponentially, until they reach zero (more precisely, are about $\sigma_1 \varepsilon_\text{mach}$). For all numerical purposes, this determines the rank of the matrix.
U,sigma,V = svd(A)
scatter(sigma,
title="Singular values",xaxis=(L"i"), yaxis=(:log10,L"\sigma_i"),leg=:none )
r = findlast(@.sigma/sigma[1] > 10*eps())
44
The rapid decrease suggests that we can get fairly good low-rank approximations.
Ak = [ U[:,1:k]*diagm(0=>sigma[1:k])*V[:,1:k]' for k=2*(1:4) ]
reshape( [ @.Gray(Ak[i]) for i=1:4 ],2,2)
2×2 Matrix{Matrix{Gray{Float64}}}: [Gray{Float64}(0.999468) Gray{Float64}(0.999468) … Gray{Float64}(0.999468) Gray{Float64}(0.999468); Gray{Float64}(0.999468) Gray{Float64}(0.999468) … Gray{Float64}(0.999468) Gray{Float64}(0.999468); … ; Gray{Float64}(0.999468) Gray{Float64}(0.999468) … Gray{Float64}(0.999468) Gray{Float64}(0.999468); Gray{Float64}(0.999468) Gray{Float64}(0.999468) … Gray{Float64}(0.999468) Gray{Float64}(0.999468)] … [Gray{Float64}(0.999987) Gray{Float64}(0.999987) … Gray{Float64}(0.999987) Gray{Float64}(0.999987); Gray{Float64}(0.999987) Gray{Float64}(0.999987) … Gray{Float64}(0.999987) Gray{Float64}(0.999987); … ; Gray{Float64}(0.999987) Gray{Float64}(0.999987) … Gray{Float64}(0.999987) Gray{Float64}(0.999987); Gray{Float64}(0.999987) Gray{Float64}(0.999987) … Gray{Float64}(0.999987) Gray{Float64}(0.999987)] [Gray{Float64}(0.999322) Gray{Float64}(0.999322) … Gray{Float64}(0.999322) Gray{Float64}(0.999322); Gray{Float64}(0.999322) Gray{Float64}(0.999322) … Gray{Float64}(0.999322) Gray{Float64}(0.999322); … ; Gray{Float64}(0.999322) Gray{Float64}(0.999322) … Gray{Float64}(0.999322) Gray{Float64}(0.999322); Gray{Float64}(0.999322) Gray{Float64}(0.999322) … Gray{Float64}(0.999322) Gray{Float64}(0.999322)] [Gray{Float64}(0.99999) Gray{Float64}(0.99999) … Gray{Float64}(0.99999) Gray{Float64}(0.99999); Gray{Float64}(0.99999) Gray{Float64}(0.99999) … Gray{Float64}(0.99999) Gray{Float64}(0.99999); … ; Gray{Float64}(0.99999) Gray{Float64}(0.99999) … Gray{Float64}(0.99999) Gray{Float64}(0.99999); Gray{Float64}(0.99999) Gray{Float64}(0.99999) … Gray{Float64}(0.99999) Gray{Float64}(0.99999)]
Consider how little data is needed to reconstruct these images. For rank 8, for instance, we have 8 left and right singular vectors plus 8 singular values, for a compression ratio of better than 25:1.
compression = 8*(m+n+1) / (m*n)
0.03336666666666667
This matrix describes the votes on bills in the 111th session of the United States Senate. (The data set was obtained from voteview.com.) Each row is one senator and each column is a vote item.
vars = load("voting.jld")
A = vars["A"]
m,n = size(A)
(112, 696)
If we visualize the votes (white is "yea," black is "nay," and gray is anything else), we can see great similarity between many rows, reflecting party unity.
heatmap(A,color=:viridis,
title="Votes in 111th U.S. Senate",xlabel="bill",ylabel="senator")
We use singular value "energy" to quantify the decay rate of the values.
U,sigma,V = svd(A)
tau = cumsum(sigma.^2) / sum(sigma.^2)
scatter(tau[1:16],label="",
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 $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];
Rep = vec(vars["Rep"]); Dem = vec(vars["Dem"]); Ind = vec(vars["Ind"]);
scatter(x1[Dem],x2[Dem],color=:blue,label="D",
xaxis=("partisanship"),yaxis=("bipartisanship"),title="111th US Senate in 2D" )
scatter!(x1[Rep],x2[Rep],color=:red,label="R")
scatter!(x1[Ind],x2[Ind],color=:yellow,label="I")