From 10061b32c17fc3a6c8b925845909c64e8de48b4b Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Mon, 19 Oct 2026 14:26:38 -0600 Subject: [PATCH 1/2] Expand the README and documentation, following ndfilters. Both now have: - A runnable quickstart, which fits a Gaussian to a thousand noisy lines, each starting at the peak of its line, and evaluates the best-fit models on a finer grid. The documentation plots three of them, and the README shows that figure from Read the Docs. - A guide to writing a model, with an exponential decay as the example. The documentation fits it to a hundred decays with bounds and plots the fitted scales against the truth. - The differences from `scipy.optimize.least_squares()`: the Jacobian, noise and priors, bounds, the convergence test, and what isn't supported. - A "Custom models" item in the list of features. - A citation section pointing to `CITATION.cff`, with a BibTeX entry. The Zenodo DOI will be added after the first release. The documentation also lists the functions of the package. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PZRXiFhGzzrF87vsFg73B9 --- README.md | 136 +++++++++++++++++++++++++++----- docs/index.rst | 206 ++++++++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 310 insertions(+), 32 deletions(-) diff --git a/README.md b/README.md index 9f06046..81ce2fd 100644 --- a/README.md +++ b/README.md @@ -14,10 +14,11 @@ once, similar to but compiled and parallelized using [Numba](https://numba.readthedocs.io/en/stable/). -Its main function, `ndoptimize.least_squares()`, fits a model to every -spectrum of a spectral cube, or any other batch of small datasets, using the -Levenberg-Marquardt method. Compared to calling -`scipy.optimize.least_squares()` in a loop, it offers: +Its main function, +[`ndoptimize.least_squares()`](https://ndoptimize.readthedocs.io/en/latest/_autosummary/ndoptimize.least_squares.html), +fits a model to every spectrum of a spectral cube, or any other batch of small +datasets, using the Levenberg-Marquardt method. +Compared to calling `scipy.optimize.least_squares()` in a loop, it offers: - **Speed.** The fits run in parallel on every CPU core, roughly half a million fits of a five-parameter model to 39 points per second on a @@ -30,44 +31,139 @@ Levenberg-Marquardt method. Compared to calling - **Missing data.** Data points which are NaN are left out of their fit. - **Broadcasting.** The data, noise, initial guess, bounds and priors are broadcast against each other, so each fit can have its own. +- **Custom models.** A model is a short Numba function which computes its + values and their derivatives. Every model shares the same fitting kernel, + which is compiled once and cached on disk. + +## Differences from `scipy.optimize.least_squares` + +For the same objective, the two find the same minimum, to a thousandth of a +standard deviation in the tests, but they differ in the following ways. + +- **Derivatives.** A model computes its own Jacobian. There is no + finite-difference approximation. +- **Noise and priors.** The variance of the data, the correlation of its + noise, and the Gaussian priors are arguments, and the result includes the + covariance of the parameters. `scipy.optimize.least_squares()` minimizes the + sum of squares of whatever residual it is given, and returns the Jacobian at + the solution, from which the covariance can be computed. +- **Bounds.** Steps are clipped to the bounds, and a parameter which the + gradient pushes against its bound is held fixed for that step. SciPy's + `"lm"` method does not support bounds, and its `"trf"` and `"dogbox"` + methods handle them with trust regions. +- **Convergence.** A single tolerance on the Newton decrement, the decrease in + chi-squared that a Gauss-Newton step would give, replaces `ftol`, `xtol` + and `gtol`. +- **Not supported:** robust loss functions and sparse Jacobians. + +The full documentation is hosted on +[Read the Docs](https://ndoptimize.readthedocs.io/en/latest/). ## Installation +`ndoptimize` is published on PyPI and can be installed using `pip`. + ```bash pip install ndoptimize ``` -## Example +## Quickstart + +Fit a Gaussian to each of a thousand noisy spectral lines at once, and +evaluate the best-fit models on a finer grid. ```python import numpy as np import ndoptimize +rng = np.random.default_rng(0) + +# a thousand noisy Gaussian lines with random centers velocity = np.linspace(-100, 100, 41) -data = ... # spectra, shape (..., 41) -variance = ... # their noise variance, the same shape +center = rng.normal(0, 10, size=(1000, 1)) +sigma = 0.5 +data = 5 * np.exp(-np.square((velocity - center) / 15) / 2) +data = data + rng.normal(0, sigma, size=data.shape) + +# start each fit at the peak of its line +guess = np.stack([ + data.max(axis=-1), # amplitude + velocity[data.argmax(axis=-1)], # center + np.full(len(data), 10), # width +], axis=-1) +# fit every line at once result = ndoptimize.least_squares( - model=ndoptimize.models.gaussian_linear, + model=ndoptimize.models.gaussian, x=velocity[:, np.newaxis], data=data, - variance=variance, - guess=[5, 0, 15, 0, 1], # amplitude, center, width, slope, intercept + variance=sigma**2, + guess=guess, ) -result.parameters # shape (..., 5) -result.covariance # shape (..., 5, 5) +result.parameters # shape (1000, 3) +result.covariance # shape (1000, 3, 3) -# the best-fit model at each point of each spectrum, shape (..., 41) +# the best-fit model of every line on a finer grid, shape (1000, 201) +velocity_fine = np.linspace(-100, 100, 201) fit = ndoptimize.evaluate( - model=ndoptimize.models.gaussian_linear, - x=velocity[:, np.newaxis], + model=ndoptimize.models.gaussian, + x=velocity_fine[:, np.newaxis], parameters=result.parameters, ) ``` -Models are short Numba functions which compute the model and its Jacobian, -compiled with the `ndoptimize.model` decorator, which also names their -parameters. -See the [documentation](https://ndoptimize.readthedocs.io/en/latest) for -details. +![quickstart](https://ndoptimize.readthedocs.io/en/latest/_images/index_0_0.png) + +## Writing a model + +A model is a function of the points `x`, an array of shape +`(num_points, num_inputs)`, and the `parameters`, which writes the model into +`f` and its derivatives with respect to the parameters into `jacobian`. +The +[`ndoptimize.model`](https://ndoptimize.readthedocs.io/en/latest/_autosummary/ndoptimize.model-decorator.html) +decorator compiles it and names its parameters. + +```python +import math +import numpy as np +import ndoptimize + +@ndoptimize.model(parameters=("amplitude", "scale")) +def exponential( + x: np.ndarray, + parameters: np.ndarray, + f: np.ndarray, + jacobian: np.ndarray, +) -> None: + """An exponential decay, amplitude * exp(-x / scale).""" + amplitude, scale = parameters[0], parameters[1] + for i in range(x.shape[0]): + e = math.exp(-x[i, 0] / scale) + f[i] = amplitude * e + jacobian[i, 0] = e + jacobian[i, 1] = amplitude * e * x[i, 0] / scale**2 +``` + +The compiled model is cached on disk, so a model defined in the interactive +interpreter, rather than in a file or a notebook, needs `cache=False`. +[`ndoptimize.models`](https://ndoptimize.readthedocs.io/en/latest/_autosummary/ndoptimize.models.html) +has some common models: a Gaussian, a Gaussian plus a straight line, and two +Gaussians plus a straight line. + +## Citation + +If you use ndoptimize in your research, please cite it. +The citation metadata is kept in [`CITATION.cff`](https://github.com/sun-data/ndoptimize/blob/main/CITATION.cff), +which the "Cite this repository" button on GitHub can export as BibTeX or APA. +Please include the version of ndoptimize that you used, +which is given by `importlib.metadata.version("ndoptimize")`. + +```bibtex +@software{ndoptimize, + author = {Smart, Roy T.}, + title = {ndoptimize}, + version = {X.Y.Z}, + url = {https://github.com/sun-data/ndoptimize}, +} +``` diff --git a/docs/index.rst b/docs/index.rst index 524fe86..deaad47 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -9,19 +9,57 @@ just-in-time compiler. Its main function, :func:`ndoptimize.least_squares`, fits a model to every spectrum of a spectral cube, or any other batch of small datasets, using the Levenberg-Marquardt method. -Each fit has its own independent variable, data, noise, initial guess, -bounds and Gaussian priors, and the noise of neighboring data points may be -correlated, as it is in resampled data, which makes the fit generalized least -squares. -The fits run in parallel on every CPU core, roughly half a million small -fits per second on a 48-thread machine. +Compared to calling :func:`scipy.optimize.least_squares` in a loop, it +offers: -Models are short Numba functions which compute the model and its Jacobian, -compiled with the :func:`ndoptimize.model` decorator, which also names their -parameters. -Some common ones are provided in :mod:`ndoptimize.models`, and -:func:`ndoptimize.evaluate` evaluates a model for many sets of parameters at -once, such as the best-fit ones. +- **Speed.** The fits run in parallel on every CPU core, roughly half a + million fits of a five-parameter model to 39 points per second on a + 48-thread machine. +- **Correlated noise.** The noise of neighboring data points may be + correlated, as it is in resampled data, which makes the fit generalized + least squares, with parameter uncertainties to match. +- **Priors and bounds.** Each parameter of each fit may have a Gaussian prior + and bounds, to keep faint fits stable. +- **Missing data.** Data points which are NaN are left out of their fit. +- **Broadcasting.** The data, noise, initial guess, bounds and priors are + broadcast against each other, so each fit can have its own. +- **Custom models.** A model is a short Numba function which computes its + values and their derivatives. Every model shares the same fitting kernel, + which is compiled once and cached on disk. + +The package provides: + +- :func:`ndoptimize.least_squares`, which fits a model to many sets of data + at once. +- :func:`ndoptimize.evaluate`, which evaluates a model for many sets of + parameters at once, such as the best-fit ones. +- :func:`ndoptimize.model`, a decorator which compiles a function as a model, + an :class:`ndoptimize.Model`, and names its parameters. +- :mod:`ndoptimize.models`, some common models: a Gaussian, a Gaussian plus a + straight line, and two Gaussians plus a straight line. + + +Differences from :func:`scipy.optimize.least_squares` +===================================================== + +For the same objective, the two find the same minimum, to a thousandth of a +standard deviation in the tests, but they differ in the following ways. + +- **Derivatives.** A model computes its own Jacobian. There is no + finite-difference approximation. +- **Noise and priors.** The variance of the data, the correlation of its + noise, and the Gaussian priors are arguments, and the result includes the + covariance of the parameters. :func:`scipy.optimize.least_squares` + minimizes the sum of squares of whatever residual it is given, and returns + the Jacobian at the solution, from which the covariance can be computed. +- **Bounds.** Steps are clipped to the bounds, and a parameter which the + gradient pushes against its bound is held fixed for that step. SciPy's + ``"lm"`` method does not support bounds, and its ``"trf"`` and + ``"dogbox"`` methods handle them with trust regions. +- **Convergence.** A single tolerance on the Newton decrement, the decrease + in chi-squared that a Gauss-Newton step would give, replaces ``ftol``, + ``xtol`` and ``gtol``. +- **Not supported:** robust loss functions and sparse Jacobians. Installation @@ -31,6 +69,150 @@ Installation pip install ndoptimize +Quickstart +========== + +Fit a Gaussian to each of a thousand noisy spectral lines at once, and +evaluate the best-fit models on a finer grid. + +.. jupyter-execute:: + + import numpy as np + import matplotlib.pyplot as plt + import ndoptimize + + rng = np.random.default_rng(0) + + # a thousand noisy Gaussian lines with random centers + velocity = np.linspace(-100, 100, 41) + center = rng.normal(0, 10, size=(1000, 1)) + sigma = 0.5 + data = 5 * np.exp(-np.square((velocity - center) / 15) / 2) + data = data + rng.normal(0, sigma, size=data.shape) + + # start each fit at the peak of its line + guess = np.stack([ + data.max(axis=-1), # amplitude + velocity[data.argmax(axis=-1)], # center + np.full(len(data), 10), # width + ], axis=-1) + + # fit every line at once + result = ndoptimize.least_squares( + model=ndoptimize.models.gaussian, + x=velocity[:, np.newaxis], + data=data, + variance=sigma**2, + guess=guess, + ) + + # the best-fit model of every line on a finer grid + velocity_fine = np.linspace(-100, 100, 201) + fit = ndoptimize.evaluate( + model=ndoptimize.models.gaussian, + x=velocity_fine[:, np.newaxis], + parameters=result.parameters, + ) + + # plot the first three lines and their fits + fig, ax = plt.subplots(constrained_layout=True) + for i in range(3): + ax.scatter(velocity, data[i], color=f"C{i}", s=10); + ax.plot(velocity_fine, fit[i], color=f"C{i}"); + ax.set_xlabel("velocity (km/s)"); + ax.set_ylabel("intensity"); + +The result holds the best-fit parameters of every line, their covariance, +chi-squared, and the convergence state of each fit, as described in +:class:`ndoptimize.LeastSquaresResult`. + +.. jupyter-execute:: + + result.parameters.shape, result.covariance.shape + + +Writing a model +=============== + +A model is a function of the points ``x``, an array of shape +``(num_points, num_inputs)``, and the ``parameters``, which writes the model +into ``f`` and its derivatives with respect to the parameters into +``jacobian``. +The :func:`ndoptimize.model` decorator compiles it and names its parameters. +Here is an exponential decay, fit to a hundred noisy decays, with bounds which +keep the amplitude and the scale positive. + +.. jupyter-execute:: + + import math + + @ndoptimize.model(parameters=("amplitude", "scale")) + def exponential( + x: np.ndarray, + parameters: np.ndarray, + f: np.ndarray, + jacobian: np.ndarray, + ) -> None: + """An exponential decay, amplitude * exp(-x / scale).""" + amplitude, scale = parameters[0], parameters[1] + for i in range(x.shape[0]): + e = math.exp(-x[i, 0] / scale) + f[i] = amplitude * e + jacobian[i, 0] = e + jacobian[i, 1] = amplitude * e * x[i, 0] / scale**2 + + time = np.linspace(0, 10, 50) + scale = rng.uniform(1, 5, size=(100, 1)) + data = 3 * np.exp(-time / scale) + rng.normal(0, 0.1, size=(100, time.size)) + + result = ndoptimize.least_squares( + model=exponential, + x=time[:, np.newaxis], + data=data, + variance=0.1**2, + guess=[1, 1], + lower=[0, 0.1], + ) + + # compare the fitted scales to the truth + fig, ax = plt.subplots(constrained_layout=True) + ax.errorbar( + x=scale[:, 0], + y=result.parameters[:, 1], + yerr=np.sqrt(result.covariance[:, 1, 1]), + fmt="o", + markersize=3, + ); + ax.plot([1, 5], [1, 5], color="gray"); + ax.set_xlabel("true scale"); + ax.set_ylabel("fitted scale"); + +The compiled model is cached on disk, so a model defined in the interactive +interpreter, rather than in a file or a notebook, needs ``cache=False``. + + +Citation +======== + +If you use :mod:`ndoptimize` in your research, please cite it. +The citation metadata is kept in +`CITATION.cff `_, +which the "Cite this repository" button on the +`GitHub page `_ +can export as BibTeX or APA. +Please include the version of :mod:`ndoptimize` that you used, +which is given by ``importlib.metadata.version("ndoptimize")``. + +.. code-block:: bibtex + + @software{ndoptimize, + author = {Smart, Roy T.}, + title = {ndoptimize}, + version = {X.Y.Z}, + url = {https://github.com/sun-data/ndoptimize}, + } + + API Reference ============= From 8b65ce68c72c198b5d17f3f1e13226971d08de6c Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Mon, 19 Oct 2026 14:51:14 -0600 Subject: [PATCH 2/2] Fix the review findings in the README and documentation. - Say that each fit can have its own independent variable, and that the correlation of the noise is shared by every fit. - Add the differences from `scipy.optimize.least_squares()` in model arguments (extra data goes in more columns of `x`), in `chi2` versus `cost`, and in counting accepted steps versus evaluations. - Call the convergence quantity the decrease in chi-squared that a Gauss-Newton step predicts, which is half the square of the Newton decrement, here and in the `tolerance` docstring. `max_iterations` and `iterations` now say they count accepted steps. - Put the note on `cache=False` before the model example, which fails when pasted into the interactive interpreter. - Rename the noise `sigma` of the quickstart to `noise`, since sigma is the width of the Gaussian model, and note the shape of `x`. - Say that the README figure comes from the documentation, and pin its cell first in `index.rst`, since its name depends on the cell order. - Fix the description of the bounds of the exponential example and the grammar of the definition of a model. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PZRXiFhGzzrF87vsFg73B9 --- README.md | 39 ++++++++++++++++++----------- docs/index.rst | 48 ++++++++++++++++++++++++------------ ndoptimize/_least_squares.py | 26 ++++++++++--------- 3 files changed, 72 insertions(+), 41 deletions(-) diff --git a/README.md b/README.md index 81ce2fd..ede654f 100644 --- a/README.md +++ b/README.md @@ -29,8 +29,9 @@ Compared to calling `scipy.optimize.least_squares()` in a loop, it offers: - **Priors and bounds.** Each parameter of each fit may have a Gaussian prior and bounds, to keep faint fits stable. - **Missing data.** Data points which are NaN are left out of their fit. -- **Broadcasting.** The data, noise, initial guess, bounds and priors are - broadcast against each other, so each fit can have its own. +- **Broadcasting.** The independent variable, data, noise variance, initial + guess, bounds and priors are broadcast against each other, so each fit can + have its own. The correlation of the noise is shared by every fit. - **Custom models.** A model is a short Numba function which computes its values and their derivatives. Every model shares the same fitting kernel, which is compiled once and cached on disk. @@ -51,9 +52,16 @@ standard deviation in the tests, but they differ in the following ways. gradient pushes against its bound is held fixed for that step. SciPy's `"lm"` method does not support bounds, and its `"trf"` and `"dogbox"` methods handle them with trust regions. -- **Convergence.** A single tolerance on the Newton decrement, the decrease in - chi-squared that a Gauss-Newton step would give, replaces `ftol`, `xtol` - and `gtol`. +- **Convergence.** A single tolerance on the decrease in chi-squared that a + Gauss-Newton step predicts replaces `ftol`, `xtol` and `gtol`. +- **Model arguments.** A model takes only the points `x` and the parameters, + so any other data it needs, such as the edges of each bin, goes in more + columns of `x`. SciPy passes extra data through `args` and `kwargs`. +- **Results.** `chi2` is the sum of the squares of the whitened residuals, + without the priors, while SciPy's `cost` is half the sum of the squares of + whatever residual it is given. `max_iterations` and `iterations` count + accepted steps, while SciPy's `max_nfev` and `nfev` count evaluations of + the residual. - **Not supported:** robust loss functions and sparse Jacobians. The full documentation is hosted on @@ -81,9 +89,9 @@ rng = np.random.default_rng(0) # a thousand noisy Gaussian lines with random centers velocity = np.linspace(-100, 100, 41) center = rng.normal(0, 10, size=(1000, 1)) -sigma = 0.5 +noise = 0.5 data = 5 * np.exp(-np.square((velocity - center) / 15) / 2) -data = data + rng.normal(0, sigma, size=data.shape) +data = data + rng.normal(0, noise, size=data.shape) # start each fit at the peak of its line guess = np.stack([ @@ -95,9 +103,9 @@ guess = np.stack([ # fit every line at once result = ndoptimize.least_squares( model=ndoptimize.models.gaussian, - x=velocity[:, np.newaxis], + x=velocity[:, np.newaxis], # shape (num_points, num_inputs) data=data, - variance=sigma**2, + variance=noise**2, guess=guess, ) @@ -113,16 +121,21 @@ fit = ndoptimize.evaluate( ) ``` +The documentation plots three of the lines with their fits: + ![quickstart](https://ndoptimize.readthedocs.io/en/latest/_images/index_0_0.png) ## Writing a model -A model is a function of the points `x`, an array of shape -`(num_points, num_inputs)`, and the `parameters`, which writes the model into -`f` and its derivatives with respect to the parameters into `jacobian`. +A model is a function, `model(x, parameters, f, jacobian)`, which writes its +value at the points `x`, an array of shape `(num_points, num_inputs)`, into +`f`, and its derivatives with respect to the `parameters` into `jacobian`. The [`ndoptimize.model`](https://ndoptimize.readthedocs.io/en/latest/_autosummary/ndoptimize.model-decorator.html) decorator compiles it and names its parameters. +The compiled model is cached on disk, so define it in a file or a notebook. +A model defined in the interactive interpreter needs +`@ndoptimize.model(..., cache=False)`. ```python import math @@ -145,8 +158,6 @@ def exponential( jacobian[i, 1] = amplitude * e * x[i, 0] / scale**2 ``` -The compiled model is cached on disk, so a model defined in the interactive -interpreter, rather than in a file or a notebook, needs `cache=False`. [`ndoptimize.models`](https://ndoptimize.readthedocs.io/en/latest/_autosummary/ndoptimize.models.html) has some common models: a Gaussian, a Gaussian plus a straight line, and two Gaussians plus a straight line. diff --git a/docs/index.rst b/docs/index.rst index deaad47..824b7dc 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -21,8 +21,9 @@ offers: - **Priors and bounds.** Each parameter of each fit may have a Gaussian prior and bounds, to keep faint fits stable. - **Missing data.** Data points which are NaN are left out of their fit. -- **Broadcasting.** The data, noise, initial guess, bounds and priors are - broadcast against each other, so each fit can have its own. +- **Broadcasting.** The independent variable, data, noise variance, initial + guess, bounds and priors are broadcast against each other, so each fit can + have its own. The correlation of the noise is shared by every fit. - **Custom models.** A model is a short Numba function which computes its values and their derivatives. Every model shares the same fitting kernel, which is compiled once and cached on disk. @@ -56,9 +57,18 @@ standard deviation in the tests, but they differ in the following ways. gradient pushes against its bound is held fixed for that step. SciPy's ``"lm"`` method does not support bounds, and its ``"trf"`` and ``"dogbox"`` methods handle them with trust regions. -- **Convergence.** A single tolerance on the Newton decrement, the decrease - in chi-squared that a Gauss-Newton step would give, replaces ``ftol``, - ``xtol`` and ``gtol``. +- **Convergence.** A single tolerance on the decrease in chi-squared that a + Gauss-Newton step predicts replaces ``ftol``, ``xtol`` and ``gtol``. +- **Model arguments.** A model takes only the points ``x`` and the + parameters, so any other data it needs, such as the edges of each bin, goes + in more columns of ``x``. SciPy passes extra data through ``args`` and + ``kwargs``. +- **Results.** :attr:`~ndoptimize.LeastSquaresResult.chi2` is the sum of the + squares of the whitened residuals, without the priors, while SciPy's + ``cost`` is half the sum of the squares of whatever residual it is given. + ``max_iterations`` and :attr:`~ndoptimize.LeastSquaresResult.iterations` + count accepted steps, while SciPy's ``max_nfev`` and ``nfev`` count + evaluations of the residual. - **Not supported:** robust loss functions and sparse Jacobians. @@ -75,6 +85,11 @@ Quickstart Fit a Gaussian to each of a thousand noisy spectral lines at once, and evaluate the best-fit models on a finer grid. +.. + The README embeds the figure of the next cell from Read the Docs as + _images/index_0_0.png, a name which jupyter-sphinx derives from the order + of the cells, so keep it the first jupyter-execute cell of this page. + .. jupyter-execute:: import numpy as np @@ -86,9 +101,9 @@ evaluate the best-fit models on a finer grid. # a thousand noisy Gaussian lines with random centers velocity = np.linspace(-100, 100, 41) center = rng.normal(0, 10, size=(1000, 1)) - sigma = 0.5 + noise = 0.5 data = 5 * np.exp(-np.square((velocity - center) / 15) / 2) - data = data + rng.normal(0, sigma, size=data.shape) + data = data + rng.normal(0, noise, size=data.shape) # start each fit at the peak of its line guess = np.stack([ @@ -100,9 +115,9 @@ evaluate the best-fit models on a finer grid. # fit every line at once result = ndoptimize.least_squares( model=ndoptimize.models.gaussian, - x=velocity[:, np.newaxis], + x=velocity[:, np.newaxis], # shape (num_points, num_inputs) data=data, - variance=sigma**2, + variance=noise**2, guess=guess, ) @@ -134,13 +149,17 @@ chi-squared, and the convergence state of each fit, as described in Writing a model =============== -A model is a function of the points ``x``, an array of shape -``(num_points, num_inputs)``, and the ``parameters``, which writes the model -into ``f`` and its derivatives with respect to the parameters into +A model is a function, ``model(x, parameters, f, jacobian)``, which writes its +value at the points ``x``, an array of shape ``(num_points, num_inputs)``, +into ``f``, and its derivatives with respect to the ``parameters`` into ``jacobian``. The :func:`ndoptimize.model` decorator compiles it and names its parameters. +The compiled model is cached on disk, so define it in a file or a notebook. +A model defined in the interactive interpreter needs +``@ndoptimize.model(..., cache=False)``. + Here is an exponential decay, fit to a hundred noisy decays, with bounds which -keep the amplitude and the scale positive. +keep the amplitude from going negative and the scale at least 0.1. .. jupyter-execute:: @@ -187,9 +206,6 @@ keep the amplitude and the scale positive. ax.set_xlabel("true scale"); ax.set_ylabel("fitted scale"); -The compiled model is cached on disk, so a model defined in the interactive -interpreter, rather than in a file or a notebook, needs ``cache=False``. - Citation ======== diff --git a/ndoptimize/_least_squares.py b/ndoptimize/_least_squares.py index 28e3605..97b75d0 100644 --- a/ndoptimize/_least_squares.py +++ b/ndoptimize/_least_squares.py @@ -342,7 +342,8 @@ def _levenberg_marquardt( max_iterations The maximum number of steps per fit, not negative. tolerance - The Newton decrement below which a fit has converged. + The predicted decrease in chi-squared below which a fit has + converged. num_threads The number of chunks to split the fits into, the result of :func:`numba.get_num_threads`, which disables caching if it is called @@ -440,9 +441,9 @@ def _levenberg_marquardt( at_upper = p[k] >= upper[i, k] and g[k] > 0 free[k] = a[k, k] != 0 and not (at_lower or at_upper) - # the Newton decrement of the free parameters, the decrease - # in chi-squared a Gauss-Newton step would give if the model - # were linear + # the decrease in chi-squared that a Gauss-Newton step of the + # free parameters predicts, half the square of the Newton + # decrement for k in range(num_parameter): rhs[k] = g[k] if free[k] else 0 for q in range(num_parameter): @@ -576,7 +577,7 @@ class LeastSquaresResult: """The number of data points used by each fit.""" iterations: np.ndarray - """The number of steps taken by each fit.""" + """The number of accepted steps taken by each fit.""" converged: np.ndarray """Whether each fit converged.""" @@ -756,13 +757,16 @@ def least_squares( It must describe a positive-definite noise covariance. If :obj:`None` (the default), the noise is uncorrelated. max_iterations - The maximum number of steps for each fit. + The maximum number of accepted steps for each fit. + A step which does not decrease :math:`\chi^2` is retried with more + damping, which evaluates the model again without counting as a step. tolerance - The convergence tolerance on the Newton decrement, - :math:`g^T A^{-1} g`, where :math:`g = J^T C^{-1} r` is the gradient - and :math:`A = J^T C^{-1} J` is the normal matrix, which is the - decrease in :math:`\chi^2` a Gauss-Newton step would give if the - model were linear. + The convergence tolerance on :math:`g^T A^{-1} g`, the decrease in + :math:`\chi^2` that a Gauss-Newton step predicts, where + :math:`g = J^T C^{-1} r` is minus half the gradient of + :math:`\chi^2` and :math:`A = J^T C^{-1} J` is the normal matrix, + both including the priors. + This is half the square of the Newton decrement of :math:`\chi^2`. The default, :math:`10^{-8}`, puts each parameter within about :math:`10^{-4}` of its uncertainty of the minimum.