LU factorization as a recurrence system

LU factorization — Gaussian elimination A = L·U with L unit-lower-triangular and U upper-triangular — is the workhorse direct solver. This page derives it in two steps: first the classic right-looking elimination as a SARE (issue #42, step 1), then a uniformization of the pivot broadcast into a nearest-neighbour propagation (step 2), whose free schedule makes the propagation constraints visible.

Caution

The executable here factors without pivoting. Partial pivoting is a data-dependent row permutation and cannot be a static SURE — see Partial pivoting below.

The Schur elimination as a recurrence

Right-looking (Doolittle) LU eliminates one column per step. Let a(i,j,k) be entry (i,j) after eliminating with pivots 0..k. Step k applies the rank-1 Schur update to the trailing submatrix:

$$a(i,j,k) \;=\; a(i,j,k-1) \;-\; \underbrace{\frac{a(i,k,k-1)}{a(k,k,k-1)}}_{l(i,k)\ \text{multiplier}}\; \underbrace{a(k,j,k-1)}_{u(k,j)\ \text{pivot row}} .$$

Each cell (i,j) is updated for k = 0 .. min(i,j)-1, over the strict trailing domain k ≤ i-1, k ≤ j-1.

A SARE — but with a linear schedule

The three operand taps are affine, not uniform: the pivot a(k,k,k-1), the multiplier column a(i,k,k-1), and the pivot row a(k,j,k-1) are read via index maps that project onto row/column k — a broadcast to the whole trailing submatrix. So LU is a SARE, like the MGS QR.

But here the affine dependence is benign. Every consumer (i,j,k) sits strictly below-and-right of the pivot (k,k,k-1) (the domain guarantees i,j > k), so the dependence vector (i-k, j-k, 1) is always forward. Unlike QR’s reduction-broadcast, a linear schedule exists: τ = [1,1,2] is legal (and so is the free schedule). This is why LU tiles far more readily than QR.

Extracting L and U

Each cell finalizes at k = min(i,j)-1, so the factors leave on super-diagonal faces:

• U(i,j) for j ≥ i leaves at k = i-1 (the face i-k = 1);
• L(i,j) for i > j is the multiplier a(i,j,j-1)/a(j,j,j-1), on j-k = 1.

The first row of U (i=0) is A’s first row and the first column of L (j=0) is A[:,0]/A[0,0] — both read from the k=-1 seed halo (they are never eliminated). The regression test reconstructs the full L, U and verifies L·U = A.

Partial pivoting is not a SURE

Partial pivoting selects, at each step, the row with the largest |a(i,k)| (a per- column iamax) and swaps it to the pivot. That choice is made at runtime from the values, and the resulting row permutation P depends on the data. A SURE’s index space is fixed at compile time, so it cannot express a data-dependent permutation — A = PLU with real pivoting is outside the model. (Numerically, the spec above requires the leading minors to be non-singular, as its bundled matrix is.)

Step 2 — uniformizing the pivot into a propagation

The pivot broadcast is exactly the affine dependence uniformization turns into a pipeline. For a single elimination step (k = 0) the pivot row and multiplier column are the first row/column of the input A, so they can seed nearest-neighbour propagations legally (input confluences read tensors):

$$\begin{aligned} pr(i,j,k) &= pr(i-1,j,k) &&\text{pivot row } A(0,j)\ \text{propagates down } +i\\ mc(i,j,k) &= mc(i,j-1,k) &&\text{multiplier col } A(i,0)\ \text{propagates across } +j\\ r(i,j,k) &= r(i,j,k-1) - g\,\bigl(mc(i,j-1,k)/pv(i-1,j,k)\bigr)\,pr(i-1,j,k) \end{aligned}$$

Every dependence is now a constant offset — a pure SURE (τ = [1,1,1]), applying the rank-1 Schur term exactly once via ger’s inject → add → drain. This is docs/SURE/lu_propagate.sure.

Why only one step. For k > 0 the pivot row a(k,:,k-1) is a computed value, and seeding a propagation from a recurrence variable is an internal confluence the DSL does not have — dfactl rejects it with “input expressions may not read recurrence variables”. This is the same gap that keeps the MGS QR a SARE. So the full multi-step LU stays a SARE; only the single step uniformizes.

What the free schedule shows. Because the pivot must hop i rows down and the multiplier j columns across, cell (i,j) cannot fire until t ≈ max(i,j) — the broadcast’s “everyone at once” becomes a diagonal propagation wavefront. That staircase is the cost the propagation trades for uniformity (it now tiles with halo exchange).

Reduced dependency graph

The reduced dependency graph (RDG) draws the recurrence as one node per variable and one arc per dependence. LU has a single variable a (the matrix eliminated in place), so its RDG is one node with four self-loops — and this is where the SURE/SARE line becomes visible:

arc	dependence map	reads
a → a	translation [0,0,1]ᵀ	a(i,j,k-1), the previous Schur value
a → a	affine (i,j,k) ↦ (i,k,k-1)	the multiplier column a(i,k,k-1)
a → a	affine (i,j,k) ↦ (k,j,k-1)	the pivot row a(k,j,k-1)
a → a	affine (i,j,k) ↦ (k,k,k-1)	the pivot a(k,k,k-1)

The first arc is a translation vector; the other three are matrix maps — each projects every cell of the trailing submatrix onto row/column k, a broadcast. Those three affine arcs are exactly what make LU a SARE, and they are the arcs a neighbour-pivoting reformulation would uniformize into translations. In the graph the affine arcs render dashed, carrying their matrix; the uniform Schur arc renders solid, carrying [0,0,1]ᵀ. (This is the SARE worked example on the RDG page.)

---

The executable specs, their schedules, and the animations of both the SARE and the propagation wavefront are on the reference page: LU decomposition under SURE Algorithms → Linear Algebra.