This solves a nonlinear least squares problem using Levenberg-Marquardt. More...

#include <connor_solver_map.h>

Inheritance diagram for FullPhysics::ConnorSolverMAP:

Collaboration diagram for FullPhysics::ConnorSolverMAP:

Public Types
enum	status_t { SUCCESS, CONTINUE, STALLED, ERROR, UNTRIED }
	Enum type for the status of the iterative solver. More...

Public Member Functions
	ConnorSolverMAP (const boost::shared_ptr< NLLSMaxAPosteriori > &NLLS_MAP, const boost::shared_ptr< ConvergenceCheck > &Convergence_check, int max_cost_function_calls, bool vrbs=false, double Gamma_initial=0.0, const std::string &Fname_test_data="")
	and optionally an initial value for gamma. More...

virtual	~ConnorSolverMAP ()

virtual std::vector< blitz::Array< double, 1 > >	accepted_points () const
	Returns a vector (std) of accepted points. More...

virtual void	add_observer (Observer< IterativeSolver > &Obs)
	Add an observer. More...

void	add_observer_and_keep_reference (boost::shared_ptr< Observer< IterativeSolver > > &Obs)
	Add an observer and keep a reference to it. More...

boost::shared_ptr< ConvergenceCheck >	convergence_check () const
	The convergence check object. More...

virtual std::vector< double >	cost_at_accepted_points () const
	Returns a vector (std) of cost function values at accepted points. More...

virtual FitStatistic	fit_statistic () const
	Return fit results for solution to last problem solved. More...

double	gamma_last_step () const
	Levenberg-Marquardt parameter for last step we processed. More...

virtual std::vector< blitz::Array< double, 1 > >	gradient_at_accepted_points () const
	Returns a vector (std) of gradients evaluated at accepted points. More...

virtual int	num_accepted_steps () const
	Returns the number of the accepted steps. More...

int	number_divergent () const
	Number of divergent steps for the last problem solved. More...

int	number_iteration () const
	Number of iterations for the last problem solved. More...

int	outcome_flag () const
	Outcome flag. This is an integer version of FitStatistic::OUTCOME. More...

void	print (std::ostream &Os) const
	Prints description of object. More...

std::string	print_to_string () const
	Print to string. More...

virtual const boost::shared_ptr< NLLSProblem >	problem () const

virtual void	remove_observer (Observer< IterativeSolver > &Obs)
	Remove an observer. More...

void	solve ()
	This solves the least squares problem starting at the initial guess. More...

virtual status_t	status () const
	Returns a value of IterativeSolver::status_t type. More...

virtual const char *const	status_str () const
	Returns the string version of the solver status. More...

blitz::Array< double, 1 >	x_update () const
	Return the a priori of the last problem solved. More...

Protected Member Functions
void	add_observer_do (Observer< IterativeSolver > &Obs, IterativeSolver &t)
	Add an observer. More...

void	add_observer_do (Observer< IterativeSolver > &Obs)

void	clean_dead_ptr ()
	Remove any dead pointers. More...

void	do_inversion ()
	This does an inversion step. More...

void	notify_update_do (const IterativeSolver &Self)
	Function to call to notify Observers of a state change. More...

void	record_accepted_point (const blitz::Array< double, 1 > &point)
	Called to record an accepted point. More...

void	record_cost_at_accepted_point (double cost)
	Called to record the cost function value at an accepted point. More...

void	record_gradient_at_accepted_point (const blitz::Array< double, 1 > &gradient)
	For recording the gradient of the cost function evaluated at an accepted point. More...

void	remove_observer_do (Observer< IterativeSolver > &Obs, IterativeSolver &t)
	Remove an observer. More...

void	remove_observer_do (Observer< IterativeSolver > &Obs)

Protected Attributes
boost::shared_ptr< ConvergenceCheck >	convergence_check_
	The convergence check object. More...

blitz::Array< double, 1 >	dx
	This is $d{x}_{i + 1}$ the update to $\mathbf{x}_i$ , $\mathbf{y} - \mathbf{F}(\mathbf{x}_i)$ . More...

std::string	fname_test_data
	If this isn't an empty string, save to this file in the first iteration. More...

FitStatistic	fstat
	Results from last fit step. More...

double	gamma
	Levenberg-Marquardt $\gamma$ parameter. More...

double	gamma_initial
	Initial value of gamma. More...

double	gamma_last_step_
	Stash gamma from last step we processed. More...

const boost::shared_ptr< MaxAPosteriori >	map

int	max_cost_f_calls

std::list< boost::weak_ptr< Observer< IterativeSolver > > >	olist

boost::shared_ptr< NLLSProblem >	P

std::vector< boost::shared_ptr< Observer< IterativeSolver > > >	ref_list

boost::shared_ptr< NLLSProblemScaled >	scaled_p

status_t	stat

bool	verbose

Static Protected Attributes
static double	rcond = 1e-12
	Factor to determine if we treat a singular factor as 0. More...

Detailed Description

This solves a nonlinear least squares problem using Levenberg-Marquardt.

This is the code that was originally written by Brian Connor.

We've rewritten the code in C++. This in entirely a matter a convenience, the wrapper code for calling the older Fortran code was bigger than the routine we were calling, so it made more sense just to duplicate the algorithm. We can revisit this if there is a reason to move this back to Fortran.

This solver is able to handle rank deficient Jacobians (e.g., we have parameters that are ignored).

At each iteration, this solves

$\left((1+\gamma) \mathbf{S}_a^{-1} + \mathbf{K}_i^T \mathbf{S}_{\epsilon}^{-1} \mathbf{K}_i\right) d\mathbf{x}_{i+1} = \left[\mathbf{K}_i^T \mathbf{S}_{\epsilon}^{-1} \left(\mathbf{y} - \mathbf{F}(\mathbf{x}_i)\right) + \mathbf{S}_a^{-1}(\mathbf{x}_a - \mathbf{x}_i)\right]$

This class has a "do_inversion" private member that is a bit difficult to test. As an aid to testing, the constructor takes a optional file name. If that name is passed, then we dump the state of this class in the first iteration of the solver to that file. We can then rerun the do_inversion by passing this file name to "test_do_inversion", which sets up the state, calls, do_inversion, and returns the results. This is entirely for use in unit testing, you wouldn't call this in normal operation.

This class is an Observable, any registered Observer will get notified after each iteration of the solver. This can be used to do things like add iteration output or extra logging.

Definition at line 52 of file connor_solver_map.h.

Member Enumeration Documentation

◆ status_t

enum FullPhysics::IterativeSolver::status_t

inherited

Enum type for the status of the iterative solver.

Enumerator
SUCCESS	solve method called and a solution found
CONTINUE	solve method called but did not converge to a solution
STALLED	solve method was called but stalled
ERROR	solve method called but an error was encountered
UNTRIED	solve method not called yet

Definition at line 50 of file iterative_solver.h.

Constructor & Destructor Documentation

◆ ConnorSolverMAP()

FullPhysics::ConnorSolverMAP::ConnorSolverMAP	(	const boost::shared_ptr< NLLSMaxAPosteriori > &	NLLS_MAP,
		const boost::shared_ptr< ConvergenceCheck > &	Convergence_check,
		int	max_cost_function_calls,
		bool	vrbs = `false`,
		double	Gamma_initial = `0.0`,
		const std::string &	Fname_test_data = `""`
	)

inline

and optionally an initial value for gamma.

See the comments above this class for the "save_test_data" argument.

Definition at line 59 of file connor_solver_map.h.

◆ ~ConnorSolverMAP()

virtual FullPhysics::ConnorSolverMAP::~ConnorSolverMAP ( )

inlinevirtual

Definition at line 69 of file connor_solver_map.h.

Member Function Documentation

◆ accepted_points()

virtual std::vector< blitz::Array<double, 1> > FullPhysics::IterativeSolver::accepted_points ( ) const

inlinevirtualinherited

Returns a vector (std) of accepted points.

This method returns a std vector of accepted points in the parameter space. The initial starting point is always an accepted point. Then the second accepted point is the point obtained after taking the first accepted step from the initial (first) point. The third accepted point is the point obtained after taking the second accepted step from the second accepted point and so on.

In other words, if the initial point and all the accepted points after taking the accepted steps are recorded correctly, then

accepted_points()[0] is the initial starting point,
accepted_points()[1] is the point obtained after taking the first accepted step from the initial point,
accepted_points()[2] is the point obtained after taking the second accepted step from accepted_points()[1] point,
...
accepted_points()[num_accepted_steps()] is the last accepted point obtained after the last accepted step.

Therefore, if the recording of the accepted points is done correctly, and num_accepted_steps() returns 2, then

accepted_points()[1] - accepted_points()[0] is the first accepted step, and
accepted_points()[2] - accepted_points()[1] is the second accepted step

Returns: A vector of accepted points

Definition at line 133 of file iterative_solver.h.

◆ add_observer()

virtual void FullPhysics::IterativeSolver::add_observer ( Observer< IterativeSolver > & Obs )

inlinevirtualinherited

Add an observer.

Implements FullPhysics::Observable< IterativeSolver >.

Definition at line 84 of file iterative_solver.h.

◆ add_observer_and_keep_reference()

void FullPhysics::Observable< IterativeSolver >::add_observer_and_keep_reference ( boost::shared_ptr< Observer< IterativeSolver > > & Obs )

inlineinherited

Add an observer and keep a reference to it.

See the discussion in the Observer class description for details.

Definition at line 107 of file observer.h.

◆ add_observer_do() [1/2]

void FullPhysics::Observable< IterativeSolver >::add_observer_do	(	Observer< IterativeSolver > &	Obs,
		IterativeSolver &	t
	)

inlineprotectedinherited

Add an observer.

Definition at line 148 of file observer.h.

◆ add_observer_do() [2/2]

void FullPhysics::Observable< IterativeSolver >::add_observer_do ( Observer< IterativeSolver > & Obs )

inlineprotectedinherited

Definition at line 159 of file observer.h.

◆ clean_dead_ptr()

void FullPhysics::Observable< IterativeSolver >::clean_dead_ptr ( )

inlineprotectedinherited

Remove any dead pointers.

Definition at line 196 of file observer.h.

◆ convergence_check()

boost::shared_ptr<ConvergenceCheck> FullPhysics::ConnorSolverMAP::convergence_check ( ) const

inline

The convergence check object.

Definition at line 85 of file connor_solver_map.h.

◆ cost_at_accepted_points()

virtual std::vector<double> FullPhysics::IterativeSolver::cost_at_accepted_points ( ) const

inlinevirtualinherited

Returns a vector (std) of cost function values at accepted points.

This method returns a std vector of cost function values computed at the accepted points. In other words, if the accepted points and the cost function values at these points are recorded correctly, then

cost_at_accepted_points()[0] is the value of the cost function at accepted_points()[0]
cost_at_accepted_points()[1] is the value of the cost function at accepted_points()[1]
...
and finally cost_at_accepted_points()[num_accepted_steps()] is the value of the cost function at accepted_points()[num_accepted_steps()]

Returns: A vector of cost function values at accepted points

Definition at line 157 of file iterative_solver.h.

◆ do_inversion()

void ConnorSolverMAP::do_inversion ( )

protected

This does an inversion step.

For this we solve the system

$\left((1+\gamma) \mathbf{S}_a^{-1} + \mathbf{K}_i^T \mathbf{S}_{\epsilon}^{-1} \mathbf{K}_i\right) d\mathbf{x}_{i+1} = \left[\mathbf{K}_i^T \mathbf{S}_{\epsilon}^{-1} \left(\mathbf{y} - \mathbf{F}(\mathbf{x}_i)\right) + \mathbf{S}_a^{-1}(\mathbf{x}_a - \mathbf{x}_i)\right]$

Note, to improve numerical stability we scale the system by

$\mathbf{N} = \sqrt{diag \mathbf{S}_a}$ and solve the system

$\mathbf{N}^T \left((1+\gamma) \mathbf{S}_a^{-1} + \mathbf{K}_i^T \mathbf{S}_{\epsilon}^{-1} \mathbf{K}_i\right) \mathbf{N} \left( \mathbf{N}^{-1} d\mathbf{x}_{i+1} \right) = \mathbf{N}^T \left[\mathbf{K}_i^T \mathbf{S}_{\epsilon}^{-1} \left(\mathbf{y} - \mathbf{F}(\mathbf{x}_i)\right) + \mathbf{S}_a^{-1}(\mathbf{x}_a - \mathbf{x}_i)\right]$

To help with testing for convergence, we also calculate:

$d \sigma^2 = \left(d\mathbf{x}_{i+1}\right)^T \left[\mathbf{K}_i^T \mathbf{S}_{\epsilon}^{-1} \left(\mathbf{y} - \mathbf{F}(\mathbf{x}_i)\right) + \mathbf{S}_a^{-1}(\mathbf{x}_a - \mathbf{x}_i)\right]$

This calculates the values kt_se_m1_k, dx, and fstat.

Definition at line 166 of file connor_solver_map.cc.

◆ fit_statistic()

virtual FitStatistic FullPhysics::ConnorSolverMAP::fit_statistic ( ) const

inlinevirtual

Return fit results for solution to last problem solved.

Definition at line 120 of file connor_solver_map.h.

◆ gamma_last_step()

double FullPhysics::ConnorSolverMAP::gamma_last_step ( ) const

inline

Levenberg-Marquardt parameter for last step we processed.

Definition at line 78 of file connor_solver_map.h.

◆ gradient_at_accepted_points()

virtual std::vector< blitz::Array<double, 1> > FullPhysics::IterativeSolverDer::gradient_at_accepted_points ( ) const

inlinevirtualinherited

Returns a vector (std) of gradients evaluated at accepted points.

This method returns a std vector of gradients computed at the accepted points. In other words, if the accepted points and the computed gradients at these points are recorded correctly, then

gradient_at_accepted_points()[0] is the gradient of the cost function evaluated at accepted_points()[0]
gradient_at_accepted_points()[1] is the gradient of the cost function evaluated at accepted_points()[1]
...
and finally gradient_at_accepted_points()[num_accepted_steps()] is the gradient of the cost function evaluated at accepted_points()[num_accepted_steps()]

Definition at line 62 of file iterative_solver_der.h.

◆ notify_update_do()

void FullPhysics::Observable< IterativeSolver >::notify_update_do ( const IterativeSolver & Self )

inlineprotectedinherited

Function to call to notify Observers of a state change.

The object should pass itself to this function, so it can be passed to the Observers.

Definition at line 121 of file observer.h.

◆ num_accepted_steps()

virtual int FullPhysics::IterativeSolver::num_accepted_steps ( ) const

inlinevirtualinherited

Returns the number of the accepted steps.

Returns: Number of the accepted steps

Definition at line 96 of file iterative_solver.h.

◆ number_divergent()

int FullPhysics::ConnorSolverMAP::number_divergent ( ) const

inline

Number of divergent steps for the last problem solved.

Definition at line 99 of file connor_solver_map.h.

◆ number_iteration()

int FullPhysics::ConnorSolverMAP::number_iteration ( ) const

inline

Number of iterations for the last problem solved.

Definition at line 92 of file connor_solver_map.h.

◆ outcome_flag()

int FullPhysics::ConnorSolverMAP::outcome_flag ( ) const

inline

Outcome flag. This is an integer version of FitStatistic::OUTCOME.

Definition at line 106 of file connor_solver_map.h.

◆ print()

void FullPhysics::ConnorSolverMAP::print ( std::ostream & Os ) const

inlinevirtual

Prints description of object.

Reimplemented from FullPhysics::NLLSSolver.

Definition at line 123 of file connor_solver_map.h.

◆ print_to_string()

std::string FullPhysics::Printable< IterativeSolver >::print_to_string ( ) const

inlineinherited

Print to string.

This is primarily useful for SWIG wrappers to this class, e.g. a to_s method in ruby.

Definition at line 31 of file printable.h.

◆ problem()

virtual const boost::shared_ptr<NLLSProblem> FullPhysics::NLLSSolver::problem ( ) const

inlinevirtualinherited

Definition at line 49 of file nlls_solver.h.

◆ record_accepted_point()

void FullPhysics::IterativeSolver::record_accepted_point ( const blitz::Array< double, 1 > & point )

inlineprotectedinherited

Called to record an accepted point.

This method is called to record an accepted point. It is the responsibility of the implementer of the solve() method to record the accepted points. The accepted points must be recorded in the same order that they are achieved.

Parameters

[in] point an accepted point in the parameter space

Definition at line 246 of file iterative_solver.h.

◆ record_cost_at_accepted_point()

void FullPhysics::IterativeSolver::record_cost_at_accepted_point ( double cost )

inlineprotectedinherited

Called to record the cost function value at an accepted point.

This method is called to record the cost function value at an accepted point. It is the responsibility of the implementer of the solve() method to record the cost function values at the accepted points. The cost values must be recorded in the same order that they are evaluated.

Parameters

[in] cost cost function value at an accepted point in the parameter space

Definition at line 269 of file iterative_solver.h.

◆ record_gradient_at_accepted_point()

void FullPhysics::IterativeSolverDer::record_gradient_at_accepted_point ( const blitz::Array< double, 1 > & gradient )

inlineprotectedinherited

For recording the gradient of the cost function evaluated at an accepted point.

This method is called to record the gradient of the cost function evaluated at an accepted point. It is the responsibility of the implementer of the solve() method to record the gradients evaluated at the accepted points. The gradients must be recorded in the same order that they are evaluated.

Parameters

[in] gradient gradient of the cost function evaluated at an accepted point in the parameter space

Definition at line 91 of file iterative_solver_der.h.

◆ remove_observer()

virtual void FullPhysics::IterativeSolver::remove_observer ( Observer< IterativeSolver > & Obs )

inlinevirtualinherited

Remove an observer.

Implements FullPhysics::Observable< IterativeSolver >.

Definition at line 87 of file iterative_solver.h.

◆ remove_observer_do() [1/2]

void FullPhysics::Observable< IterativeSolver >::remove_observer_do	(	Observer< IterativeSolver > &	Obs,
		IterativeSolver &	t
	)

inlineprotectedinherited

Remove an observer.

Definition at line 173 of file observer.h.

◆ remove_observer_do() [2/2]

void FullPhysics::Observable< IterativeSolver >::remove_observer_do ( Observer< IterativeSolver > & Obs )

inlineprotectedinherited

Definition at line 181 of file observer.h.

◆ solve()

void ConnorSolverMAP::solve ( )

virtual

This solves the least squares problem starting at the initial guess.

We don't directly return the solution and related matrices, you can query this object for them after the solver has completed.

Implements FullPhysics::IterativeSolver.

Definition at line 56 of file connor_solver_map.cc.

◆ status()

virtual status_t FullPhysics::IterativeSolver::status ( ) const

inlinevirtualinherited

Returns a value of IterativeSolver::status_t type.

This method returns the status of the solver. The status of the solver is initialized to IterativeSolver::UNTRIED, then it must be set to one of the following values by the implemented version of solve() method:

Please, read the comments on IterativeSolver::status_t type and its possible values.

Returns: Solver status

Definition at line 190 of file iterative_solver.h.

◆ status_str()

const char *const IterativeSolver::status_str ( ) const

virtualinherited

Returns the string version of the solver status.

If the method status() returns

then status_str() will return

"UNTRIED",
"SUCCESS",
"CONTINUE",
"STALLED", or
"ERROR"

respectively.

Returns: Solver status in string form

Definition at line 14 of file iterative_solver.cc.

◆ x_update()

blitz::Array<double, 1> FullPhysics::ConnorSolverMAP::x_update ( ) const

inline

Return the a priori of the last problem solved.

Definition at line 114 of file connor_solver_map.h.

Member Data Documentation

◆ convergence_check_

boost::shared_ptr<ConvergenceCheck> FullPhysics::ConnorSolverMAP::convergence_check_

protected

The convergence check object.

Definition at line 172 of file connor_solver_map.h.

◆ dx

blitz::Array<double, 1> FullPhysics::ConnorSolverMAP::dx

protected

This is $d{x}_{i + 1}$ the update to $\mathbf{x}_i$ , $\mathbf{y} - \mathbf{F}(\mathbf{x}_i)$ .

Definition at line 165 of file connor_solver_map.h.

◆ fname_test_data

std::string FullPhysics::ConnorSolverMAP::fname_test_data

protected

If this isn't an empty string, save to this file in the first iteration.

Definition at line 135 of file connor_solver_map.h.

◆ fstat

FitStatistic FullPhysics::ConnorSolverMAP::fstat

protected

Results from last fit step.

Definition at line 155 of file connor_solver_map.h.

◆ gamma

double FullPhysics::ConnorSolverMAP::gamma

protected

Levenberg-Marquardt $\gamma$ parameter.

This gets updated in each iteration in "solve"

Definition at line 143 of file connor_solver_map.h.

◆ gamma_initial

double FullPhysics::ConnorSolverMAP::gamma_initial

protected

Initial value of gamma.

Definition at line 177 of file connor_solver_map.h.

◆ gamma_last_step_

double FullPhysics::ConnorSolverMAP::gamma_last_step_

protected

Stash gamma from last step we processed.

Definition at line 149 of file connor_solver_map.h.

◆ map

const boost::shared_ptr<MaxAPosteriori> FullPhysics::ConnorSolverMAP::map

protected

Definition at line 182 of file connor_solver_map.h.

◆ max_cost_f_calls

int FullPhysics::IterativeSolver::max_cost_f_calls

protectedinherited

Definition at line 229 of file iterative_solver.h.

◆ olist

std::list<boost::weak_ptr<Observer<IterativeSolver > > > FullPhysics::Observable< IterativeSolver >::olist

protectedinherited

Definition at line 200 of file observer.h.

◆ P

boost::shared_ptr<NLLSProblem> FullPhysics::NLLSSolver::P

protectedinherited

Definition at line 64 of file nlls_solver.h.

◆ rcond

double ConnorSolverMAP::rcond = 1e-12

staticprotected

Factor to determine if we treat a singular factor as 0.

This is Rcond in solve_least_squares routine.

Definition at line 179 of file connor_solver_map.h.

◆ ref_list

std::vector<boost::shared_ptr<Observer<IterativeSolver > > > FullPhysics::Observable< IterativeSolver >::ref_list

protectedinherited

Definition at line 201 of file observer.h.

◆ scaled_p

boost::shared_ptr<NLLSProblemScaled> FullPhysics::ConnorSolverMAP::scaled_p

protected

Definition at line 183 of file connor_solver_map.h.

◆ stat

status_t FullPhysics::IterativeSolver::stat

protectedinherited

Definition at line 232 of file iterative_solver.h.

◆ verbose

bool FullPhysics::IterativeSolver::verbose

protectedinherited

Definition at line 231 of file iterative_solver.h.

The documentation for this class was generated from the following files:

Implementation/connor_solver_map.h
Implementation/connor_solver_map.cc

Public Types

Public Member Functions

Protected Member Functions

Protected Attributes

Static Protected Attributes

Detailed Description

Member Enumeration Documentation

◆ status_t

Constructor & Destructor Documentation

◆ ConnorSolverMAP()

◆ ~ConnorSolverMAP()

Member Function Documentation

◆ accepted_points()

◆ add_observer()

◆ add_observer_and_keep_reference()

◆ add_observer_do() [1/2]

◆ add_observer_do() [2/2]

◆ clean_dead_ptr()

◆ convergence_check()

◆ cost_at_accepted_points()

◆ do_inversion()

◆ fit_statistic()

◆ gamma_last_step()

◆ gradient_at_accepted_points()

◆ notify_update_do()

◆ num_accepted_steps()

◆ number_divergent()

◆ number_iteration()

◆ outcome_flag()

◆ print()

◆ print_to_string()

◆ problem()

◆ record_accepted_point()

◆ record_cost_at_accepted_point()

◆ record_gradient_at_accepted_point()

◆ remove_observer()

◆ remove_observer_do() [1/2]

◆ remove_observer_do() [2/2]

◆ solve()

◆ status()

◆ status_str()

◆ x_update()

Member Data Documentation

◆ convergence_check_

◆ dx

◆ fname_test_data

◆ fstat

◆ gamma

◆ gamma_initial

◆ gamma_last_step_

◆ map

◆ max_cost_f_calls

◆ olist

◆ P

◆ rcond

◆ ref_list

◆ scaled_p

◆ stat

◆ verbose