Doc: added change-of-variable subsection in model-fitting section.

doc/src/Model.rst

@@ -13,7 +13,7 @@ SExtractor can fit models to the images of detected objects since version 2.8. T
 .. math::
   :label: loss_func
 
-  \lambda(\boldsymbol{q}) = \sum_i \left(g\left(\frac{p_i - \tilde{m}_i(\boldsymbol{q})}{\sigma_i}\right)\right)^2 + \sum_j \left(\frac{q_j - \mu_j}{\rho_j}\right)^2
+  \lambda(\boldsymbol{q}) = \sum_i \left(g\left(\frac{p_i - \tilde{m}_i(\boldsymbol{q})}{\sigma_i}\right)\right)^2 + \sum_j \left(\frac{f_j(Q_j) - \mu_j}{s_j}\right)^2
 
 with respect to components of the model parameter vector :math:`\boldsymbol{q}`. :math:`\boldsymbol{q}` comprises parameters describing the shape of the model and the model pixel coordinates :math:`\boldsymbol{x}`.
 
@@ -62,21 +62,72 @@ The vector :math:`\tilde{\boldsymbol{m}}(\boldsymbol{q})` is obtained by convolv
 where :math:`h` is a 2-dimensional interpolant (interpolating function), :math:`\boldsymbol{x}_i` is the coordinate vector of image pixel :math:`i`, :math:`\boldsymbol{x}_j` the coordinate vector of model sample :math:`j`, and :math:`\eta` is the image-to-model sampling step ratio (sampling factor) which is by default defined by the PSF model sampling.
 We adopt a Lánczos-4 function :cite:`duchon1979` as interpolant.
 
+Change of variables
+~~~~~~~~~~~~~~~~~~~
+
+Many model parameters are valid only over a restricted domain.
+Fluxes, for instance, cannot be negative. 
+In order to avoid invalid values and also to facilitate convergence, a change of variables is applied individually to each model parameter:
+
+.. math::
+  :label: change_of_variables
+
+  q_j = f_j(a_j, b_j, Q_j).
+
+The "model" variable :math:`q_j` is bounded by the lower limit :math:`a_j` and the upper limit :math:`b_j` by construction.
+The "engine" variable :math:`Q_j` can take any value, and is actually the parameter that is being adjusted in the fit, although it does not have any physical meaning.
+In |SExtractor| three different types of changes of variables :math:`f_j()` are applied, depending on the parameter (:numref:`change_of_variable_table`).
+
+.. _change_of_variable_table:
+
+.. list-table:: Types of changes of variables applied to model parameters
+  :header-rows: 1
+
+  * - Type
+    - Model :math:`\stackrel{f^{-1}}{\to}` Engine
+    - Engine :math:`\stackrel{f}{\to}` Model
+    - Examples
+  * - Unbounded (linear)
+    - :math:`Q_j = q_j`
+    - :math:`q_j = Q_j`
+    - | :param:`SPHEROID_POSANGLE`
+      | :param:`DISK_POSANGLE`
+  * - Bounded linear
+    - :math:`Q_j = \ln \frac{q_j - a_j}{b_j - q_j}`
+    - :math:`q_j = \frac{b_j - a_j}{1 + \exp -Q_j} + a_j`
+    - | :param:`XMODEL_IMAGE`
+      | :param:`SPHEROID_SERSICN`
+  * - Bounded logarithmic
+    - :math:`Q_j = \ln \frac{\ln q_j - \ln a_j}{\ln b_j - \ln q_j}`
+    - :math:`q_j = a_j \frac{\ln b_j - \ln a_j}{1 + \exp -Q_j}`
+    - | :param:`FLUX_SPHEROID`
+      | :param:`DISK_ASPECT`
+
+In practice, this approach works well, and was found to be much more reliable than a box constrained algorithm :cite:`Kanzow2004375`.
+
 Regularization
 ~~~~~~~~~~~~~~
 
+Although minimizing the (modified) weighted sum of least squares gives a solution that fits best the data, it does not necessarily correspond to the most probable solution given what we know about celestial objects.
+The discrepancy is particularly significant in very faint (|SNR| :math:`\le 20`) and barely resolved galaxies, for which there is, for example, a tendency to overestimate the elongation, known as the "noise bias" in the weak-lensing community :cite:`2004MNRAS_353_529H,2012MNRAS_424_2757M,2012MNRAS_425_1951R,2012MNRAS_427_2711K`.
+To mitigate this issue, |SExtractor| implements a simple `Tikhonov regularization <https://en.wikipedia.org/wiki/Tikhonov_regularization>`_ scheme on some engine parameters, in the form of an additional penalty term in :eq:`loss_func`.
+This term acts as a Gaussian prior on the selected *engine* parameters. However for the associated *model* parameters, the change of variable can make the (improper) prior far from Gaussian.
+
 .. _model_minimization_def:
 
 Minimization
 ~~~~~~~~~~~~
 
 Minimization of the loss function :math:`\lambda(\boldsymbol{q})` is carried out using the `Levenberg-Marquardt algorithm <http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm>`_, and more specifically the |LevMar|_ implementation :cite:`lourakis04LM`.
+The library approximates the Jacobian matrix of the model from finite differences using Broyden's :cite:`Broyden1965ACo` rank one updates.
 The fit is done inside a disk which diameter is scaled to include the isophotal footprint of the object, plus the FWHM of the PSF, plus a 20 % margin.
 The number of iterations is returned in the :param:`NITER_MODEL` measurement parameter.
 It is generally a few tens.
 The final value of the modified chi square term in :eq:`loss_func`, divided by the number of degrees of freedom, is returned in :param:`CHI2_MODEL`.
 The :param:`FLAGS_MODEL` parameter flags various issues which may happen during the fitting process (see the flags section for details on how flags are managed in |SExtractor|):
 
+.. _flags_model_table:
+
 .. csv-table:: :param:`FLAGS_MODEL` flag description
   :header: "Value", "Meaning"
   :widths: 3 60
@@ -87,7 +138,7 @@ The :param:`FLAGS_MODEL` parameter flags various issues which may happen during
   8, "at least one of the fitted parameters hits the lower bound"
   16, "at least one of the fitted parameters hits the upper bound"
 
-:math:`1\,\sigma` error estimates are provided for most measurement parameters; they are obtained by marginalizing the full covariance matrix of the fit.
+:math:`1\,\sigma` error estimates are provided for most measurement parameters; they are obtained from the full covariance matrix of the fit, which is itself computed by inverting the approximate `Hessian matrix <https://en.wikipedia.org/wiki/Hessian_matrix>`_ of :math:`\lambda(\boldsymbol{q})` at the solution.
 
 .. _models_def:
 

doc/src/conf.py

@@ -81,7 +81,7 @@ language = None
 # There are two options for replacing |today|: either, you set today to some
 # non-false value, then it is used:
 #
-today = 'Wed Dec 13 2017'
+today = 'Fri Dec 29 2017'
 #
 # Else, today_fmt is used as the format for a strftime call.
 #

doc/src/keys.rst

@@ -43,6 +43,9 @@
 .. |SExtractor| replace:: :program:`SExtractor`
 .. _SExtractor: http://astromatic.net/software/sextractor
 
+.. |SNR| replace:: :abbr:`SNR (Signal-to-Noise Ratio)`
+.. _SNR: https://en.wikipedia.org/wiki/Signal-to-noise_ratio
+
 .. |SWarp| replace:: :program:`SWarp`
 .. _SWarp: http://astromatic.net/software/swarp
 

doc/src/references.bib

@@ -12,6 +12,15 @@
   adsnote = {Provided by the SAO/NASA Astrophysics Data System}
 }
 
+@ARTICLE{Broyden1965ACo,
+       AUTHOR = "C. Broyden",
+       TITLE = "A Class of Methods for Solving Nonlinear Simultaneous Equations",
+       JOURNAL = "Mathematics of Computation",
+       VOLUME = "19",
+       YEAR = "1965",
+       PAGES = "577-593"
+}
+
 @ARTICLE{2002AA_395_1077C,
    author = {{Calabretta}, M.~R. and {Greisen}, E.~W.},
     title = "{Representations of celestial coordinates in FITS}",
@@ -66,6 +75,21 @@
   adsnote = {Provided by the SAO/NASA Astrophysics Data System}
 }
 
+@ARTICLE{2004MNRAS_353_529H,
+   author = {{Hirata}, C.~M. and {Mandelbaum}, R. and {Seljak}, U. and {Guzik}, J. and {Padmanabhan}, N. and {Blake}, C. and {Brinkmann}, J. and {Bud{\'a}vari}, T. and {Connolly}, A. and {Csabai}, I. and {Scranton}, R. and {Szalay}, A.~S.},
+    title = "{Galaxy-galaxy weak lensing in the Sloan Digital Sky Survey: intrinsic alignments and shear calibration errors}",
+  journal = {MNRAS},
+   eprint = {astro-ph/0403255},
+ keywords = {gravitational lensing, galaxies: haloes},
+     year = 2004,
+    month = sep,
+   volume = 353,
+    pages = {529-549},
+      doi = {10.1111/j.1365-2966.2004.08090.x},
+   adsurl = {http://adsabs.harvard.edu/abs/2004MNRAS.353..529H},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
 @ARTICLE{1987AA_183_177I,
    author = {{Infante}, L.},
     title = "{A faint object processing software - Description and testing}",
@@ -79,6 +103,36 @@
   adsnote = {Provided by the SAO/NASA Astrophysics Data System}
 }
 
+@ARTICLE{2012MNRAS_427_2711K,
+   author = {{Kacprzak}, T. and {Zuntz}, J. and {Rowe}, B. and {Bridle}, S. and  {Refregier}, A. and {Amara}, A. and {Voigt}, L. and {Hirsch}, M.},
+    title = "{Measurement and calibration of noise bias in weak lensing galaxy shape estimation}",
+  journal = {MNRAS},
+archivePrefix = "arXiv",
+   eprint = {1203.5049},
+ primaryClass = "astro-ph.CO",
+ keywords = {gravitational lensing: weak, methods: data analysis, methods: statistical, techniques: image processing, cosmology: observations},
+     year = 2012,
+    month = dec,
+   volume = 427,
+    pages = {2711-2722},
+      doi = {10.1111/j.1365-2966.2012.21622.x},
+   adsurl = {http://adsabs.harvard.edu/abs/2012MNRAS.427.2711K},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@ARTICLE{Kanzow2004375,
+author = {{Kanzow}, C. and {Yamashita}, N. and {Fukushima}, M.},
+title = "{Levenberg-Marquardt methods with strong local convergence properties for solving nonlinear equations with convex constraints}",
+journal = "Journal of Computational and Applied Mathematics ",
+volume = 172,
+number = 2,
+pages = {375-397},
+year = 2004,
+issn = "0377-0427",
+doi = "http://dx.doi.org/10.1016/j.cam.2004.02.013",
+url = "http://www.sciencedirect.com/science/article/pii/S0377042704001256",
+}
+
 @ARTICLE{1980ApJS_43_305K,
    author = {{Kron}, R.~G.},
     title = "{Photometry of a complete sample of faint galaxies}",
@@ -105,6 +159,23 @@
   adsnote = {Provided by the SAO/NASA Astrophysics Data System}
 }
 
+@ARTICLE{2012MNRAS_424_2757M,
+   author = {{Melchior}, P. and {Viola}, M.},
+    title = "{Means of confusion: how pixel noise affects shear estimates for weak gravitational lensing}",
+  journal = {MNRAS},
+archivePrefix = "arXiv",
+   eprint = {1204.5147},
+ primaryClass = "astro-ph.IM",
+ keywords = {gravitational lensing: weak, techniques: image processing },
+     year = 2012,
+    month = aug,
+   volume = 424,
+    pages = {2757-2769},
+      doi = {10.1111/j.1365-2966.2012.21381.x},
+   adsurl = {http://adsabs.harvard.edu/abs/2012MNRAS.424.2757M},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
 @ARTICLE{1856MNRAS_17_12P,
    author = {{Pogson}, N.},
     title = "{Magnitudes of Thirty-six of the Minor Planets for the first day of each month of the year 1857}",
@@ -117,6 +188,23 @@
   adsnote = {Provided by the SAO/NASA Astrophysics Data System}
 }
 
+@ARTICLE{2012MNRAS_425_1951R,
+   author = {{Refregier}, A. and {Kacprzak}, T. and {Amara}, A. and {Bridle}, S. and  {Rowe}, B.},
+    title = "{Noise bias in weak lensing shape measurements}",
+  journal = {MNRAS},
+archivePrefix = "arXiv",
+   eprint = {1203.5050},
+ primaryClass = "astro-ph.CO",
+ keywords = {gravitational lensing: weak, methods: statistical, techniques: image processing, cosmology: observations, dark energy, dark matter },
+     year = 2012,
+    month = sep,
+   volume = 425,
+    pages = {1951-1957},
+      doi = {10.1111/j.1365-2966.2012.21483.x},
+   adsurl = {http://adsabs.harvard.edu/abs/2012MNRAS.425.1951R},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
 @INPROCEEDINGS{2002ASPC_281_228B,
    author = {{Bertin}, E. and {Mellier}, Y. and {Radovich}, M. and {Missonnier}, G. and 
 	{Didelon}, P. and {Morin}, B.},

doc/theme/css/custom.css

@@ -92,3 +92,7 @@ caption:hover a.headerlink {
 	visibility: visible;
 }
 
+/* Remove space before muliline cells in tables */
+td div.line-block {
+	margin-left: 0px !important;
+}