Electrical Networks and Graph Theory, in depth: Kirchhoff's laws as a Laplacian system, effective resistance, Thomson, Rayleigh and Foster, harmonic functions and solving at scale

From Kirchhoff's laws to the Laplacian

Model a circuit as an undirected graph: vertices are junctions and each edge (u, v) is a resistor with resistance r_uv, or equivalently conductance c_uv = 1/r_uv. Give every vertex a potential φ(v). Three rules fix everything:

• Ohm's law: the current from u to v along their edge is (φ(u) - φ(v)) · c_uv.
• Kirchhoff's current law: at every vertex, the current flowing out equals the current injected there from outside, b(v). Interior junctions have b(v) = 0.
• Kirchhoff's voltage law: drops around a cycle sum to zero. Writing currents as potential differences satisfies it automatically.

Substituting Ohm into the current law at vertex v gives the sum over neighbours u of c_uv(φ(v) - φ(u)) = b(v). In matrix form that is L φ = b, where the Laplacian L has the weighted degree on the diagonal and -c_uv off it. L is symmetric and positive semidefinite, and on a connected graph its null space is exactly the constant vectors, because adding the same voltage everywhere changes no current. Two consequences follow. The injections must sum to zero, since current that goes in must come out, and the potentials are fixed only up to a constant, so you ground one vertex by setting its potential to 0 and deleting its row and column. The reduced matrix is positive definite and the solve is unique.

Worked example: a five-resistor bridge

Take four junctions S, A, B and T with resistors S-A of 1 Ω, S-B of 2 Ω, A-T of 2 Ω, B-T of 1 Ω, and a 1 Ω edge between A and B, the middle edge of a Wheatstone bridge. Inject 1 A at S, take it out at T, and ground T. The current law at A says (φA - φS)/1 + φA/2 + (φA - φB)/1 = 0, at B it says (φB - φS)/2 + φB/1 + (φB - φA)/1 = 0, and at S the outflow is 1. Solving gives φS = 1.4, φA = 0.8 and φB = 0.6 volts.

Check the current law at A: 0.6 A arrives from S, 0.4 A leaves toward T and 0.2 A crosses to B. B receives 0.4 A from S plus 0.2 A from A and sends 0.6 A to T. The effective resistance R_eff(S, T) is the potential difference needed to push one ampere, here 1.4 Ω: the whole network behaves like a single 1.4 Ω resistor between S and T.

Solving it in code

The code builds L from an edge list, grounds the sink, solves, and reads currents off the potentials. The dense solve is right for a few thousand vertices; the section on scale replaces it.

import numpy as np

def laplacian(n, edges):
    """edges: (u, v, resistance in ohms). Conductance c = 1/r."""
    L = np.zeros((n, n))
    for u, v, r in edges:
        c = 1.0 / r
        L[u, u] += c; L[v, v] += c
        L[u, v] -= c; L[v, u] -= c
    return L

def solve_potentials(n, edges, s, t, current=1.0):
    """Inject `current` amps at s, extract at t. Returns node potentials and edge currents."""
    L = laplacian(n, edges)
    b = np.zeros(n); b[s] = current; b[t] = -current
    keep = [i for i in range(n) if i != t]           # ground t: phi[t] = 0
    phi = np.zeros(n)
    phi[keep] = np.linalg.solve(L[np.ix_(keep, keep)], b[keep])
    flows = [(u, v, (phi[u] - phi[v]) / r) for u, v, r in edges]
    return phi, flows

def effective_resistance(n, edges, s, t):
    phi, _ = solve_potentials(n, edges, s, t)
    return phi[s] - phi[t]

def harmonic_extension(n, edges, fixed):
    """fixed: {node: value}. Every other node gets the weighted average of its neighbours."""
    L = laplacian(n, edges)
    free = [i for i in range(n) if i not in fixed]
    bnd = list(fixed)
    x = np.zeros(n)
    x[bnd] = [fixed[i] for i in bnd]
    rhs = -L[np.ix_(free, bnd)] @ x[bnd]
    x[free] = np.linalg.solve(L[np.ix_(free, free)], rhs)
    return x

S, A, B, T = 0, 1, 2, 3
edges = [(S, A, 1), (S, B, 2), (A, T, 2), (B, T, 1), (A, B, 1)]
phi, flows = solve_potentials(4, edges, S, T)    # phi = [1.4, 0.8, 0.6, 0.0]

It reproduces the hand solution: potentials 1.4, 0.8, 0.6 and 0 volts, currents 0.6, 0.4, 0.4, 0.6 and 0.2 amps, and the net current at A and B is zero to within 3e-16. The same number comes from the Moore-Penrose pseudoinverse: R_eff(s, t) = (e_s - e_t)^T L⁺ (e_s - e_t), which also gives 1.4. Use the pseudoinverse when you need resistances between many pairs in a small graph, since one O(n³) factorisation answers every query; use grounded solves for a few pairs in a large graph.

Effective resistance and its properties

Four facts make effective resistance a working tool rather than a curiosity. All four were checked on the example.

• Thomson's principle. Among all ways to route one unit of flow from s to t, the electrical flow minimises the energy, the sum of r_e · f_e², and the minimum equals R_eff(s, t). In the example the energy is 0.36 + 0.32 + 0.32 + 0.36 + 0.04 = 1.4. Any other routing of the same ampere dissipates more.
• Rayleigh monotonicity. Lowering any edge's resistance can only lower R_eff, and raising it can only raise it. Delete the middle edge and the network becomes two parallel 3 Ω paths, R_eff = 1.5 Ω; halve the middle edge to 0.5 Ω and it drops to 1.375 Ω. Monotonicity lets you bound a hard network by an easier one: delete edges for an upper bound, short vertices together for a lower bound.
• It is a metric. R_eff satisfies the triangle inequality, so it defines the resistance distance on a graph. Unlike shortest-path distance, it shrinks when many parallel routes exist, so it rewards redundancy.
• Foster's theorem. On a connected graph with n vertices, the sum over edges of c_e · R_eff(e) equals n - 1. The example gives exactly 3.0. Each term is the probability that the edge appears in a random spanning tree, so the theorem says a spanning tree has n - 1 edges; see random spanning trees.

Series and parallel rules are special cases of the linear system: resistances add in series and conductances add in parallel. The Wheatstone example is the smallest network that cannot be reduced by those two rules alone, which is why it needs the solve.

Potentials are random-walk probabilities

Ground t at 0 and hold s at 1 volt. At every other vertex the current law says the potential is the conductance-weighted average of its neighbours' potentials: the function is harmonic off the boundary. Now consider a random walk that moves along each edge with probability proportional to its conductance. The probability that the walk from v reaches s before t satisfies exactly the same averaging equation with the same boundary values, and the solution is unique, so the two are equal.

In the example, harmonic_extension with S fixed at 1 and T at 0 gives A = 0.571 and B = 0.429: a walk from A reaches S before T with probability 4/7. On a path of four unit resistors with the ends fixed at 1 and 0 the values are 0.75, 0.5 and 0.25, the gambler's-ruin answer. The same identity gives the commute time between s and t as 2m · R_eff(s, t), derived in random walks.

Machine learning uses this directly. In graph-based semi-supervised learning (Zhu, Ghahramani and Lafferty, 2003), labelled nodes are the fixed boundary, edge weights encode similarity, and the harmonic extension assigns every unlabelled node a score in [0, 1]. It is one linear solve, with no training loop.

Solving at scale

A dense solve costs O(n³) time and n² memory, which stops at a few thousand vertices. Real networks are sparse, and L x can be computed edge by edge, so conjugate gradient applies: it needs only matrix-vector products and converges because L is positive semidefinite. Starting from zero with a right-hand side that sums to zero, every iterate stays orthogonal to the constant vector, so no grounding is needed.

def laplacian_matvec(n, us, vs, cs):
    """y = L x without building L: one pass over the edge arrays."""
    def mv(x):
        d = cs * (x[us] - x[vs])               # current on each edge
        y = np.zeros(n)
        np.add.at(y, us, d)
        np.add.at(y, vs, -d)
        return y
    return mv

def cg_laplacian(mv, b, tol=1e-8, max_iter=10000):
    """Conjugate gradient for L x = b on a connected graph; returns the mean-zero solution."""
    b = b - b.mean()                           # b must be orthogonal to the all-ones vector
    x = np.zeros_like(b); r = b.copy(); p = r.copy(); rs = r @ r
    for it in range(max_iter):
        Ap = mv(p)
        alpha = rs / (p @ Ap)
        x += alpha * p; r -= alpha * Ap
        rs_new = r @ r
        if np.sqrt(rs_new) <= tol * np.linalg.norm(b):
            return x - x.mean(), it + 1
        p = r + (rs_new / rs) * p; rs = rs_new
    raise RuntimeError("CG did not converge: graph disconnected or badly conditioned?")

Two lessons sit in that table. Iterations roughly double each time the side doubles, because the condition number of a grid Laplacian grows with the square of the side and CG needs iterations proportional to its square root. And the corner-to-corner resistance grows only logarithmically, by about 0.88 Ω per doubling, because the number of parallel routes grows as the grid does. On the five-resistor network the same CG code returns 1.4 Ω in two iterations.

For large or badly conditioned graphs, precondition. Jacobi scaling by the degree is free; incomplete Cholesky and algebraic multigrid cut iterations sharply on meshes. In theory, Spielman and Teng showed that Laplacian systems can be solved in nearly linear time, and Christiano and colleagues (2011) built approximate max flow on repeated electrical-flow solves. For exact max flow in practice, augmenting-path methods such as Dinic's algorithm remain the standard.

Where the system appears

• Chip power grids. IR-drop analysis models a chip's power delivery mesh as millions of resistors with current sinks at the cells and solves L φ = b to find where the supply voltage sags. It is this article's system with very large n, solved with preconditioned iterative methods.
• DC power flow. The standard linear approximation of an AC power grid sets line flows proportional to phase-angle differences, giving the same Laplacian system with line susceptances as weights.
• Graph sparsification. Spielman and Srivastava sample each edge with probability proportional to c_e · R_eff(e) and reweight it; Foster's theorem says those weights sum to n - 1, so O(n log n / ε²) samples give a sparse graph whose Laplacian approximates the original within a factor of 1 ± ε.
• Robustness metrics. The sum of effective resistances over all vertex pairs, the Kirchhoff index, is lower in networks with redundant paths, so it is a better resilience score than the average hop count.

Failure modes

• Disconnected graph. Each component adds a zero eigenvalue, the grounded matrix is singular, and CG stalls. Find components first and solve each separately; resistance between components is infinite.
• Injections that do not sum to zero. L φ = b then has no solution. The CG code silently projects b, which hides the bug, so assert the sum in the caller.
• Zero or tiny resistances. Conductance becomes infinite or huge. Merge zero-resistance edges into one vertex before solving, and expect slow convergence when conductances span many orders of magnitude.
• Negative weights. L stops being positive semidefinite and every theorem above fails. Signed graphs need different tools.
• Directed edges. Resistors are symmetric. A directed graph has a non-symmetric Laplacian, and effective resistance does not carry over directly.
• Resistance versus conductance. Mixing the two inverts the weights and quietly produces plausible numbers. Name variables by unit and check one hand-solvable case, such as two resistors in series.

What to do next

1. Run the code on the five-resistor network and check the potentials and the 1.4 Ω result by hand.
2. Verify Thomson, Rayleigh and Foster numerically on a random weighted graph of 50 vertices.
3. Use harmonic_extension to label a small similarity graph from three seed labels and compare with a nearest-neighbour baseline.
4. Scale the CG solver to a 500 by 500 grid, then add Jacobi preconditioning and compare iteration counts.
5. Read the matrix-tree theorem and random spanning trees to see the same Laplacian count and sample spanning trees.