Generic least square fitting class. More...

#include <LeastSquares.h>

Collaboration diagram for Isis::LeastSquares:

Public Types
enum	SolveMethod { SVD, QRD, SPARSE }

Public Member Functions
	LeastSquares (Isis::BasisFunction &basis, bool sparse=false, int sparseRows=0, int sparseCols=0, bool jigsaw=false)
	Creates a LeastSquares Object.

	~LeastSquares ()
	Destroys the LeastSquares object.

void	AddKnown (const std::vector< double > &input, double expected, double weight=1.0)
	Invoke this method for each set of knowns.

std::vector< double >	GetInput (int row) const
	This method returns the data at the given row.

double	GetExpected (int row) const
	This method returns the expected value at the given row.

int	Rows () const
	This methods returns the number of rows in the matrix.

int	Solve (Isis::LeastSquares::SolveMethod method=SVD)
	After all the data has been registered through AddKnown, invoke this method to solve the system of equations.

double	Evaluate (const std::vector< double > &input)
	Invokes the BasisFunction Evaluate method.

std::vector< double >	Residuals () const
	Returns a vector of residuals (errors).

double	Residual (int i) const
	Returns the ith residual.

void	Weight (int index, double weight)
	Reset the weight for the ith known.

int	Knowns () const
	The number of knowns (or times AddKnown was invoked) linear combination of the variables.

double	GetSigma0 ()

int	GetDegreesOfFreedom ()

void	Reset ()

void	ResetSparse ()

bool	SparseErrorPropagation ()
	Error propagation for sparse least-squares solution.

std::vector< double >	GetEpsilons () const

void	SetParameterWeights (const std::vector< double > weights)

void	SetNumberOfConstrainedParameters (int n)

const gmm::row_matrix < gmm::rsvector< double > > &	GetCovarianceMatrix () const

Private Member Functions
void	SolveSVD ()
	After all the data has been registered through AddKnown, invoke this method to solve the system of equations.

void	SolveQRD ()
	After all the data has been registered through AddKnown, invoke this method to solve the system of equations with a QR decomposition of A = QR.

void	SolveCholesky ()

int	SolveSparse ()
	Solve using sparse class.

void	FillSparseA (const std::vector< double > &data)
	Invoke this method for each set of knowns for sparse solutions.

bool	ApplyParameterWeights ()

Private Attributes
std::vector< double >	p_xSparse
	sparse solution vector

std::vector< double >	p_epsilonsSparse
	sparse vector of total parameter corrections

std::vector< double >	p_parameterWeights
	vector of parameter weights

gmm::row_matrix< gmm::rsvector < double > >	p_sparseA
	design matrix 'A'

gmm::row_matrix< gmm::rsvector < double > >	p_normals
	normal equations matrix 'N'

gmm::dense_matrix< double >	p_ATb
	right-hand side vector

gmm::SuperLU_factor< double >	p_SLU_Factor
	decomposed normal equations matrix

bool	p_jigsaw

bool	p_sparse

bool	p_solved
	Boolean value indicating solution is complete.

int	p_currentFillRow

int	p_sparseRows

int	p_sparseCols

int	p_constrainedParameters
	constrained parameters

int	p_degreesOfFreedom
	degrees of freedom (redundancy)

double	p_sigma0
	sigma nought - reference variance

std::vector< std::vector < double > >	p_input
	A vector of the input variables to evaluate.

std::vector< double >	p_expected
	A vector of the expected values when solved.

std::vector< double >	p_sqrtWeight
	A vector of the square roots of the weights for each known value.

std::vector< double >	p_residuals
	A vector of the residuals (or difference between expected and solved values).

Isis::BasisFunction *	p_basis
	Pointer to the BasisFunction object.

Detailed Description

Generic least square fitting class.

This class can be used to solved systems of linear equations through least squares fitting. The solution is derived through singular value decomposition or QR decomposition. For example:

$x + y = 3$

$-2x + 3y = 1$

$2x - y = 2$

Is a simple system of equations that can be solved using this class.

Isis::BasisFunction b("Linear",2,2);
Isis::LeastSquares lsq(b);
std::vector<double> one;
one.push_back(1.0);
one.push_back(1.0);
lsq.AddKnown(one,3.0);
std::vector<double> two;
two.push_back(-2.0);
two.push_back(3.0);
lsq.AddKnown(two,1.0);
std::vector<double> tre;
tre.push_back(2.0);
tre.push_back(-1.0);
lsq.AddKnown(tre,2.0);
lsq.Solve();
double out1 = lsq.Evaluate(one);
double out2 = lsq.Evaluate(two);
double out3 = lsq.Evaluate(tre);

See Also: PolynomialUnivariate; PolynomialBiVariate; BasisFunction

Author: 2004-06-24 Jeff Anderson

History:

2004-06-24 Jeff Anderson, Original version

2007-11-16 Tracie Sucharski and Debbie A. Cook Added SolveQRD method for a faster solve

2008-02-06 Tracie Sucharski, Renamed Solve to SolveSVD and make private. Add public Solve, with a new parameter, method which defaults to SVD. This calls the correct solution method. This was done for documentation purposes-clarifies the default solves by SVD.

2008-04-16 Debbie Cook / Tracie Sucharski, Added SolveSparse.

2008-06-09 Tracie Sucharski, Added conditional compilations for Solaris. We could not get SuperLu to build under Solaris due to a confilict with Complex definitions.

2009-04-06 Tracie Sucharski, Added return value to SolveSparse, which indicates a column containing all zeroes which causes superlu to bomb.

2009-12-21 Jeannie Walldren, Changed p_weight to p_sqrtweight. Modified Weight() and AddKnown() to take the square root of the given weight and add it to the p_sqrtweight vector.

2010-04-20 Debbie A. Cook, Replaced SparseReset with Reset to reset all solution methods

2010-11-22 Debbie A. Cook - Merged with Ken Edmundson version

2013-12-29 Jeannie Backer - Improved error messages. Fixes #962.

Definition at line 115 of file LeastSquares.h.

Member Enumeration Documentation

enum Isis::LeastSquares::SolveMethod

Enumerator
SVD	Singular Value Decomposition.
QRD	QR Decomposition.
SPARSE	Sparse.

Definition at line 128 of file LeastSquares.h.

Constructor & Destructor Documentation

Isis::LeastSquares::LeastSquares	(	Isis::BasisFunction &	basis,
		bool	sparse = `false`,
		int	sparseRows = `0`,
		int	sparseCols = `0`,
		bool	jigsaw = `false`
	)

Creates a LeastSquares Object.

Parameters

basis A BasisFunction. This parameter allows for the least squares fitting to be applied to arbitrary equations.

Definition at line 41 of file LeastSquares.cpp.

References _FILEINFO_, p_ATb, p_basis, p_epsilonsSparse, p_normals, p_parameterWeights, p_sigma0, p_solved, p_sparseA, p_xSparse, and Isis::IException::Programmer.

Isis::LeastSquares::~LeastSquares ( )

Destroys the LeastSquares object.

Definition at line 83 of file LeastSquares.cpp.

Member Function Documentation

void Isis::LeastSquares::AddKnown	(	const std::vector< double > &	data,
		double	result,
		double	weight = `1.0`
	)

Invoke this method for each set of knowns.

Given our example in the description, we have three knowns and expecteds. They are

$(1,1) = 3$

$(-2,3) = 1$

$(2,-1) = 2$

Parameters

data	A vector of knowns.
result	The expected value for the knowns.
weight	(Default = 1.0) How strongly to weight this known. Weight less than 1 increases residual for this known, while weight greater than 1 decreases the residual for this known.

Exceptions

Isis::IException::Programmer - Number of elements in data does not match basis requirements

History:

2008-04-22 Tracie Sucharski, Fill sparse matrix.

2009-12-21 Jeannie Walldren, Modified code to add the square root of the weight to the vector p_sqrtweight.

Definition at line 117 of file LeastSquares.cpp.

References _FILEINFO_, FillSparseA(), Isis::BasisFunction::Name(), p_basis, p_expected, p_input, p_sqrtWeight, Isis::IException::Programmer, and Isis::BasisFunction::Variables().

Referenced by Isis::SurfaceModel::AddTriplet(), Isis::ReseauDistortionMap::SetFocalPlane(), Isis::VimsGroundMap::SetGround(), Isis::SpicePosition::SetPolynomial(), Isis::SpiceRotation::SetPolynomial(), Isis::VimsSkyMap::SetSky(), Isis::ReseauDistortionMap::SetUndistortedFocalPlane(), and Isis::Affine::Solve().

double Isis::LeastSquares::Evaluate ( const std::vector< double > & data )

Invokes the BasisFunction Evaluate method.

Parameters

data	The input variables to evaluate.

Returns: The evaluation for the input variable.

Exceptions

Isis::IException::Programmer - Unable to evaluate until a solution has been computed

Definition at line 732 of file LeastSquares.cpp.

References _FILEINFO_, Isis::BasisFunction::Evaluate(), p_basis, p_solved, and Isis::IException::Programmer.

Referenced by Isis::SurfaceModel::Evaluate(), Isis::ReseauDistortionMap::SetFocalPlane(), Isis::VimsGroundMap::SetGround(), Isis::VimsSkyMap::SetSky(), and Isis::ReseauDistortionMap::SetUndistortedFocalPlane().

void Isis::LeastSquares::FillSparseA ( const std::vector< double > & data )

private

Invoke this method for each set of knowns for sparse solutions.

The A sparse matrix must be filled as we go or we will quickly run out of memory for large solutions. So, expand the basis function, apply weights (which is done in the Solve method for the non-sparse case.

Parameters

data	A vector of knowns.

History:

2008-04-22 Tracie Sucharski - New method for sparse solutions.

2009-12-21 Jeannie Walldren - Changed variable name from p_weight to p_sqrtweight.

Definition at line 161 of file LeastSquares.cpp.

References Isis::BasisFunction::Expand(), p_basis, p_sparseA, p_sqrtWeight, and Isis::BasisFunction::Term().

Referenced by AddKnown().

double Isis::LeastSquares::GetExpected ( int row ) const

This method returns the expected value at the given row.

Parameters

row

Returns: double

Definition at line 199 of file LeastSquares.cpp.

References _FILEINFO_, p_expected, Isis::IException::Programmer, Rows(), and Isis::toString().

std::vector< double > Isis::LeastSquares::GetInput ( int row ) const

This method returns the data at the given row.

Parameters

row

Returns: std::vector<double>

Definition at line 184 of file LeastSquares.cpp.

References _FILEINFO_, p_input, Isis::IException::Programmer, Rows(), and Isis::toString().

int Isis::LeastSquares::Knowns ( ) const

inline

The number of knowns (or times AddKnown was invoked) linear combination of the variables.

Returns: int The number of knowns

Definition at line 145 of file LeastSquares.h.

References p_expected.

Referenced by Isis::VimsGroundMap::SetGround(), and Isis::VimsSkyMap::SetSky().

double Isis::LeastSquares::Residual ( int i ) const

Returns the ith residual.

That is, the difference between the evaluation of a known with the solution against the expected value. There is one residual for each time AddKnown was invoked.

Parameters

i	The number of times AddKnown was invoked to be evaluated.

Returns: The output value of the residual.

Exceptions

Isis::IException::Programmer - Unable to return residuals until a solution has been computed

Definition at line 769 of file LeastSquares.cpp.

References _FILEINFO_, p_residuals, p_solved, and Isis::IException::Programmer.

std::vector< double > Isis::LeastSquares::Residuals ( ) const

Returns a vector of residuals (errors).

That is, the difference between the evaluation of a known with the solution against the expected value.

Returns: The vector of residuals.

Exceptions

Isis::IException::Programmer - Unable to return residuals until a solution has been computed

Definition at line 749 of file LeastSquares.cpp.

References _FILEINFO_, p_residuals, p_solved, and Isis::IException::Programmer.

int Isis::LeastSquares::Rows ( ) const

This methods returns the number of rows in the matrix.

Returns: int

Definition at line 213 of file LeastSquares.cpp.

References p_input.

Referenced by GetExpected(), GetInput(), and Solve().

int Isis::LeastSquares::Solve ( Isis::LeastSquares::SolveMethod method = SVD )

After all the data has been registered through AddKnown, invoke this method to solve the system of equations.

You can then use the Evaluate and Residual methods freely.

History:

2008-04-16 Debbie Cook / Tracie Sucharski, Added SolveSparse.

2009-04-08 Tracie Sucharski - Added return value which will pass on what is returned from SolveSparse which is a column number of a column that contained all zeros.

2010-12-12 Debbie A. Cook Fixed "no data" test for SPARSE case

Definition at line 231 of file LeastSquares.cpp.

References _FILEINFO_, p_solved, QRD, Rows(), SolveQRD(), SolveSparse(), SolveSVD(), SPARSE, SVD, and Isis::IException::Unknown.

Referenced by Isis::ReseauDistortionMap::SetFocalPlane(), Isis::VimsGroundMap::SetGround(), Isis::SpicePosition::SetPolynomial(), Isis::SpiceRotation::SetPolynomial(), Isis::VimsSkyMap::SetSky(), Isis::ReseauDistortionMap::SetUndistortedFocalPlane(), Isis::SurfaceModel::Solve(), and Isis::Affine::Solve().

void Isis::LeastSquares::SolveQRD ( )

private

After all the data has been registered through AddKnown, invoke this method to solve the system of equations with a QR decomposition of A = QR.

You can then use the Evaluate and Residual methods freely. The QR decomposition is only slightly less reliable than the SVD, but much faster.

History:: 2009-12-21 Jeannie Walldren - Changed variable name from p_weight to p_sqrtweight.

Definition at line 371 of file LeastSquares.cpp.

References _FILEINFO_, Isis::BasisFunction::Coefficients(), Isis::BasisFunction::Evaluate(), Isis::BasisFunction::Expand(), p_basis, p_expected, p_input, p_residuals, p_solved, p_sqrtWeight, Isis::BasisFunction::SetCoefficients(), Isis::BasisFunction::Term(), and Isis::IException::Unknown.

Referenced by Solve().

int Isis::LeastSquares::SolveSparse ( )

private

Solve using sparse class.

After all the data has been registered through AddKnown, invoke this method to solve the system of equations Nx = b, where

N = ATPA b = ATPl, and

N is the "normal equations matrix; A is the so-called "design" matrix; P is the observation weight matrix (typically diagonal); l is the "computed - observed" column vector;

The solution is achieved using a sparse matrix formulation of the LU decomposition of the normal equations.

You can then use the Evaluate and Residual methods freely.

History:

2008-04-16 Debbie Cook / Tracie Sucharski, New method

2008-04-23 Tracie Sucharski, Fill sparse matrix as we go in AddKnown method rather than in the solve method, otherwise we run out of memory very quickly.

2009-04-08 Tracie Sucharski - Added return value which is a column number of a column that contained all zeros.

2009-12-21 Jeannie Walldren - Changed variable name from p_weight to p_sqrtweight.

2010 Ken Edmundson

2010-11-20 Debbie A. Cook Merged Ken Edmundson verion with system version

2011-03-17 Ken Edmundson Corrected computation of residuals

Definition at line 461 of file LeastSquares.cpp.

References p_ATb, p_basis, p_constrainedParameters, p_degreesOfFreedom, p_epsilonsSparse, p_expected, p_normals, p_parameterWeights, p_residuals, p_sigma0, p_SLU_Factor, p_solved, p_sparseA, p_sqrtWeight, p_xSparse, and Isis::BasisFunction::SetCoefficients().

Referenced by Solve().

void Isis::LeastSquares::SolveSVD ( )

private

After all the data has been registered through AddKnown, invoke this method to solve the system of equations.

You can then use the Evaluate and Residual methods freely.

History:: 2009-12-21 Jeannie Walldren - Changed variable name from p_weight to p_sqrtweight.

Definition at line 268 of file LeastSquares.cpp.

References _FILEINFO_, Isis::BasisFunction::Coefficients(), Isis::BasisFunction::Evaluate(), Isis::BasisFunction::Expand(), p_basis, p_degreesOfFreedom, p_expected, p_input, p_residuals, p_sigma0, p_solved, p_sqrtWeight, Isis::BasisFunction::SetCoefficients(), Isis::BasisFunction::Term(), Isis::toString(), and Isis::IException::Unknown.

Referenced by Solve().

bool Isis::LeastSquares::SparseErrorPropagation ( )

Error propagation for sparse least-squares solution.

Computes the variance-covariance matrix of the parameters. This is the inverse of the normal equations matrix, scaled by the reference variance (also called variance of unit weight, etc).

History:: 2009-11-19 Ken Edmundson, New method

Notes:

1) The SLU_Factor (Super LU Factor) has already been factorised in each iteration but there is no gmm++ method to get the inverse of the sparse matrix implementation. so we have to get it ourselves. This is don by solving the factorized normals repeatedly with right-hand sides that are the columns of the identity matrix (which is how gmm - or anybody else would do it).

2) We should create our own matrix library, probably wrapping the gmm++ library (and perhaps other(s) that may have additional desired functionality). The inverse function should be part of the library, along with the capacity for triangular storage and other decomposition techniques - notably Cholesky.

3) The LU decomposition can be be performed even if the normal equations are singular. But, we should always be dealing with a positive-definite matrix (for the bundle adjustment). Cholesky is faster, more efficient, and requires a positive-definite matrix, so if it fails, it is an indication of a singular matrix - bottom line - we should be using Cholesky. Or a derivative of Cholesky, i.e., UTDU (LDLT).

4) As a consequence of 3), we should be checking for positive-definite state of the normals, perhaps via the matrix determinant, prior to solving. There is a check in place now that checks to see if a column of the design matrix (or basis?) is all zero. This is equivalent - if a set of vectors contains the zero vector, then the set is linearly dependent, and the matrix is not positive-definite. In Jigsaw, the most likely cause of the normals being non-positive-definite probably is failure to establish the datum (i.e. constraining a minimum of seven parameters - usually 3 coordinates of two points and 1 of a third).

Definition at line 649 of file LeastSquares.cpp.

References p_ATb, p_normals, p_sigma0, p_SLU_Factor, and p_xSparse.

void Isis::LeastSquares::Weight	(	int	index,
		double	weight
	)

Reset the weight for the ith known.

This weight will not be used unless the system is resolved using the Solve method.

Parameters

index	The position in the array to assign the given weight value
weight	A weight factor to apply to the ith known. A weight less than one increase the residual for this known while a weight greater than one decrease the residual for this known.