Linear least-squares unmixing with analytical uncertainties

README.md

@@ -103,6 +103,56 @@ Some documented example scripts are given in the directory `examples/`, and are
 python examples/unmix_mwe.py
 ```
 
+## Choosing a solver
+
+Two solvers are available. They differ in what they treat as "misfit", which changes both the
+answer and what you can say about its uncertainty.
+
+| | `SampleNetworkUnmixer` | `LinearSampleNetworkUnmixer` |
+|---|---|---|
+| Misfit | relative (log-ratio surrogate) | absolute (least-squares) |
+| Best for | data spanning orders of magnitude | data with low log-variance, e.g. isotopic |
+| Forward model | convex program | exactly invertible matrix `d = Mc` |
+| Regularization | shrinks towards the mean *observation* | shrinks towards the mean *model* |
+| Uncertainty | Monte Carlo (`solve_montecarlo`) | closed form, `C_c = R C_d R^T` |
+| Resolution | not available | resolution matrix and effective DOF |
+
+The linear solver is documented in Appendix A of the paper. It is the faster and more
+informative of the two when its assumptions hold, because the mixing matrix is square and
+invertible, so the estimator, its covariance, and its resolution are all available in closed
+form:
+
+```python
+problem = funmixer.LinearSampleNetworkUnmixer(sample_network, use_regularization=True)
+solution = problem.solve(
+    element_data,
+    regularization_strength=1.0,
+    data_covariance=10.0,          # a 10% relative error on the observations
+)
+solution.upstream_preds        # recovered source concentrations
+solution.upstream_std          # their 1-sigma uncertainties, no Monte Carlo needed
+solution.downstream_covariance # covariance of the modelled observations
+solution.effective_dof         # how many degrees of freedom the data actually constrain
+```
+
+See `examples/unmix_linear_mwe.py`.
+
+**Check the diagnostics before trusting a linear solution.** Inverting `M` amounts to
+differencing each site against its upstream neighbours, and the noise amplification at each
+site is `Q_i/q_i`, the ratio of total upstream flux to the flux the sub-basin itself generates.
+Where sub-basin areas are very uneven this is large, and the unregularized inversion will
+produce negative or physically impossible concentrations. `solution.amplification`,
+`solution.condition_number`, `solution.clamped_nodes` and `solution.unconstrained_preds` all
+report on this. On the `Mg` example data, the unregularized inversion clamps 10 of 63 sites at
+zero and would otherwise return concentrations up to 145% by mass; at `lambda = 1` none clamp.
+Regularization is not optional for this kind of data.
+
+Note that `regularization_strength` is **not** comparable between the two solvers: the linear
+one weights a squared penalty on deviations from the mean model, the non-linear one an
+unsquared norm of deviations from the mean observation. Within the linear solver, lambda is
+dimensionless (the data are mean-normalised internally), so it does transfer between elements
+and datasets.
+
 ## Cite 
 
 If you use this please cite the paper, which is published at *Water Resources Research*.

docs/.gitignore

@@ -0,0 +1,10 @@
+# LaTeX build artefacts
+*.aux
+*.bcf
+*.blg
+*.fdb_latexmk
+*.fls
+*.log
+*.out
+*.run.xml
+*.bbl

docs/linear_unmixing_note.bib

@@ -0,0 +1,139 @@
+@article{barnes_using_2024,
+  title   = {Using Convex Optimization to Efficiently Apportion Tracer and Pollutant Sources from Point Concentration Observations},
+  author  = {Barnes, Richard and Lipp, Alex G.},
+  journal = {Water Resources Research},
+  year    = {2024},
+  doi     = {10.1029/2023WR036159},
+}
+
+@book{menke_geophysical_2012,
+  title     = {Geophysical Data Analysis: Discrete Inverse Theory},
+  author    = {Menke, William},
+  edition   = {3},
+  publisher = {Academic Press},
+  address   = {Boston},
+  year      = {2012},
+  isbn      = {9780123971609},
+}
+
+@book{aster_parameter_2018,
+  title     = {Parameter Estimation and Inverse Problems},
+  author    = {Aster, Richard C. and Borchers, Brian and Thurber, Clifford H.},
+  edition   = {3},
+  publisher = {Elsevier},
+  year      = {2018},
+  isbn      = {9780128046517},
+}
+
+@book{tarantola_inverse_2005,
+  title     = {Inverse Problem Theory and Methods for Model Parameter Estimation},
+  author    = {Tarantola, Albert},
+  publisher = {Society for Industrial and Applied Mathematics},
+  address   = {Philadelphia},
+  year      = {2005},
+  isbn      = {9780898715729},
+}
+
+@book{hansen_rank-deficient_1998,
+  title     = {Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion},
+  author    = {Hansen, Per Christian},
+  publisher = {Society for Industrial and Applied Mathematics},
+  address   = {Philadelphia},
+  year      = {1998},
+  isbn      = {9780898714036},
+}
+
+@article{hoerl_ridge_1970,
+  title   = {Ridge Regression: Biased Estimation for Nonorthogonal Problems},
+  author  = {Hoerl, Arthur E. and Kennard, Robert W.},
+  journal = {Technometrics},
+  volume  = {12},
+  number  = {1},
+  pages   = {55--67},
+  year    = {1970},
+  doi     = {10.1080/00401706.1970.10488634},
+}
+
+@article{tikhonov_solution_1963,
+  title   = {Solution of Incorrectly Formulated Problems and the Regularization Method},
+  author  = {Tikhonov, Andrey N.},
+  journal = {Soviet Mathematics Doklady},
+  volume  = {4},
+  pages   = {1035--1038},
+  year    = {1963},
+}
+
+@article{elden_algorithms_1977,
+  title   = {Algorithms for the Regularization of Ill-Conditioned Least Squares Problems},
+  author  = {Eld{\'e}n, Lars},
+  journal = {BIT Numerical Mathematics},
+  volume  = {17},
+  number  = {2},
+  pages   = {134--145},
+  year    = {1977},
+  doi     = {10.1007/BF01932285},
+}
+
+@book{golub_matrix_2013,
+  title     = {Matrix Computations},
+  author    = {Golub, Gene H. and Van Loan, Charles F.},
+  edition   = {4},
+  publisher = {Johns Hopkins University Press},
+  address   = {Baltimore},
+  year      = {2013},
+  isbn      = {9781421407944},
+}
+
+@book{boyd_convex_2004,
+  title     = {Convex Optimization},
+  author    = {Boyd, Stephen and Vandenberghe, Lieven},
+  publisher = {Cambridge University Press},
+  year      = {2004},
+  isbn      = {9780521833783},
+}
+
+@article{diamond_cvxpy_2016,
+  title   = {{CVXPY}: A Python-Embedded Modeling Language for Convex Optimization},
+  author  = {Diamond, Steven and Boyd, Stephen},
+  journal = {Journal of Machine Learning Research},
+  volume  = {17},
+  number  = {83},
+  pages   = {1--5},
+  year    = {2016},
+}
+
+@inproceedings{blondes_practical_2016,
+  title     = {A Practical Guide to the Use of Major Elements, Trace Elements, and Isotopes in Compositional Data Analysis: Applications for Deep Formation Brine Geochemistry},
+  author    = {Blondes, M. S. and Engle, M. A. and Geboy, N. J.},
+  booktitle = {Compositional Data Analysis},
+  series    = {Springer Proceedings in Mathematics \& Statistics},
+  publisher = {Springer International Publishing},
+  pages     = {13--29},
+  year      = {2016},
+  isbn      = {978-3-319-44811-4},
+}
+
+@book{aitchison_statistical_1986,
+  title     = {The Statistical Analysis of Compositional Data},
+  author    = {Aitchison, John},
+  publisher = {Chapman and Hall},
+  year      = {1986},
+}
+
+@book{cormen_introduction_2009,
+  title     = {Introduction to Algorithms},
+  author    = {Cormen, Thomas H. and Leiserson, Charles E. and Rivest, Ronald L. and Stein, Clifford},
+  edition   = {3},
+  publisher = {MIT Press},
+  year      = {2009},
+  isbn      = {9780262033848},
+}
+
+@book{horn_matrix_2012,
+  title     = {Matrix Analysis},
+  author    = {Horn, Roger A. and Johnson, Charles R.},
+  edition   = {2},
+  publisher = {Cambridge University Press},
+  year      = {2012},
+  isbn      = {9780521548236},
+}