simplex: hypersparse FTRAN via activation graph (#20)

ftran mirrors btran's sparse path: for a mostly-zero rhs (nnz*10 <= m)
route through ftranSparse, which fires only the pivots reachable through
the activation graph. Skipped pivots would compute exactly 0, so results
match the dense path bitwise.

buildBtranGraph now also builds ftranOwner (row -> resolving ftran seq),
covering the forward singletons, the shared kernel seq, and the reversed
backward singletons.

Adds TestFactorFtranSparse: 200-trial property test over unit vectors and
1-3 nonzero spikes, checking B*x == v on the sparse path.

Co-authored-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Author: andig

simplex/factor.go

@@ -47,6 +47,10 @@ type factor struct {
 	rdrHead []int32 // row -> [rdrHead[r], rdrHead[r+1]) into rdrList
 	rdrList []int32 // pivot seqs whose column reads that row
 
+	// ftran activation: a written row activates the one pivot that reads it
+	// as its pivot row; ftranOwner maps row -> that ftran seq (kernel shares one).
+	ftranOwner []int32
+
 	ws *factorWS // shared per-LP solve scratch (solver is single-threaded)
 }
 
@@ -336,11 +340,36 @@ func (f *factor) buildBtranGraph() {
 	}
 	head[0] = 0
 	f.rdrHead, f.rdrList = head, list
+
+	// ftran solve order: fwd forward (seq i), kernel (seq nf), bwd reversed
+	// (bwd[i] fires at seq nf+1+(nb-1-i)); own each row by its resolving seq.
+	owner := make([]int32, m)
+	for i, tp := range f.fwd {
+		owner[tp.row] = int32(i)
+	}
+	for _, r := range f.kRows {
+		owner[r] = int32(nf)
+	}
+	for i, tp := range f.bwd {
+		owner[tp.row] = int32(nf + 1 + (nb - 1 - i))
+	}
+	f.ftranOwner = owner
 }
 
 // ftranFactor solves B*x = v in place: v becomes x indexed by basis
 // position (x[pos] = multiplier of basic column pos).
 func (f *factor) ftran(v []float64) {
+	nnz := 0
+	for _, val := range v[:f.m] {
+		if val != 0 {
+			nnz++
+		}
+	}
+	// same activation-graph payoff threshold as btran: mostly-zero rhs
+	if nnz*10 <= f.m {
+		f.ftranSparse(v)
+		return
+	}
 	x := f.ws.xScratch // by basis position
 	clear(x)
 	// forward: row singletons
@@ -390,6 +419,83 @@ func (f *factor) ftran(v []float64) {
 	copy(v, x)
 }
 
+// ftranSparse is ftran firing only activation-reachable pivots; skipped
+// pivots would compute exactly 0, so results match the dense path bitwise.
+func (f *factor) ftranSparse(v []float64) {
+	nf, nb := len(f.fwd), len(f.bwd)
+	kernelSeq := int32(nf)
+	x := f.ws.xScratch // by basis position
+	clear(x)
+	f.ws.gen++
+	gen := f.ws.gen
+	act := f.ws.actGen
+	for r := range f.m {
+		if v[r] != 0 {
+			act[f.ftranOwner[r]] = gen
+		}
+	}
+	// a fired pivot writes its column's rows; activate each written row's owner
+	mark := func(pos int32) {
+		for _, r := range f.tblRow[f.cols[pos]] {
+			act[f.ftranOwner[r]] = gen
+		}
+	}
+	// forward: row singletons
+	for i, tp := range f.fwd {
+		if act[i] != gen {
+			continue
+		}
+		xc := v[tp.row] / tp.a
+		x[tp.pos] = xc
+		if xc != 0 {
+			cr, cv := f.tblRow[f.cols[tp.pos]], f.tblVal[f.cols[tp.pos]]
+			for k, r := range cr {
+				v[r] -= cv[k] * xc
+			}
+			v[tp.row] = 0 // exactly resolved
+			mark(tp.pos)
+		}
+	}
+	// kernel: dense solve when reachable (matches btranSparse's solveT)
+	if k := len(f.kRows); k > 0 && act[kernelSeq] == gen {
+		kv := f.ws.kScratch[:k]
+		for ki, r := range f.kRows {
+			kv[ki] = v[r]
+		}
+		f.klu.solve(kv)
+		for ki := range k {
+			s := kv[ki]
+			pos := f.kPos[ki]
+			x[pos] = s
+			if s != 0 {
+				cr, cv := f.tblRow[f.cols[pos]], f.tblVal[f.cols[pos]]
+				for kk, r := range cr {
+					v[r] -= cv[kk] * s
+				}
+				mark(int32(pos))
+			}
+		}
+	}
+	// backward: col singletons in reverse
+	for i := nb - 1; i >= 0; i-- {
+		seq := int32(nf + 1 + (nb - 1 - i))
+		if act[seq] != gen {
+			continue
+		}
+		tp := f.bwd[i]
+		xc := v[tp.row] / tp.a
+		x[tp.pos] = xc
+		if xc != 0 {
+			cr, cv := f.tblRow[f.cols[tp.pos]], f.tblVal[f.cols[tp.pos]]
+			for k, r := range cr {
+				v[r] -= cv[k] * xc
+			}
+			mark(tp.pos)
+		}
+	}
+	copy(v, x)
+}
+
 // btran solves B^T*y = w in place: w is indexed by basis position on entry,
 // y by row on exit. Ops are the adjoints of ftran's, in reverse order.
 func (f *factor) btran(w []float64) {

simplex/factor_test.go

@@ -54,6 +54,48 @@ func multiply(colRow [][]int32, colVal [][]float64, m int, x []float64) []float6
 	return v
 }
 
+// TestFactorFtranSparse drives the mostly-zero rhs that routes through
+// ftranSparse (unit vectors and 1-3 nonzeros), checking B*x == v.
+func TestFactorFtranSparse(t *testing.T) {
+	rng := rand.New(rand.NewSource(13))
+	for trial := range 200 {
+		m := 5 + rng.Intn(80)
+		colRow, colVal := randomBasis(rng, m)
+		f := factorize(m, identCols(m), colRow, colVal, nil)
+		if f == nil {
+			continue
+		}
+		v := make([]float64, m)
+		want := make([]float64, m)
+		// unit vector on even trials (the alpha/tau case), else 1-3 spikes
+		if trial%2 == 0 {
+			v[rng.Intn(m)] = 1
+		} else {
+			for n := 1 + rng.Intn(3); n > 0; n-- {
+				v[rng.Intn(m)] = rng.NormFloat64()
+			}
+		}
+		copy(want, v)
+		nnz := 0
+		for _, x := range v {
+			if x != 0 {
+				nnz++
+			}
+		}
+		if nnz*10 > m {
+			continue // would route dense; not the path under test
+		}
+		f.ftran(v)
+		back := multiply(colRow, colVal, m, v)
+		for i := range back {
+			if math.Abs(back[i]-want[i]) > 1e-7 {
+				t.Fatalf("trial %d m=%d sparse ftran: B*x[%d]=%g want %g (kernel %d fwd %d bwd %d)",
+					trial, m, i, back[i], want[i], len(f.kRows), len(f.fwd), len(f.bwd))
+			}
+		}
+	}
+}
+
 func TestFactorFtranBtran(t *testing.T) {
 	rng := rand.New(rand.NewSource(7))
 	for trial := range 50 {