In [1]:
using Plots,LaTeXStrings
default(markersize=3,linewidth=1.5)
using LightGraphs,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.
In [2]:
g = watts_strogatz(200,4,0.06);
gplot(g)
Out[2]:
The adjacency matrix for this graph reveals the connections as mostly local (i.e., the nonzeros are near the diagonal).
In [3]:
A = adjacency_matrix(g,Float64);
spy(A,m=2,color=:blues,title="Adjacency matrix",leg=:none)
Out[3]:
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.
In [4]:
img = testimage("peppers")
Out[4]:
The details vary by image type, but for the most part an image is an array of color values.
In [5]:
@show size(img),eltype(img)
Out[5]:
The elements here have four values, for red, green, blue, and alpha (opacity). We can convert each of those "planes" into an ordinary matrix.
In [6]:
R = red.(img)
R[1:5,1:5]
Out[6]:
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).
In [7]:
G = Gray.(img)
Out[7]:
In [8]:
A = @. gray(Gray(img))
A[1:5,1:5]
Out[8]:
Finally, we can save an image locally for reloading later.
In [9]:
save("peppers.png",Gray.(img))
In [10]:
load("peppers.png")
Out[10]:
The eigvals command will return just the eigenvalues, as a vector.
In [11]:
A = pi*ones(2,2)
lambda = eigvals(A)
Out[11]:
If we also want the eigenvectors (returned as the matrix $V$), we use eigen.
In [12]:
lambda,V = eigen(A)
Out[12]:
We can check the fact that this is an EVD.
In [13]:
D = diagm(0=>lambda)
opnorm( A - V*D/V ) # "/V" is like "*inv(V)""
Out[13]:
Even if the matrix is not diagonalizable, eigen will run successfully, but the matrix ${V}$ will not be invertible.
In [14]:
lambda,V = eigen([1 1;0 1])
@show rank(V);
We will confirm the Bauer-Fike theorem on a triangular matrix. These tend to be far from normal.
In [15]:
n = 15;
lambda = 1:n
A = triu( ones(n)*lambda' );
A[1:5,1:5]
Out[15]:
The Bauer-Fike theorem provides an upper bound on the condition number of these eigenvalues.
In [16]:
lambda,V = eigen(A)
@show cond(V);
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.
In [17]:
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
Let's start with a known set of eigenvalues and an orthogonal eigenvector basis.
In [18]:
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.
In [19]:
Q,R = qr(A)
A = R*Q;
It turns out that this is a similarity transformation, so the eigenvalues are unchanged.
In [20]:
sort( eigvals(A) )
Out[20]:
What's remarkable is that if we repeat the transformation many times, the process converges to $D$.
In [21]:
for k = 1:40
Q,R = qr(A)
A = R*Q
end
A
Out[21]:
We verify some of the fundamental SVD properties using the built-in Julia command svd.
In [22]:
A = [i^j for i=1:5, j=0:3]
Out[22]:
In [23]:
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.
In [24]:
@show size(U),opnorm(U'*U - I);
In [25]:
@show size(V),opnorm(V'*V - I);
In [26]:
sigma
Out[26]:
In [27]:
@show opnorm(A),sigma[1];
In [28]:
@show cond(A), sigma[1]/sigma[end];
The following matrix is not hermitian.
In [29]:
A = [0 2; -2 0]
Out[29]:
It has an eigenvalue decomposition with a unitary matrix of eigenvectors, though, so it is normal.
In [30]:
lambda,V = eigen(A)
opnorm( V'*V - I )
Out[30]:
The eigenvalues are pure imaginary.
In [31]:
lambda
Out[31]:
The singular values are the complex magnitudes of the eigenvalues.
In [32]:
svdvals(A)
Out[32]:
We construct a real symmetric matrix with known eigenvalues by using the QR factorization to produce a random orthogonal set of eigenvectors.
In [33]:
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$.
In [34]:
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
We construct a symmetric matrix with a known EVD.
In [35]:
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.
In [36]:
R = x -> (x'*A*x)/(x'*x);
R(V[:,7])
Out[36]:
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.
In [37]:
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)
Out[37]:
We make an image from some text, then reload it as a matrix.
In [38]:
plot([],[],leg=:none,annotations=(0.5,0.5,text("Hello world",44,:center,:middle)),
grid=:none,frame=:none)
Out[38]:
In [39]:
savefig("hello.png")
img = load("hello.png")
A = @. Float64(Gray(img));
@show m,n = size(A);
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.
In [40]:
U,sigma,V = svd(A)
scatter(sigma,
title="Singular values",xaxis=(L"i"), yaxis=(:log10,L"\sigma_i"),leg=:none )
Out[40]:
In [41]:
r = findlast(@.sigma/sigma[1] > 10*eps())
Out[41]:
The rapid decrease suggests that we can get fairly good low-rank approximations.
In [42]:
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)
Out[42]:
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.
In [43]:
compression = 8*(m+n+1) / (m*n)
Out[43]:
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.
In [44]:
vars = load("voting.jld")
A = vars["A"]
m,n = size(A)
Out[44]:
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.
In [45]:
heatmap(A,color=:viridis,
title="Votes in 111th U.S. Senate",xlabel="bill",ylabel="senator")
Out[45]:
We use singular value "energy" to quantify the decay rate of the values.
In [46]:
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")
Out[46]:
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.
In [47]:
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")
Out[47]:
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.
In [48]:
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")
Out[48]: