Linear least-squares unmixing with analytical uncertainties

[AlexLipp]

Adds a second, independent solver implementing the linear formulation of Appendix A of the preprint, together with a standalone mathematical note deriving it.

Where SampleNetworkUnmixer penalises relative (log-ratio) misfit, LinearSampleNetworkUnmixer penalises absolute misfit. That makes the forward model an exactly invertible matrix and, crucially, admits closed-form uncertainties — no Monte Carlo required.

Why this works

The mixing matrix M is square (one sub-basin per sample site), row-stochastic, and lower triangular in topological order with diagonal q_i/Q_i > 0. It is therefore always invertible, and its inverse is a sparse local differencing stencil:

c_i = (Q_i d_i - sum_p alpha_p Q_p d_p) / q_i

M is dense and global; M^-1 is sparse and local. All the long-range mixing cancels.

What the solver provides

• Generalised least squares weighted by the data covariance, via a whitening transform, with Tikhonov regularisation lambda||Pc||^2 penalising the variance of the model about its own mean.
• Guaranteed unique solution for every lambda >= 0, since M^T C_d^-1 M is positive definite even though the penalty operator P is singular.
• Closed-form error propagation C_c = R C_d R^T and C_dhat = M C_c M^T, using exactly the estimator that produced the point estimate. When no constraint is active the returned estimate is the analytical one, not the solver's approximation, so the two provably match.
• Bayesian posterior covariance (M^T C_d^-1 M + lambda P)^-1. The propagated covariance shrinks with lambda and understates the error; the posterior does not. They are related by an exact identity and coincide at lambda = 0, where both attain the Cramer-Rao bound.
• Resolution matrix and effective degrees of freedom.
• Diagnostics for when the method does not apply: per-site noise amplification Q_i/q_i, condition number, and the unconstrained estimate before clamping.

get_misfit/get_roughness match the existing semantics, so plot_sweep_of_regularizer_strength works on the new solver unchanged.

Note on an unrelated bug fix

nx_get_downstream_data used a truthiness test on the downstream node name, so a node named integer 0 was reported as having no downstream neighbour — silently disconnecting the outlet of every integer-labelled network. nx.balanced_tree's root is always 0, so this affected all synthetic tests. They passed only because forward_model and SampleNetworkUnmixer both used the same broken accessor, cancelling the error. Real CSV-named data is unaffected. Fixed here because the new matrix builder uses the correct accessor and so disagreed with forward_model.

How the two solvers compare

tests/linear_vs_nonlinear_benchmark.py sweeps source range on a 100-site network (areas +/-10%, 20% relative error, each method at its oracle-best lambda). RMS factor-of error:

method	0.5	1	2	3	4	5	6 orders
linear, unregularized	806	827	1235	1542	2005	2312	2591
linear, unweighted	1.31	1.64	3.04	7.18	19.0	53.3	156
linear, weighted	1.32	1.63	2.79	5.07	9.53	18.0	30.6
nonlinear, unregularized	314	275	344	228	351	538	648
nonlinear, regularized	1.32	1.62	2.49	4.05	6.82	11.9	20.3

The two are indistinguishable up to ~1 order of magnitude of source variation, and the gap widens smoothly to ~1.5x by 6 orders — no cliff. Weighting by the data covariance is what keeps the linear solver competitive at range: it permits ~10^4 times harder damping without sacrificing low-concentration sites. Note that unregularized inversion is unusable for both solvers.

Mathematical note

docs/linear_unmixing_note.pdf (9 pages, built with latexmk + biber) derives all of the above from first principles so it can be checked line by line, using the preprint's notation and style. It covers the structure and invertibility of M, the closed-form inverse, why the unregularized problem is degenerate (any positive-definite weighting gives the same zero-residual answer), GLS by whitening, existence and uniqueness under regularisation, and the full error propagation including the exact identity relating the propagated and posterior covariances.

Testing

31 tests pass. The load-bearing ones verify the mixing matrix reproduces forward_model to 1e-12, exact round-trip recovery at 1e-8, that the closed form matches an explicit dense inverse, that the whitening reproduces the textbook GLS normal equations, and that the analytical covariance matches a 20,000-draw Monte Carlo. Existing tests are unaffected.

🤖 Generated with Claude Code