In [1]:
# Pkg.add("MatrixMarket")
using MatrixMarket
In [2]:
varinfo(MatrixMarket)
Out[2]:
In [3]:
#############################################################
# Download and parse every file from the NIST Matrix Market #
#############################################################
#Convenience function to emulate the behavior of gunzip
using GZip
function gunzip(fname)
destname, ext = splitext(fname)
if ext != ".gz"
error("gunzip: $fname: unknown suffix -- ignored")
end
open(destname, "w") do f
GZip.open(fname) do g
write(f, string(g))
end
end
destname
end
#Download and parse master list of matrices
if !isfile("matrices.html")
download("http://math.nist.gov/MatrixMarket/matrices.html", "matrices.html")
end
matrixmarketdata = Any[]
open("matrices.html") do f
for line in readlines(f)
if occursin("""<A HREF="/MatrixMarket/data/""",line)
collectionname, setname, matrixname = split(split(line, '"')[2], '/')[4:6]
matrixname = split(matrixname, '.')[1]
push!(matrixmarketdata, (collectionname, setname, matrixname) )
end
end
end
#Download one matrix at random
n = rand(1:length(matrixmarketdata))
for (collectionname, setname, matrixname) in matrixmarketdata[n:n]
fn = string(collectionname, '_', setname, '_', matrixname)
mtxfname = string(fn, ".mtx")
if !isfile(mtxfname)
url = "ftp://math.nist.gov/pub/MatrixMarket2/$collectionname/$setname/$matrixname.mtx.gz"
gzfname = string(fn, ".mtx.gz")
try
download(url, gzfname)
catch
continue
end
gunzip(gzfname)
end
end
In [4]:
readdir()
Out[4]:
In [5]:
A=mmread("Harwell-Boeing_bcsstruc1_bcsstk11.mtx")
Out[5]:
In [6]:
size(A)
Out[6]:
In [7]:
using LinearAlgebra
In [8]:
issymmetric(A)
Out[8]:
In [9]:
cond(Matrix(A))
Out[9]:
In [10]:
using Plots
using SparseArrays
In [12]:
# Plot of a small random matrix, for illustration
B=rand(10,10)
spy(sparse(B))
heatmap(B)
Out[12]:
In [13]:
# Now our matrix - Wait!
heatmap(A)
Out[13]:
In [14]:
using MatrixDepot
In [15]:
varinfo(MatrixDepot)
Out[15]:
In [16]:
mdinfo()
Out[16]:
In [17]:
listnames(:eigen)
Out[17]:
In [18]:
mdlist(:sparse & :symmetric)
Out[18]:
In [19]:
A=matrixdepot("fiedler",100)
Out[19]:
In [20]:
cond(A)
Out[20]:
In [25]:
include("ModuleB.jl")
Out[25]:
In [26]:
ModuleB.myPowerMethod(A,1e-10)
Out[26]:
In [27]:
eigvals(A)
Out[27]:
In [28]:
mdlist(:illcond & :symmetric)
Out[28]:
google (and others)
Some programs:
You can also try function eigs() or the Power method.
We will try example from the Moler's paper.
In [25]:
i = vec([ 2 6 3 4 4 5 6 1 1])
j = vec([ 1 1 2 2 3 3 3 4 6])
Out[25]:
In [26]:
G=sparse(i,j,1.0)
Out[26]:
In [27]:
typeof(G)
Out[27]:
In [28]:
full(G)
Out[28]:
In [29]:
c=sum(G,1)
n=size(G,1)
for j=1:n
if c[j]>0
G[:,j]=G[:,j]/c[j]
end
end
In [30]:
full(G)
Out[30]:
In [31]:
p=0.85
δ=(1-p)/n
Out[31]:
In [32]:
z = ((1-p)*(c.!=0) + (c.==0))/n
Out[32]:
In [33]:
A=p*G+ones(n)*z
Out[33]:
In [34]:
sum(A,1)
Out[34]:
Let us start a random walk from the vector $x_0=\begin{bmatrix} 1/n \\ 1/n \\ \vdots \\ 1/n \end{bmatrix}$.
The subsequent vectors are
$$ x_1=A\cdot x_0 \\ x_2=A\cdot x_1 \\ x_3=A\cdot x_2\\ \vdots $$When the vector stabilizes:
$$ A\cdot x\approx x, $$then $x[i]$ is the rank of the page $i$.
In [35]:
function myPageRank(G::SparseMatrixCSC{Float64,Int64},steps::Int)
p=0.85
c=sum(G,1)/p
n=size(G,1)
for i=1:n
G.nzval[G.colptr[i]:G.colptr[i+1]-1]./=c[i]
end
e=ones(n)
x=e/n
z = vec(((1-p)*(c.!=0) + (c.==0))/n)
for j=1:steps
x=G*x+(z⋅x)
end
x/norm(x,1)
end
Out[35]:
In [36]:
fieldnames(G)
Out[36]:
In [37]:
G
Out[37]:
In [38]:
G.colptr
Out[38]:
In [39]:
G.nzval
Out[39]:
In [40]:
# Starting vector
x=ones(n)/n
Out[40]:
In [42]:
myPageRank(G,15)
Out[42]:
This is somewhat larger test problem.
In [43]:
W=readdlm("web-Stanford.txt",Int)
Out[43]:
In [44]:
?sparse;
In [45]:
S=sparse(W[:,2],W[:,1],1.0)
Out[45]:
In [46]:
@time x100=myPageRank(S,100);
In [47]:
x101=myPageRank(S,101);
In [48]:
maximum(abs,(x101-x100)./x101)
Out[48]:
In [49]:
# Ranks
sort(x100,rev=true)
Out[49]:
In [50]:
# Pages
sortperm(x100,rev=true)
Out[50]:
In [ ]: