Graph-based semi-supervised learning (SSL) propagates a few labels over a similarity graph by minimizing a Dirichlet-type energy. The standard quadratic (
p=2) energy reduces to a single graph-Laplacian solve, but it degenerates exactly where SSL is most useful when labels are scarce: gathering more unlabeled data drives the
p=2 estimate to a near-constant function whenever
d≥2 (Nadler-Srebro-Zhou). Well-posedness requires the nonlinear
p-Laplacian energy with
p>d. Existing solvers reduce this to a sequence of weighted Laplacian solves, but their reference implementations use a direct sparse factorization or ichol-preconditioned CG instead. Plugging a near-linear Laplacian solver is not straightforward: at large
p the conductance weights degenerate near flat-gradient edges, making the system nearly singular and causing stagnation without a damped outer iteration. We close this gap. Recasting
p-Laplacian SSL as a source-form nonlinear Laplacian flow
Bρp(B⊤x)=b and solving by damped chord-Newton continuation in
p, every linearized system stays well-conditioned and can be delegated to a near-linear Laplacian engine. On size-scaled graph families the wall-clock is empirically
m0.96-
m1.02 per family (approximate Cholesky default), and a pooled fit across 228 SuiteSparse graphs gives
m1.19 vs.\
m1.45 for direct factorization; the solver handles a
6.8×107-edge social network in minutes. Memory is the binding constraint: Cholesky fill reaches
10-
280× the graph nonzeros vs.\ our
O(m) hierarchy. Against the released FCL solver we are
1.5-
14× faster at matched accuracy. On MNIST
10-NN,
p=3 scores
64% at one label per class vs.\
36% for
p=2. Code: https://github.com/orenlivne/np.