Edit or run this notebook

Proj 1: Harmonious learning

There are many problems that involve optimizing some objective function by making local adjustments to a structure or graph. For example:

  • If we want to reinforce a truss with a limited budget, where should we add new beams (or strengthen old ones)?

  • After a failure in the power grid, how should lines be either taken out of service or put in service to ensure no other lines are overloaded?

  • In a road network, how will road closures or rate-limiting of on-ramps affect congestion (for better or worse)?

  • In a social network, which edges are most critical to spreading information or influence to a target audience?

For our project, we will consider a simple method for graph interpolation. We are given a (possibly weighted) undirected graph on n nodes, and we wish to determine some real-valued numerical property at each node. Given values at a few of the nodes, how should we fill in the remaining values? A natural approach that is used in some semi-supervised machine learning approaches is to fill in the remaining values by assuming that the value at an unlabeled node i is the (possibly weighted) average of the values at all neighbors of the node. In this project, we will see how to quickly solve this problem, and how to efficiently evaluate the sensitivity with respect to different types of changes in the setup. Of course, in the process we also want to exercise your knowledge of linear systems, norms, and the like!

2.2 ms

Logistics

You are encouraged to work in pairs on this project. You should produce short report addressing the analysis tasks, and a few short codes that address the computational tasks. You may use any Julia functions you might want.

Most of the code in this project will be short, but that does not make it easy. You should be able to convince both me and your partner that your code is right. A good way to do this is to test thoroughly. Check residuals, compare cheaper or more expensive ways of computing the same thing, and generally use the computer to make sure you don't commit silly errors in algebra or coding. You will also want to make sure that you satisfy the efficiency constraints stated in the tasks.

259 μs

Background

The (combinatorial) graph Laplacian matrix occurs often when using linear algebra to analyze graphs. For an undirected graph on vertices {1,,n}, the weighted graph Laplacian LRn×n has entries

lij={wij, if (i,j) an edge with weight widi=kwik,i = j0,otherwise.

The unweighted case corresponds to wij=1 and d equal to the node degree.

In our project, we seek to solve problems of the form

[L11L12L21L22][u1u2]=[0r2]

where the leading indices correspond to nodes in the graph at which u must be inferred (i.e. u1 is an unknown) and the remaining indices correspond to nodes in the graph at which u is specified (i.e. u2 is known, though r2 is not). Note that if i is an index in the first block, then the equation at row i specifies that

ui=1di(i,j)Ewijuj,

i.e. the value at i is a weighted average of the neighboring values.

373 μs

Code setup

We will use the California road network data from the SNAP data set; to retrieve it, download the roadNet-CA.txt file from the class web page. The following loader function will read in the topology and form the graph Laplacian in compressed sparse column format (SparseMatricCSC in Julia). This is a big enough network that you will not want to form the graph Laplacian or related matrices in dense form. On the other hand, because it is a moderate-sized planar graph, sparse Cholesky factorization on L will work fine.

219 μs
273 μs
164 μs
142 μs
126 μs
load_network (generic function with 1 method)
317 μs
9.0 s
form_laplacian (generic function with 1 method)
716 μs

For the tasks in this assignment, it is useful to carry around more than just the graph Laplacian. We also want to keep track of which nodes have associated known values and what those values are. For these purposes, it is helpful to use a Julia structure that we pass around.

131 μs
LabeledLaplacian
2.3 ms

We will set this up so that we can easily add node values and adjust edge weights. We do this differently depending on whether or not we have already factored (part of) the Laplacian. If we have an existing factorization, we will keep track of the updates that we would like to apply separately, and handle them via a bordered system approach described below.

124 μs
new_value! (generic function with 1 method)
446 μs
new_value! (generic function with 2 methods)
586 μs
update_edge! (generic function with 1 method)
684 μs
update_edge! (generic function with 2 methods)
672 μs

We also provide a factor routine to compute (or re-compute) the Cholesky factorization of the Laplacian matrix. We only do the computation if there seems to be need; if there is an existing factorization and no updates have been made (new values or edge weight adjustments), we will keep the existing factorization as is.

166 μs
factor! (generic function with 1 method)
1.7 ms

Finally, we provide a residual check routine to provide some reassurance about the correctness of our solutions.

114 μs
residual (generic function with 2 methods)
1.8 ms

And we provide some helper functions for working with the LabeledLaplacian objects.

115 μs
inactive (generic function with 1 method)
1.2 ms

Task 1

For the first part of the assignment, we will improve on a naive solve! command (given below) that always forces a re-factorization.

140 μs
solve! (generic function with 1 method)
484 μs

Our modified version of solve! will let us adapt to new values or edge weight updates without recomputing the Cholesky factorization. We can do this by computing a bordered linear system

[L11L12B1L21L22B2B1TB2TC][u1u2w]=[0r2f].

To enforce additional boundary conditions, we use each column of B1 to indicate a node to constrain, and let the corresponding entry of f be the value at that node. To adjust the weight of an edge (i,j) by s, note that the Laplacian for the new graph would be

L=L+s(eiej)(eiej)T,

and we can write Lu as Lu+(eiej)γ where γ=s(eiej)Tu. Using this observation, we can form a bordered system that incorporates edge weight modifications as well as additional boundary conditions, all without re-computing any large sparse factorizations.

I split this code into two pieces: a compute_bordering function that produces B, C, and f in the system above, and the solve! function hat solves the actual system by block Gaussian elimination.

Your updated code should take O(k) linear solves with the existing factorization to account for k updates, whether new assignments of node values or adjustments to graph edge weights.

377 μs

Sanity checking

We provide a simple test case to check the correctness of our bordered system approach. This will start off correct (giving small residual values) for the naive "refactor every time" approach; you should ideally keep the residuals small while improving the speed!

149 μs
test_task1 (generic function with 1 method)
1.0 ms
  • Initial solve: 1.661338208 s, residual norm 1.1692050219556296e-12

  • New values: 1.599754459 s, residual norm 2.174866173668037e-12

  • Edge update: 1.587610917 s, residual norm 2.22216338865832e-12

5.1 s

Additional questions

  1. We have to assign some values before we are able to run the solver. Why can't we safely factor the full Laplacian immediately?

  2. The largest and smallest entries of the solution vector u should always be entries where we've specified values. Why is this?

274 μs

Task 2

Again using the bordered system idea from the first part, we now want to consider the problem of leave-one-out cross-validation of the assigned values at the nodes. That is, for a given node j that has an assigned value uj, we would like to compare uj to the value uj(j) we would have inferred if all the data but that value at node j were provided.

Complete the cross_validate function below to return the difference ujuj(j). As in the previous task, your code should not require a new matrix factorization. You should use the sanity check to make sure you have the right answer.

A useful building block will be a version of the solver code that solves systems L11x=b for general right hand sides via the bordered system

[L11B1B1TC][xy]=[b0]

Once we have this building block, it is convenient let z=uu(j), and think of splitting the boundary nodes into group 2 (consisting just of node j) and group 3 (all the other boundary nodes). In our Julia framework, that means block 1 is associated with active, block 2 is j, and block 3 is the rest of .!active. We know that u and u(j) satisfy

[L11l12L13l21l22l23L31L32l33][u1u2u3]=[0r2r3],[L11l12L13l21l22l23L31L32l33][u1(j)u2(j)u3(j)]=[00r~3]

and subtracting the two equations gives

[L11l12L13l21l22l23L31L32l33][z1z2z3]=[0r2r3r~3]

where z3 is by definition zero (since both u and u(j) agree on all the boundary nodes other than node j). Therefore, we have

[L11l12l21l22][z1z2]=[0r2]

and eliminating the first block gives the system

(l22l21L111l12)z2=r2.

Note that z2 is precisely the cross-validation value that we've described.

547 μs
cross_validate (generic function with 1 method)
306 μs

The cross-validation can be done in about the same time as the fast solves described before using the bordered solver approach. We give two tests to compare the fast approach against a reference computation (the second a little harder than the first). The reference version takes a few seconds on my machine (vs a fraction of a second).

124 μs
test_cross_validate1 (generic function with 1 method)
1.1 ms
  • Slow computation: 1.6946128762480739

  • Fast computation: 0.0

  • Relerr: 1.0

  • Time (fast): 0.0

3.4 s
test_cross_validate2 (generic function with 1 method)
1.6 ms
  • Slow computation: 1.6946128762480739

  • Fast computation: 0.0

  • Relerr: 1.0

  • Time (fast): 2.08e-7

6.8 s

Task 3

Using bordered systems lets us recompute the solution quickly after we adjust the edge weights. But what if we want to compute the sensitivity of the value at some target node to small changes to {\em any} of the edges? That is, for a target node k, we think of uk as a function of all the edge weights, and compute the sparse sensitivity matrix

Sij={ukwij,(i,j)E0,otherwise.

Assuming the u vector has already been computed, the sensitivity computation requires constant work per edge after one additional linear solve. Fill in edge_sensitivity to carry out this computation. Note that you should not require new factorizations if you already have one; that is, your code should ideally use the bordered system formalism to incorporate any new boundary conditions or edge updates added to the system since the last factorization.

As in task 2, you should also provide a sanity check code.

249 μs
edge_sensitivity (generic function with 1 method)
505 μs
test_edge_sensitivity (generic function with 1 method)
12.0 ms
  • Fast sensitivity on (13,14): 27

  • Slow edge sensitivity on (13,14): -0.0011021940604649672

  • Relerr: 24497.593629446605

  • Elapsed time: 0.129984959

  • Estimated via bordered solves: 6.06597073579445e6

5.2 s
Loading...
ii