This is a demo/instructions for solving equations in Laplacian and SDD matrices. The sections are:
In [1]:
using Laplacians
We first generate a SDD, Positive Definite, system, and solve it using a direct solver. This uses the amd ordering, and is very fast.
In [2]:
a = grid2(5)
la = lap(a)
la[1,1] = la[1,1] + 1
F = cholfact(la)
Out[2]:
We can now use F to solve systems in this matrix, la. It is a complex structure that encodes a cholesky factorization, but it is much more than that.
In [3]:
n = size(a)[1]
b = randn(n)
x = F \ b
norm(la*x-b)
Out[3]:
Let's poke around to see what F has inside it.
In [4]:
F
Out[4]:
In [5]:
super(typeof(F))
Out[5]:
In [6]:
isa(F,Factorization)
Out[6]:
In [7]:
fieldnames(F)
Out[7]:
Let's see how long this will take on biggish grids, and on random regular graphs (which will have a lot of fill)
In [8]:
n = 500;
a = grid2(n)
la = lap(a)
la[1,1] = la[1,1] + 1
@time F = cholfact(la);
For a 500-by-500 grid, it took 1.5 seconds. We will now see that to use the solver, it takes 0.4 seconds.
In [9]:
N = size(la)[1]
b = randn(N)
@time x = F \ b
norm(la*x-b)
Out[9]:
For a random regular graph, we hit 1.5 seconds at around 20k vertices.
In [10]:
N = 20000;
a = randRegular(N,3)
la = lap(a)
la[1,1] = la[1,1] + 1
@time F = cholfact(la);
The time required for the solve is then around 0.05 seconds.
In [11]:
b = randn(N)
@time x = F \ b
norm(la*x-b)
Out[11]:
In [12]:
n = 500;
a = grid2(n)
la = lap(a)
la[1,1] = la[1,1] + 1
N = size(la)[1]
b = randn(N)
@time x = la \ b
norm(la*x-b)
Out[12]:
Just using \ appears to be right, so it is probably using cholfact.
We solve Laplacian systems by solving a system in the induced submatrix. Here are the steps, which I will then put into a wrapper function. It works by solving in a submatrix, like this:
In [13]:
la = lap(grid2(500))
N = size(la)[1]
lasub = la[1:(N-1),1:(N-1)]
Fsub = cholfact(lasub);
In [14]:
b = randn(N)
b = b - mean(b)
bs = b[1:(N-1)]
xs = Fsub \ bs;
x = [xs;0]
x = x - mean(x)
norm(la*x-b)
Out[14]:
The following wraps a solver for SDD systems into a solver for Laplacian systems. We really need to work on the types of solver, and actually for everything else inside.
Let's see this work.
In [15]:
la = lap(a);
f = lapWrapSolver(cholfact,la)
b = randn(size(a)[1]); b = b - mean(b);
norm(la*f(b) - b)
Out[15]:
We now make two more versions: one that just takes the solver, and one that takes b as well.
In [16]:
lapChol2 = lapWrapSolver(cholfact)
Out[16]:
In [17]:
f = lapChol(la)
Out[17]:
In [18]:
norm(la*f(b) - b)
Out[18]:
In [19]:
x = lapWrapSolver(cholfact,la,b)
norm(la*x - b)
Out[19]:
In [ ]:
I really like the fact that Julia lets me type the following. It still needs reasonable types.
Here are examples of how to solve systems using the Conjugate Gradient.
In [20]:
n = 50
a = randn(n,n); a = a * a';
b = randn(n)
x = cg(a,b,maxits=100)
norm(a*x - b)
Out[20]:
In [21]:
bbig = convert(Array{BigFloat,1},b);
xbig = cg(a,bbig,maxits=100)
norm(a*xbig - bbig)
Out[21]:
In [22]:
la = lap(grid2(200))
n = size(la)[1]
b = randn(n)
b = b - mean(b);
In [23]:
@time x = cg(la,b,maxits=1000)
norm(la*x-b)
Out[23]:
In [24]:
a = mapweight(grid2(200),x->1/(rand(1)[1]));
la = lap(a)
n = size(la)[1]
b = randn(n)
b = b - mean(b);
In [25]:
@time x = cg(la,b,maxits=4000)
norm(la*x-b)
Out[25]:
Now, let's try a diagonal preconditioner.
In [26]:
d = diag(la)
pre(x) = x ./ d
@time x = pcg(la,b,pre,tol=1e-1,maxits=10^5)
norm(la*x-b)
Out[26]:
In [27]:
@time x = cg(la,b,tol=1e-1,maxits=10^5)
norm(la*x-b)
Out[27]:
It is very different for a random regular graph of the same size
In [28]:
n = 1000000
la = lap(randRegular(n,3))
b = randn(n)
b = b - mean(b);
In [29]:
@time x = cg(la,b,maxits=100)
norm(la*x-b)
Out[29]:
The following is a test of our stretch computation code. I begin by creating a grid graph with random weights, using our stretch computation, and checking that it agrees with the trace computation.
In [30]:
a = grid2(3)
a = uniformWeight(a)
a = a + a';
In [31]:
mst = kruskal(a)
Out[31]:
The following computes a matrix with entries corresponding to the nonzeros of a. For each nonzero, it puts in the stretch. So, to find the total stretch, we should sum them all and then divide by 2.
In [32]:
st = compStretches(mst,a)
Out[32]:
In [33]:
sum(st)/2
Out[33]:
We now check that we got the right answer by using the algebraic formula.
In [34]:
trace( pinv( full(lap(mst))) * lap(a) )
Out[34]:
Now, let's do a speed test on a randomly weighted grid of side 2000.
In [35]:
a = grid2(2000)
@time a = uniformWeight(a)
@time a = a + a';
@time mst = kruskal(a);
@time st = compStretches(mst,a);
Right now there are a multitude of solvers available in the repository. All of them have both calls for SDD systems and Laplacian systems.
In [36]:
SDDSolvers
Out[36]:
In [37]:
LapSolvers
Out[37]:
In [47]:
a = mapweight(grid2(100),x->1/(rand(1)[1]));
n = a.n
la = lap(a);
dval = zeros(n); dval[1] = dval[n] = 1e-3;
sdd = la + spdiagm(dval);
b = randn(n); b = b - mean(b);
The SDD solvers take in the sdd matrix alongside the tolerance, maxits and maxtime parameters.
In [49]:
for solver in SDDSolvers
println("Solver ", solver)
@time f = solver(sdd, maxits=1000, maxtime=10, tol=1e-4, verbose=false)
@time x = f(b);
println("Relative norm: ", norm(sdd * x - b) / norm(b), "\n")
end
The SDD solvers take in the sdd matrix alongside the tolerance, maxits and maxtime parameters.
In [50]:
for solver in LapSolvers
println("Solver ", solver)
@time f = solver(a, maxits=1000, maxtime=10, tol=1e-4, verbose=false)
@time x = f(b);
println("Relative norm: ", norm(sdd * x - b) / norm(b), "\n")
end
In [ ]: