Stable least-squares fitting

Menu Start Research Publications Statistics playground Teaching CV Contact

General background Least-squares problems are omnipresent, occurring in almost every project, from basic data reduction to producing plots for publications. They are mathematically simple to solve on paper, which often leads to naive software implementations that quickly suffer from numerical instabilities (e.g. when there are a lot of data points or many fit parameters). Avoiding numerical instabilities The key to numerically stable least-squares solutions is to avoid matrix inversions or (at least) to exploit clever matrix decompositions. A standard trick to use here is the Cholesky decomposition. For instances, scipy provides an easy-to-use Cholesky decomposition. Alternatively, one may also implement one's own Cholesky decomposition using Householder transformations, which is particularly stable for solving large least-squares systems. For instances, this was done for Hipparcos astrometry and also for Gaia low-resolution BP/RP spectra. Procedure In a first step, one would apply a Cholesky decomposition to the least-squares problem. Second, given the Cholesky decomposition, one can obtain an initial estimate of the solution by inverting the decomposed matrix, but this can still be unstable. Finally, one follows up on the initial estimate with an iterative gradient-descent algorithm which itself does not require any matrix inversion and is thus numerically stable. The gradient descent will thus take us to the global minimum within a few iterations.