Prerequisite
- Some knowledge of Matrix operations
- The background information provided in Note 224
- Fitting module
Etymology
From Least SQuares and Fitting
Synopsis
The LSQFit class contains the basic functions to do all the fitting described in the Note about fitting. It handles real, and complex equations;linear and non-linear (Levenberg-Marquardt) solutions;
regular (with optional constraints) or Singular Value Decomposition (SVD).
In essence they are a set of routines to generate normal equations (makeNorm()) in triangular form from a set of condition equations;
to do a Cholesky-type decomposition of the normal equations (either regular or SVD) and test its rank (invert());
to do a quasi inversion of the decomposed equations (solve()) to obtain the solution and/or the errors.
All calculations are done in place. Methods to obtain additional information about the fitting process are available.
This class can be used as a stand-alone class outside of the aips++
environment. In that case the aips.h include file
can be replaced if necessary by appropriate typedefs for Double, Float and
uInt.
The interface to the methods have standard data or standard STL iterator
arguments only. They can be used with any container having an STL
random-access iterator interface. Especially they can be used with
carrays (necessary templates provided),
aips++ Vectors (necessary templates
provided in LSQaips),
standard random access STL containers (like std::vector).
The normal operation of the class consists of the following steps:
- Create an LSQFit object.
The information that can be provided in the constructor of the object,
either directly, or indirectly using the set() commands, is
(see Note 224):
- The number of unknowns that have to be solved for (mandatory)
- The number of constraint equations you want to use explicitly (defaults to 0, but can be changed on-the-fly)
- A collinearity test factor (defaults to 1e-8)
The collinearity factor is the square of the sine of the angle between
a column in the normal equations, and the hyper-plane through
all the other columns. In special cases (e.g. fitting a polynomial
in a very narrow bandwidth window, it could be advisable to set this
factor to zero if you want a solution (whatever the truth of it maybe).
- A Levenberg-Marquardt adjustment factor (if appropriate, defaults to 1e-3)
- Create the normal equations used in solving the set of condition
equations of the user, by using the makeNorm() methods.
Separate makenorm() methods are provided for sparse condition
equations (e.g. if data for 3 antennas are provided, rather than for all 64)
- If there are user provided constraints, either limiting constraints like
the sum of the angles in a triangle is 180 degrees, or constraints to add
missing information if e.g. only differences between parameters have been
measured, they can be added to the normal
equations with the setConstraint() or
the addConstraint() methods. Lagrange multipliers will be used to
solve the extended normal equations.
- The normal equations are triangu;arised (using the collinearity factor
as a check for solvability) with the invert() method. If the
normal equations are non-solvable an error is returned, or a switch to
an SVD solution is made if indicated in the invert call.
- The solutions and adjustment errors are obtained with the
solve() method.
A non-linear loop in a Levenberg-Marquardt adjustment can be obtained
(together with convergence information), with the solveLoop()
method (see below) replacing the combination of
invert and solve.
- Non-linear loops are done by looping through the data using
makeNorm() calls, and upgrade the solution with the
solveLoop() method.
The normal equations are upgraded by changing LM factors. Upgrade depends
on the 'balanced' factor. The LM factor is either added in some way to all
diagonal elements (if balanced) or all diagonal elements are multiplied by
(1+factor) After each loop convergence can be tested
by the isReady() call; which will return False or
a non-zero code indicating the ready reason. Reasons for stopping can be:
- SOLINCREMENT: the relative change in the norm of the parameter solutions is less than (a settable, setEpsValue(), default 1e-8) value.
- DERIVLEVEL: the inf-norm of the known vector of the equations to be solved is less than the settable, setEpsDerivative(), default 1e-8, value.
- MAXITER: maximum number of iterations reached (only possible if a number is explicitly set)
- NOREDUCTION: if the Levenberg-Marquardt correction factor goes towards infinity. I.e. if no Chi2 improvement seems possible. Have to redo the solution with a different start condition for the unknowns.
- SINGULAR: can only happen due to numeric rounding, since the LM equations are always positive-definite. Best solution is to indicate SVD needed in the solveLoop call, which is cost-free
- Covariance information in various forms can be obtained with the getCovariance(), getErrors(), getChi() (or getChi2), getSD and getWeightedSD methods after a solve() or after the final loop in a non-linear solution (of course, when necessary only).
An LSQFit object can be re-used by issuing the reset() command, or set() of new values. If an unknown has not been used in the condition equations at all, the doDiagonal() will make sure a proper solution is obtained, with missing unknowns zeroed.
Most of the calculations are done in place; however, enough data is saved that it is possible to continue with the same (partial) normal equations after e.g. an interim solution.
If the normal equations are produced in separate partial sets (e.g. in a multi-processor environment) a merge() method can combine them.
It is suggested to add any possible constraint equations after the merge.
A debugIt() method provides read access to all internal information.
The member definitions are split over three files. The second one contains the templated member function definitions, to bypass the problem of duplicate definitions of non-templated members when pre-compiling them. The third contains methods for saving objects as Records or through AipsIO.
The contents can be saved in a record (toRecord<src>),
and an object can be created from a record (<src>fromRecord).
The record identifier is 'lfit'.
The object can also be saved or restored using AipsIO.
Example
See the tLSQFit.cc and tLSQaips.cc program for extensive examples.The following example will first create 2 condition equations for 3 unknowns (the third is degenerate). It will first create normal equations for a 2 unknown solution and solve; then it will create normal equations for a 3 unknown solution, and solve (note that the degenerate will be set to 0. The last one will use SVD and one condition equation.r
#include <casa/aips.h>
#include <scimath/Fitting/LSQFit.h>
#include <iostream>
int main() {
// Condition equations for x+y=2; x-y=4;
Double ce[2][3] = {{1, 1, 0}, {1, -1, 0}};
Double m[2] = {2, 4};
// Solution and error area
Double sol[3];
Double sd, mu;
uInt rank;
Bool ok;
// LSQ object
LSQFit fit(2);
// Make normal equation
for (uInt i=0; i<2; i++) fit.makeNorm(ce[i], 1.0, m[i]);
// Invert(decompose) and show
ok = fit.invert(rank);
cout << "ok? " << ok << "; rank: " << rank << endl;
// Solve and show
if (ok) {
fit.solve(sol, &sd, &mu);
for (uInt i=0; i<2; i++) cout << "Sol" << i << ": " << sol[i] << endl;
cout << "sd: "<< sd << "; mu: " << mu << endl;
};
cout << "----------" << endl;
// Retry with 3 unknowns: note auto fill of unmentioned one
fit.set(uInt(3));
for (uInt i=0; i<2; i++) fit.makeNorm(ce[i], 1.0, m[i]);
ok = fit.invert(rank);
cout << "ok? " << ok << "; rank: " << rank << endl;
if (ok) {
fit.solve(sol, &sd, &mu);
for (uInt i=0; i<3; i++) cout << "Sol" << i << ": " << sol[i] << endl;
cout << "sd: "<< sd << "; mu: " << mu << endl;
};
cout << "----------" << endl;
// Retry with 3 unknowns; but 1 condition equation and use SVD
fit.reset();
for (uInt i=0; i<1; i++) fit.makeNorm(ce[i], 1.0, m[i]);
ok = fit.invert(rank, True);
cout << "ok? " << ok << "; rank: " << rank << endl;
if (ok) {
fit.solve(sol, &sd, &mu);
for (uInt i=0; i<3; i++) cout << "Sol" << i << ": " << sol[i] << endl;
cout << "sd: "<< sd << "; mu: " << mu << endl;
};
cout << "----------" << endl;
// Without SVD it would be:
fit.reset();
for (uInt i=0; i<1; i++) fit.makeNorm(ce[i], 1.0, m[i]);
ok = fit.invert(rank);
cout << "ok? " << ok << "; rank: " << rank << endl;
if (ok) {
fit.solve(sol, &sd, &mu);
for (uInt i=0; i<3; i++) cout << "Sol" << i << ": " << sol[i] << endl;
cout << "sd: "<< sd << "; mu: " << mu << endl;
};
cout << "----------" << endl;
exit(0);
}
Which will produce the output:
ok? 1; rank: 2 Sol0: 3 Sol1: -1 sd: 0; mu: 0 ---------- ok? 1; rank: 3 Sol0: 3 Sol1: -1 Sol2: 0 sd: 0; mu: 0 ---------- ok? 1; rank: 2 Sol0: 1 Sol1: 1 Sol2: 0 sd: 0; mu: 0 ---------- ok? 0; rank: 2 ----------
Motivation
The class was written to be able to do complex, real standard and SVD solutions in a simple and fast way.
To Do
- a thorough check if all loops are optimal in the makeNorm() methods
- input of condition equations with cross covariance
Member Description
enum ReadyCode
State of the non-linear solution
explicit LSQFit(uInt nUnknowns, uInt nConstraints=0)
Construct an object with the number of unknowns and constraints, using the default collinearity factor and the default Levenberg-Marquardt adjustment factor.Assume real
LSQFit(uInt nUnknowns, const LSQReal &, uInt nConstraints=0)
Construct an object with the number of unknowns and constraints, using the default collinearity factor and the default Levenberg-Marquardt adjustment factor.Allow explicit Real specification
LSQFit(uInt nUnknowns, const LSQComplex &, uInt nConstraints=0)
Construct an object with the number of unknowns and constraints, using the default collinearity factor and the default Levenberg-Marquardt adjustment factor.Allow explicit Complex specification
LSQFit()
Default constructor (empty, only usable after a set(nUnknowns))LSQFit(const LSQFit &other)
Copy constructor (deep copy)LSQFit &operator=(const LSQFit &other)
Assignment (deep copy)~LSQFit()
Bool invert(uInt &nRank, Bool doSVD=False)
Triangularize the normal equations and determine the rank nRank of the normal equations and, in the case of an SVD solution, the constraint equations. The collinearity factor is used to determine if the system can be solved (in essence it is the square of the sine of the angle between a column in the normal equations and the plane suspended by the other columns: if too parallel, the equations are degenerate). If doSVD is given as False, False is returned if rank not maximal, else an SVD solution is done.
template <class U> void copy(const Double *beg, const Double *end, U &sol, U)
template <class U> void copy(const Double *beg, const Double *end, U &sol, U)
template <class U> void copy(const Double *beg, const Double *end, U *sol, U)
template <class U> void copy(const Double *beg, const Double *end, U *sol, U)
template <class U> void uncopy(Double *beg, const Double *end, U &sol, U)
template <class U> void uncopy(Double *beg, const Double *end, U &sol, U)
template <class U> void uncopy(Double *beg, const Double *end, U *sol, U)
template <class U> void uncopy(Double *beg, const Double *end, U *sol, U)
template <class U> void copyDiagonal(U &errors, U)
template <class U> void copyDiagonal(U &errors, U)
Copy date from beg to end; converting if necessary to complex data
template <class U> void solve(U *sol)
template <class U> void solve(std::complex<U> *sol)
template <class U> void solve(U &sol)
Solve normal equations.
The solution will be given in sol.
template <class U> Bool solveLoop(uInt &nRank, U *sol, Bool doSVD=False)
template <class U> Bool solveLoop(uInt &nRank, std::complex<U> *sol, Bool doSVD=False)
template <class U> Bool solveLoop(uInt &nRank, U &sol, Bool doSVD=False)
template <class U> Bool solveLoop(Double &fit, uInt &nRank, U *sol, Bool doSVD=False)
template <class U> Bool solveLoop(Double &fit, uInt &nRank, std::complex<U> *sol, Bool doSVD=False)
template <class U> Bool solveLoop(Double &fit, uInt &nRank, U &sol, Bool doSVD=False)
Solve a loop in a non-linear set.
The methods with the fit argument are deprecated. Use
the combination without the 'fit' parameter, and the isReady()
call. The 'fit' parameter returns
for each loop a goodness
of fit indicator. If it is >0; more loops are necessary.
If it is negative,
and has an absolute value of say less than .001, it is probably ok, and
the iterations can be stopped.
Other arguments are as for solve() and invert().
The sol is used for both input (parameter guess) and output.
template <class U, class V> void makeNorm(const V &cEq, const U &weight, const U &obs, Bool doNorm=True, Bool doKnown=True)
template <class U, class V> void makeNorm(const V &cEq, const U &weight, const U &obs, LSQFit::Real, Bool doNorm=True, Bool doKnown=True)
template <class U, class V> void makeNorm(const V &cEq, const U &weight, const std::complex<U> &obs, Bool doNorm=True, Bool doKnown=True)
template <class U, class V> void makeNorm(const V &cEq, const U &weight, const std::complex<U> &obs, LSQFit::Complex, Bool doNorm=True, Bool doKnown=True)
template <class U, class V> void makeNorm(const V &cEq, const U &weight, const std::complex<U> &obs, LSQFit::Separable, Bool doNorm=True, Bool doKnown=True)
template <class U, class V> void makeNorm(const V &cEq, const U &weight, const std::complex<U> &obs, LSQFit::AsReal, Bool doNorm=True, Bool doKnown=True)
template <class U, class V> void makeNorm(const V &cEq, const U &weight, const std::complex<U> &obs, LSQFit::Conjugate, Bool doNorm=True, Bool doKnown=True)
template <class U, class V, class W> void makeNorm(uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const U &obs, Bool doNorm=True, Bool doKnown=True)
template <class U, class V, class W> void makeNorm(uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const U &obs, LSQFit::Real, Bool doNorm=True, Bool doKnown=True)
template <class U, class V, class W> void makeNorm(uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const std::complex<U> &obs, Bool doNorm=True, Bool doKnown=True)
template <class U, class V, class W> void makeNorm(uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const std::complex<U> &obs, LSQFit::Complex, Bool doNorm=True, Bool doKnown=True)
template <class U, class V, class W> void makeNorm(uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const std::complex<U> &obs, LSQFit::Separable, Bool doNorm=True, Bool doKnown=True)
template <class U, class V, class W> void makeNorm(uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const std::complex<U> &obs, LSQFit::AsReal, Bool doNorm=True, Bool doKnown=True)
template <class U, class V, class W> void makeNorm(uInt nIndex, const W &cEqIndex, const V &cEq, const U &weight, const std::complex<U> &obs, LSQFit::Conjugate, Bool doNorm=True, Bool doKnown=True)
Make normal equations using the cEq condition equation
(with nUnknowns elements) and a weight weight,
given the known observed value obs.
doNorm and doKnown can be used to e.g. re-use existing normal equations, i.e. the condition equations, but make a new known side (i.e. new observations).
The versions with cEqIndex[] indicate which of the nUnknowns are actually present in the condition equation (starting indexing at 0); the other terms are supposed to be zero. E.g. if a 12-telescope array has an equation only using telescopes 2 and 4, the lengths of cEqIndex and cEq will be both 2, and the index will contain 1 and 3 (when telescope numbering starts at 1) or 2 and 4 (when telescope numbering starts at 0. The index is given as an iterator (and hence can be a raw pointer)
Note that the use of const U & is due to a Float->Double conversion problem on Solaris. Linux was ok.
template <class U> Bool getConstraint(uInt n, U *cEq) const
template <class U> Bool getConstraint(uInt n, std::complex<U> *cEq) const
template <class U> Bool getConstraint(uInt n, U &cEq) const
Get the n-th (from 0 to the rank deficiency, or missing rank,
see e.g. getDeficiency())
constraint equation as determined by invert() in SVD-mode in
cEq[nUnknown]. False returned for illegal n. Note
that nMissing will be equal to the number of unknowns
(nUnknowns, or double that for the complex case) minus the
rank as returned from the invert() method.
template <class U, class V> Bool setConstraint(uInt n, const V &cEq, const U &obs)
template <class U, class V> Bool setConstraint(uInt n, const V &cEq, const std::complex<U> &obs)
template <class U, class V, class W> Bool setConstraint(uInt n, uInt nIndex, const W &cEqIndex, const V &cEq, const U &obs)
template <class U, class V, class W> Bool setConstraint(uInt n, uInt nIndex, const W &cEqIndex, const V &cEq, const std::complex<U> &obs)
template <class U, class V> Bool addConstraint(const V &cEq, const U &obs)
template <class U, class V> Bool addConstraint(const V &cEq, const std::complex<U> &obs)
template <class U, class V, class W> Bool addConstraint(uInt nIndex, const W &cEqIndex, const V &cEq, const U &obs)
template <class U, class V, class W> Bool addConstraint(uInt nIndex, const W &cEqIndex, const V &cEq, const std::complex<U> &obs)
Add a new constraint equation (updating nConstraints); or set a
numbered constraint equation (0..nConstraints-1). False if illegal
number n. The constraints are equations with nUnknowns terms,
and a constant value. E.g. measuring three angles of a triangle
could lead to equation [1,1,1] with obs as
3.1415. Note that each complex constraint will be
converted into two real constraints (see
Note 224).
Bool merge(const LSQFit &other)
Bool merge(const LSQFit &other, uInt nIndex, uInt *nEqIndex)
template <class W> Bool merge(const LSQFit &other, uInt nIndex, const W &nEqIndex)
Merge other LSQFit object (i.e. the normal equation and
related information) into this. Both objects must have the
same number of unknowns, and be pure normal equations (i.e. no
invert(), solve(), solveLoop() or statistics calls//
should have been made). If merging cannot be done, False
is returned. The index case (the index is an iterator) assumes that
the normal equations to be merged are a sparse subset of the complete
matrix. The index 'vector' specifies which unknown are present. An index
outside the scope of the final equations will be skipped.
void reset()
Reset status to empty
void set(uInt nUnknowns, uInt nConstraints=0)
void set(Int nUnknowns, Int nConstraints=0)
void set(uInt nUnknowns, const LSQReal &, uInt nConstraints=0)
void set(Int nUnknowns, const LSQReal &, Int nConstraints=0)
void set(uInt nUnknowns, const LSQComplex &, uInt nConstraints=0)
void set(Int nUnknowns, const LSQComplex &, Int nConstraints=0)
Set new sizes (default is for Real)
void set(Double factor=1e-6, Double LMFactor=1e-3)
Set new factors (collinearity factor, and Levenberg-Marquardt LMFactor)void setEpsValue(Double epsval=1e-8)
Set new value solution testvoid setEpsDerivative(Double epsder=1e-8)
Set new derivative testvoid setMaxIter(uInt maxiter=0)
Set maximum number of iterationsuInt nIterations()
Get number of iterations donevoid setBalanced(Bool balanced=False)
Set the expected form of the normal equations
LSQFit::ReadyCode isReady()
const std::string &readyText()
Ask the state of the non-linear solutions
template <class U> Bool getCovariance(U *covar)
template <class U> Bool getCovariance(std::complex<U> *covar)
Get the covariance matrix (of size nUnknowns * nUnknowns)
template <class U> Bool getErrors(U *errors)
template <class U> Bool getErrors(std::complex<U> *errors)
template <class U> Bool getErrors(U &errors)
Get main diagonal of covariance function (of size nUnknowns)
uInt nUnknowns() const
Get the number of unknownsuInt nConstraints() const
Get the number of constraintsuInt getDeficiency() const
Get the rank deficiency
Double getChi() const
Double getChi2() const
Double getSD() const
Double getWeightedSD() const
Get chi^2 (both are identical); the standard deviation (per observation)
and the standard deviation per weight unit.
void debugIt(uInt &nun, uInt &np, uInt &ncon, uInt &ner, uInt &rank, Double *&nEq, Double *&known, Double *&constr, Double *&er, uInt *&piv, Double *&sEq, Double *&sol, Double &prec, Double &nonlin) const
Debug:- nun = number of unknowns
- np = total number of solved unknowns (nun+ncon)
- ncon = number of constraint equations
- ner = number of elements in chi2 vector
- rank = rank)
- nEq = normal equation (nun*nun as triangular matrix)
- known = known vector (np)
- constr = constraint matrix (ncon*nun)
- er = error info vector (ner)
- piv = pivot vector (np)
- sEq = normal solution equation (np*np triangular)
- sol = internal solution vector (np)
- prec = collinearity precision
- nonlin = current Levenberg factor-1
Bool fromRecord(String &error, const RecordInterface &in)
Create an LSQFit object from a record. An error message is generated, and False returned if an invalid record is given. A valid record will return True. Error messages are postfixed to error.
Bool toRecord(String &error, RecordInterface &out) const
Create a record from an LSQFit object. The return will be False and an error message generated only if the object does not contain a valid object. Error messages are postfixed to error.const String &ident() const
Get identification of recordenum ErrorField
Offset of fields in error_p data area.enum StateBit
Bits that can be set/referenced
void toAipsIO (AipsIO&) const
void fromAipsIO (AipsIO&)
Work areas for interim solutions and covariance
Double *rowrt(uInt i) const
Double *rowru(uInt i) const
Get pointer in rectangular array