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]:
Base.SparseMatrix.CHOLMOD.Factor{Float64}
type:          LLt
method: simplicial
maxnnz:        102
nnz:           102

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]:
5.143087378554464e-15

Let's poke around to see what F has inside it.

In [4]:
F
Out[4]:
Base.SparseMatrix.CHOLMOD.Factor{Float64}
type:          LLt
method: simplicial
maxnnz:        102
nnz:           102
In [5]:
super(typeof(F))
Out[5]:
Factorization{Float64}
In [6]:
isa(F,Factorization)
Out[6]:
true
In [7]:
fieldnames(F)
Out[7]:
1-element Array{Symbol,1}:
 :p

Testing the speed of that

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)
1.527121 seconds (205 allocations: 256.844 MB, 6.59% gc time)
  
Out[9]:
2.976278365938353e-10

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);
0.377724 seconds (19 allocations: 7.634 MB, 1.35% gc time)
  

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]:
2.069855996116447e-12

What about just using \ ?

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)
1.434061 seconds (59 allocations: 259.422 MB, 0.89% gc time)
  0.053558 seconds (19 allocations: 648.281 KB)
  
Out[12]:
3.047155955006307e-10

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)
1.980881 seconds (139.37 k allocations: 274.019 MB, 4.79% gc time)
Out[14]:
5.809541011026003e-9

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]:
2.621755646216747e-9

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]:
(anonymous function)
In [17]:
f = lapChol(la)
Out[17]:
(anonymous function)
In [18]:
norm(la*f(b) - b)
Out[18]:
2.621755646216747e-9
In [19]:
x = lapWrapSolver(cholfact,la,b)
norm(la*x - b)
Out[19]:
2.621755646216747e-9
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]:
6.530265309919312e-6
In [21]:
bbig = convert(Array{BigFloat,1},b);
xbig = cg(a,bbig,maxits=100)
norm(a*xbig - bbig)
Out[21]:
3.408504794566293816762470763584253827015914440348506316254697226173169945005103e-37
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]:
0.00019928918001516186
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)
0.816791 seconds (29.34 k allocations: 207.333 MB, 3.30% gc time)
  
Out[25]:
97.96603417977055

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)
4.563222 seconds (16.01 k allocations: 1.193 GB, 3.50% gc time)
  
Out[26]:
19.882642385705193
In [27]:
@time x = cg(la,b,tol=1e-1,maxits=10^5)
norm(la*x-b)
1.198206 seconds (484.48 k allocations: 307.148 MB, 3.73% gc time)
  
Out[27]:
19.590092408356856

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)
7.891222 seconds (27.15 k allocations: 2.024 GB, 3.64% gc time)
  
Out[29]:
0.0007745578882358395

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]:
9x9 sparse matrix with 16 Float64 entries:
	[2, 1]  =  0.718277
	[4, 1]  =  1.30955
	[1, 2]  =  0.718277
	[3, 2]  =  0.656996
	[2, 3]  =  0.656996
	[6, 3]  =  0.640175
	[1, 4]  =  1.30955
	[6, 5]  =  1.22709
	[3, 6]  =  0.640175
	[5, 6]  =  1.22709
	[9, 6]  =  0.587062
	[8, 7]  =  0.360802
	[7, 8]  =  0.360802
	[9, 8]  =  1.02741
	[6, 9]  =  0.587062
	[8, 9]  =  1.02741
2.854674 seconds (172 allocations: 328.068 MB, 4.00% gc time)

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]:
9x9 sparse matrix with 24 Float64 entries:
	[2, 1]  =  1.0
	[4, 1]  =  1.0
	[1, 2]  =  1.0
	[3, 2]  =  1.0
	[5, 2]  =  6.45151
	[2, 3]  =  1.0
	[6, 3]  =  1.0
	[1, 4]  =  1.0
	[5, 4]  =  12.0766
	[7, 4]  =  14.5912
	⋮
	[8, 5]  =  6.03562
	[3, 6]  =  1.0
	[5, 6]  =  1.0
	[9, 6]  =  1.0
	[4, 7]  =  14.5912
	[8, 7]  =  1.0
	[5, 8]  =  6.03562
	[7, 8]  =  1.0
	[9, 8]  =  1.0
	[6, 9]  =  1.0
	[8, 9]  =  1.0
In [33]:
sum(st)/2
Out[33]:
47.15492426761423

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]:
47.154924267614184

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);
  2.406615 seconds (31.98 M allocations: 1.609 GB, 15.55% gc time)
  1.107632 seconds (35 allocations: 823.610 MB, 17.57% gc time)
  8.919166 seconds (94 allocations: 1.103 GB, 7.60% gc time)
  
In [ ]:

Here is my a simple solver that uses an augmented spanning tree. It should never be too bad. It wants a positive definite system.

In [36]:
a = mapweight(grid2(500),x->1/(rand(1)[1]));
la = lap(a)
n = size(la)[1]
la[1,1] = la[1,1] + 1
@time F = augTreeSolver(la,tol=1e-1,maxits=1000)
b = randn(n)
@time x = F(b)
norm(la*x - b)
3.746014 seconds (79 allocations: 1.467 GB, 15.42% gc time)
  2.534656 seconds (1.74 M allocations: 508.631 MB, 12.85% gc time)
  
Out[36]:
49.93092418950264
In [37]:
@time F = augTreeSolver(la,tol=1e-1,maxits=1000,treeAlg=randishPrim)
b = randn(n)
@time x = F(b)
norm(la*x - b)
4.106294 seconds (5.89 k allocations: 3.254 GB, 9.19% gc time)
  1.520148 seconds (1.08 M allocations: 523.226 MB, 18.93% gc time)
  
Out[37]:
49.82364163609152

Compare to the running time of CG

In [38]:
@time y = cg(la,b,tol=1e-1,maxits=1000)
norm(la*y-b)
4.837599 seconds (6.04 k allocations: 3.597 GB, 8.72% gc time)
  
Out[38]:
519.1489086730564
In [39]:
n = 40000
la = lap(randRegular(n,3))
la[1,1] = la[1,1] + 1
@time F = augTreeSolver(la)
b = randn(n)
@time x = F(b)
norm(la*x - b)
7.434149 seconds (4.01 k allocations: 1.868 GB, 4.28% gc time)
  0.174792 seconds (121.18 k allocations: 66.463 MB, 4.97% gc time)
  
Out[39]:
20.639305575286002

Compare to the running time of a direct method.

In [40]:
@time Fc = cholfact(la)
@time x = Fc \ b
norm(la*x - b)
0.299919 seconds (2.52 k allocations: 247.258 MB, 6.75% gc time)
  8.628694 seconds (68 allocations: 1015.664 MB, 0.04% gc time)
  
Out[40]:
1.1468469752650618e-11

Now, let's try wrapping for singular systems.

In [41]:
n = 40000
la = lap(randRegular(n,3))
f = lapWrapSolver(augTreeSolver,la,tol=1e-6,maxits=1000)
Out[41]:
(anonymous function)
0.128699 seconds (19 allocations: 1.266 MB)
In [42]:
a = grid2(1000)
t = randishKruskal(a);
st = compStretches(t,a);
sum(st)/nnz(a)
Out[42]:
48.24403503503503
In [43]:
t = randishPrim(a);
st = compStretches(t,a);
sum(st)/nnz(a)
Out[43]:
32.09173773773774