ARTS  2.3.1285(git:92a29ea9-dirty)
oem.cc File Reference

Optimal inversion methods for linear and non-linear inverse problems. More...

#include "arts.h"
#include <iostream>
#include "lin_alg.h"
#include "logic.h"
#include "math.h"
#include "stdlib.h"
#include "arts_omp.h"

Go to the source code of this file.

Functions

void separator (ostream &stream, Index length)
 
void log_init_li (ostream &stream)
 Initial log message, linear. More...
 
void log_step_li (ostream &stream, Index step_number, Numeric cost, Numeric cost_x, Numeric cost_y)
 Step log message, linear. More...
 
void log_finalize_li (ostream &stream)
 Final log message, linear. More...
 
void log_init_gn (ostream &stream, Numeric tol, Index max_iter)
 Initial log message, Gauss-Newton. More...
 
void log_step_gn (ostream &stream, Index step_number, Numeric cost, Numeric cost_x, Numeric cost_y, Numeric di2)
 Step log message, Gauss-Newton. More...
 
void log_finalize_gn (ostream &stream, bool converged, Index iter, Index max_iter)
 Final log message, Gauss-Newton. More...
 
void log_init_lm (ostream &stream, Numeric tol, Index max_iter)
 Initial log message, Levenberg-Marquardt. More...
 
void log_gamma_step_lm (ostream &stream, Numeric cost, Numeric cost_x, Numeric cost_y, Numeric gamma)
 
void log_step_lm (ostream &stream, Index step_number, Numeric cost, Numeric cost_x, Numeric cost_y, Numeric gamma, Numeric di2)
 Step log message, Levenberg-Marquardt. More...
 
void log_finalize_lm (ostream &stream, bool converged, Numeric cost, Numeric cost_x, Numeric cost_y, Numeric gamma, Numeric gamma_max, Index iter, Index max_iter)
 Final log message, Levenberg-Marquardt. More...
 
template<typename Se_t >
void oem_cost_y (Numeric &cost_y, ConstVectorView y, ConstVectorView yf, Se_t SeInv, const Numeric &normfac)
 Calculation of y-part of cost function. More...
 
void oem_cost_x (Numeric &cost_x, ConstVectorView x, ConstVectorView xa, ConstMatrixView SxInv, const Numeric &normfac)
 Calculation of x-part of cost function. More...
 
void mult_outer (MatrixView A, ConstMatrixView B, ConstVectorView b)
 Multiply matrix element-wise by outer product of vector. More...
 
void scale_columns (MatrixView A, ConstMatrixView B, ConstVectorView b)
 Scale columns. More...
 
void scale_rows (MatrixView A, ConstMatrixView B, ConstVectorView b)
 Scale rows. More...
 
Index oem_linear_nform (Vector &x, Matrix &G, Matrix &J, Vector &yf, Numeric &cost_y, Numeric &cost_x, ForwardModel &F, ConstVectorView xa, ConstVectorView x_norm, ConstVectorView y, ConstMatrixView SeInv, ConstMatrixView SxInv, const Numeric &cost_start, const bool &verbose)
 Linear OEM, n-form. More...
 
template<typename Se_t , typename Sx_t >
Index oem_gauss_newton (Vector &x, Matrix &G, Matrix &J, Vector &yf, Numeric &cost_y, Numeric &cost_x, Index &iter, ForwardModel &F, ConstVectorView xa, ConstVectorView x_norm, ConstVectorView y, const Se_t &SeInv, const Sx_t &SxInv, const Index max_iter, const Numeric tol, bool verbose)
 Gauss-Newton non-linear OEM using precomputed inverses, n-form. More...
 
template<typename Se_t , typename Sx_t >
Index oem_levenberg_marquardt (Vector &x, Matrix &G, Matrix &J, Vector &yf, Numeric &cost_y, Numeric &cost_x, Index &iter, ForwardModel &F, ConstVectorView xa, ConstVectorView x_norm, ConstVectorView y, const Se_t &SeInv, const Sx_t &SxInv, Index max_iter, Numeric tol, Numeric gamma_start, Numeric gamma_decrease, Numeric gamma_increase, Numeric gamma_max, Numeric gamma_threshold, bool verbose)
 Non-linear OEM using Levenberg-Marquardt method. More...
 
void oem_linear_mform (VectorView x, MatrixView G, ConstVectorView xa, ConstVectorView yf, ConstVectorView y, ConstMatrixView K, ConstMatrixView Se, ConstMatrixView Sx)
 Linear OEM, m-form. More...
 
bool oem_gauss_newton_m_form (VectorView x, ConstVectorView y, ConstVectorView xa, ForwardModel &K, ConstMatrixView Se, ConstMatrixView Sx, Numeric tol, Index max_iter)
 Non-linear OEM using Gauss-Newton method. More...
 
bool oem_gauss_newton_n_form (VectorView x, ConstVectorView y, ConstVectorView xa, ForwardModel &K, ConstMatrixView Se, ConstMatrixView Sx, Numeric tol, Index max_iter)
 Non-linear OEM using Gauss-Newton method. More...
 

Detailed Description

Optimal inversion methods for linear and non-linear inverse problems.

Author
Simon Pfreundschuh simon.nosp@m.pf@c.nosp@m.halme.nosp@m.rs.s.nosp@m.e
Date
Fri Apr 17 16:39:25 2015

Definition in file oem.cc.

Function Documentation

◆ log_finalize_gn()

void log_finalize_gn ( ostream &  stream,
bool  converged,
Index  iter,
Index  max_iter 
)

Final log message, Gauss-Newton.

Parameters
[in]streamStream to print message to.
[in]convergedconverged flag.
[in]costValue of cost function.
[in]cost_xx-component of cost function.
[in]cost_yy-component of cost function.
[in]iterFinal step number.
[in]max_iterMaximum step number.

Definition at line 163 of file oem.cc.

References separator().

◆ log_finalize_li()

void log_finalize_li ( ostream &  stream)

Final log message, linear.

Parameters
[in]streamStream to print message to.
[in]convergedconverged flag.
[in]costValue of cost function.
[in]cost_xx-component of cost function.
[in]cost_yy-component of cost function.
[in]iterFinal step number.
[in]max_iterMaximum step number.

Definition at line 91 of file oem.cc.

References separator().

Referenced by oem_linear_nform().

◆ log_finalize_lm()

void log_finalize_lm ( ostream &  stream,
bool  converged,
Numeric  cost,
Numeric  cost_x,
Numeric  cost_y,
Numeric  gamma,
Numeric  gamma_max,
Index  iter,
Index  max_iter 
)

Final log message, Levenberg-Marquardt.

Parameters
[in]streamStream to print message to.
[in]convergedconverged flag.
[in]costValue of cost function.
[in]cost_xx-component of cost function.
[in]cost_yy-component of cost function.
[in]gammaFinal gamma value.
[in]gamma_maxMaximum gamma value.
[in]iterFinal step number.
[in]max_iterMaximum step number.

Definition at line 267 of file oem.cc.

References separator().

◆ log_gamma_step_lm()

void log_gamma_step_lm ( ostream &  stream,
Numeric  cost,
Numeric  cost_x,
Numeric  cost_y,
Numeric  gamma 
)

Definition at line 207 of file oem.cc.

◆ log_init_gn()

void log_init_gn ( ostream &  stream,
Numeric  tol,
Index  max_iter 
)

Initial log message, Gauss-Newton.

Parameters
[in]streamStream to print message to.
[in]tolTolerance.
[in]max_iterMaximum Iterations.

Definition at line 104 of file oem.cc.

References separator().

◆ log_init_li()

void log_init_li ( ostream &  stream)

Initial log message, linear.

Parameters
[in]streamStream to print message to.

Definition at line 44 of file oem.cc.

References separator().

Referenced by oem_linear_nform().

◆ log_init_lm()

void log_init_lm ( ostream &  stream,
Numeric  tol,
Index  max_iter 
)

Initial log message, Levenberg-Marquardt.

Parameters
[in]streamStream to print message to.
[in]tolTolerance.
[in]max_iterMaximum Iterations.

Definition at line 187 of file oem.cc.

References separator().

◆ log_step_gn()

void log_step_gn ( ostream &  stream,
Index  step_number,
Numeric  cost,
Numeric  cost_x,
Numeric  cost_y,
Numeric  di2 
)

Step log message, Gauss-Newton.

Parameters
[in]streamStream to print message to.
[in]step_numberCurrent step number.
[in]costCurrent value of cost function.
[in]cost_xx-component of current cost function value.
[in]cost_yy-component of current cost function value.
[in]di2Convergence criterion.

Definition at line 133 of file oem.cc.

◆ log_step_li()

void log_step_li ( ostream &  stream,
Index  step_number,
Numeric  cost,
Numeric  cost_x,
Numeric  cost_y 
)

Step log message, linear.

Parameters
[in]streamStream to print message to.
[in]step_numberCurrent step number.
[in]costCurrent value of cost function.
[in]cost_xx-component of current cost function value.
[in]cost_yy-component of current cost function value.

Definition at line 66 of file oem.cc.

Referenced by oem_linear_nform().

◆ log_step_lm()

void log_step_lm ( ostream &  stream,
Index  step_number,
Numeric  cost,
Numeric  cost_x,
Numeric  cost_y,
Numeric  gamma,
Numeric  di2 
)

Step log message, Levenberg-Marquardt.

Parameters
[in]streamStream to print message to.
[in]step_numberCurrent step number.
[in]costCurrent value of cost function.
[in]cost_xx-component of current cost function value.
[in]cost_yy-component of current cost function value.
[in]gammaCurrent gamma value.
[in]di2Convergence criterion.

Definition at line 234 of file oem.cc.

◆ mult_outer()

void mult_outer ( MatrixView  A,
ConstMatrixView  B,
ConstVectorView  b 
)

Multiply matrix element-wise by outer product of vector.

Multiplies the matrix B by the outer product b' * b of b. This is used to scale the a priori covariance matrix to avoid numerical problems. The Matrices A and B can be the same but shouldn't overlap in other ways.

Parameters
[out]AThe matrix B divided by b'*b.
[in]BThe matrix B to divide.
[in]bThe vector used to form the outer product.

Definition at line 377 of file oem.cc.

References i, is_size(), n, and ConstVectorView::nelem().

Referenced by scale_rows().

◆ oem_cost_x()

void oem_cost_x ( Numeric cost_x,
ConstVectorView  x,
ConstVectorView  xa,
ConstMatrixView  SxInv,
const Numeric normfac 
)

Calculation of x-part of cost function.

This version is suitable if no term is already at hand.

Parameters
[out]cost_xConst function value
[in]xState vector.
[in]xaA priori state.
[in]SxInvInverse of relevenaty covariance matrix
[in]normfacNormalisation factor. The cost is scaled with 1/normfac
Author
Patrick Eriksson
Date
2015-10-05

Definition at line 345 of file oem.cc.

References dx, mult(), and ConstVectorView::nelem().

◆ oem_cost_y()

template<typename Se_t >
void oem_cost_y ( Numeric cost_y,
ConstVectorView  y,
ConstVectorView  yf,
Se_t  SeInv,
const Numeric normfac 
)

Calculation of y-part of cost function.

This version is suitable if no term is already at hand.

Parameters
[out]cost_yConst function value
[in]yMeasurement vector.
[in]yfFitted measurement.
[in]SeInvInverse of relevenaty covariance matrix
[in]normfacNormalisation factor. The cost is scaled with 1/normfac
Author
Patrick Eriksson
Date
2015-10-05

Definition at line 319 of file oem.cc.

References mult(), and ConstVectorView::nelem().

◆ oem_gauss_newton()

template<typename Se_t , typename Sx_t >
Index oem_gauss_newton ( Vector x,
Matrix G,
Matrix J,
Vector yf,
Numeric cost_y,
Numeric cost_x,
Index iter,
ForwardModel &  F,
ConstVectorView  xa,
ConstVectorView  x_norm,
ConstVectorView  y,
const Se_t &  SeInv,
const Sx_t &  SxInv,
const Index  max_iter,
const Numeric  tol,
bool  verbose 
)

Gauss-Newton non-linear OEM using precomputed inverses, n-form.

Computes the optimal nonlinear inverse of a given forward model using the Gauss-Newton method as given in eq. (5.8) in Rodgers book. The forward model is given by an object implementing the interface described by the ForwardModel class. Convergence is determined using equation (5.29). The given tolerance is scaled by the problem size n. If the method doesn't converge it abords after the given maximum number of iterations.

During execution two additional n-times-m and one n-times-n matrix is allocated. In addition to that, space for 4 length-n vectors and two length-m vectors is allocated. The given Matrix and Vector views may not overlap.

The returned Index value indicates if the computation was successful or if an error occured. The following encoding for errors is used:

0 - Success
1 - ERROR: Computation didn't converge. Maximum of iterations reached.
2 - ERROR: Computation didn't converge. Maximum gamma value reached.
9 - ERROR: Evaluation of forward model failed.
Parameters
[out]xThe optimal inverse x.
[out]GThe gain matrix corresponding to x.
[out]JThe jacobian as computed in the last LM iteration.
[out]yfThe fitted measurement vector.
[out]cost_yThe cost related to the measurement vector.
[out]cost_xThe cost related to the state vector.
[out]iterThe number of iterations.
[in]FThe ForwardModel representing the forward model to invert.
[in]xaThe size-n a-priori state vector.
[in]x_normNormalization vector for the state covariance matrix.
[in]yThe length-m, input measurement vector.
[in]SeInvThe inverse of the measurement error covariance (m,m)-matrix
[in]SxInvThe inverse of the a priori covariance (n,n)-matrix
[in]max_iterTha maximum number of iterations in case of no convergence.
[in]tolThe convergence criterium before scaling by the problem size.
[in]verboseIf true, log message are printed to stdout.
Returns
Error code indicating success or failure of the method.

Definition at line 1603 of file oem.cc.

References n, ConstMatrixView::ncols(), ConstVectorView::nelem(), and ConstMatrixView::nrows().

◆ oem_gauss_newton_m_form()

bool oem_gauss_newton_m_form ( VectorView  x,
ConstVectorView  y,
ConstVectorView  xa,
ForwardModel &  K,
ConstMatrixView  Se,
ConstMatrixView  Sx,
Numeric  tol,
Index  max_iter 
)

Non-linear OEM using Gauss-Newton method.

Computes the optimal nonlinear inverse of a given forward model using the Gauss-Newton method as given in eq. (5.10) in Rodgers book. The forward model is given by an object implementing the interface described by the ForwardModel class. Convergence is determined using equation (5.33). The given tolerance is scaled by the problem size m. If the method doesn't converge it abords after the given maximum number of iterations.

During the execution, space for up to 6 additional matrices and vectors is allocated. The given Matrix and Vector views should not overlap.

Parameters
[out]xThe optimal inverse x.
[in]yThe measurement vector containing m measurements.
[in]xaThe size-n a-priori state vector.
[in]KThe ForwardModel representing the forward model to invert.
[in]SeThe measurement error covariance (m,m)-matrix
[in]SxThe a priori covariance (n,n)-matrix
[in]tolThe convergence criterium before scaling by the problem size.
max_iterTha maximum number of iterations in case of no convergence.

Definition at line 1852 of file oem.cc.

References mult(), n, ConstVectorView::nelem(), solve(), and transpose().

◆ oem_gauss_newton_n_form()

bool oem_gauss_newton_n_form ( VectorView  x,
ConstVectorView  y,
ConstVectorView  xa,
ForwardModel &  K,
ConstMatrixView  Se,
ConstMatrixView  Sx,
Numeric  tol,
Index  max_iter 
)

Non-linear OEM using Gauss-Newton method.

Computes the optimal nonlinear inverse of a given forward model using the Gauss-Newton method as given in eq. (5.9) in Rodgers book. The forward model is given by an object implementing the interface described by the ForwardModel class. Convergence is determined using equation (5.29). The given tolerance is scaled by the problem size m. If the method doesn't converge it abords after the given maximum number of iterations.

During the execution, space for up to 6 additional matrices and vectors is allocated. The given Matrix and Vector views should not overlap.

Parameters
[out]xThe optimal inverse x.
[in]yThe measurement vector containing m measurements.
[in]xaThe size-n a-priori state vector.
[in]KThe ForwardModel representing the forward model to invert.
[in]SeThe measurement error covariance (m,m)-matrix
[in]SxThe a priori covariance (n,n)-matrix
[in]tolThe convergence criterium before scaling by the problem size.
max_iterTha maximum number of iterations in case of no convergence.
Returns
True if the method has converged, false otherwise.

Definition at line 1941 of file oem.cc.

References dx, inv(), mult(), n, ConstVectorView::nelem(), solve(), and transpose().

◆ oem_levenberg_marquardt()

template<typename Se_t , typename Sx_t >
Index oem_levenberg_marquardt ( Vector x,
Matrix G,
Matrix J,
Vector yf,
Numeric cost_y,
Numeric cost_x,
Index iter,
ForwardModel &  F,
ConstVectorView  xa,
ConstVectorView  x_norm,
ConstVectorView  y,
const Se_t &  SeInv,
const Sx_t &  SxInv,
Index  max_iter,
Numeric  tol,
Numeric  gamma_start,
Numeric  gamma_decrease,
Numeric  gamma_increase,
Numeric  gamma_max,
Numeric  gamma_threshold,
bool  verbose 
)

Non-linear OEM using Levenberg-Marquardt method.

Inverts a given non-linear forward model using the Levenberg-Marquardt method. The implementation follows Eq. (5.36) in Rodger's book.

Communication with the forward model is performed in the same way as in oem_gauss_newton().

The start value for gamma is given by the parameter gamma_start. If a new x value is found, gamma is decreased by a factor gamma_scale_dec. If the ost is increased, gamma is increased by a factor gamma_scale_inc. If gamma falls below gamma_threshold, it is set to zero. If no lower cost can be obtained with gamma = 0, gamma is set to 1. If gamma becomes larger than gamma_max and the cost can not be reduced, the iteration is stopped.

During the execution, space for two n-times-n and one n-times-m matrices is allocated as well as space for 5 length-n vectors and two length-m vectors.

The returned Index value indicates if the computation was successful or if an error occured. The following encoding for errors is used:

0 - Success
1 - ERROR: Computation didn't converge. Maximum of iterations reached.
2 - ERROR: Computation didn't converge. Maximum gamma value reached.
9 - ERROR: Evaluation of forward model failed.
Parameters
[out]xThe optimal estimator of the state vector.
[out]yfThe fitted state vector as computed in the second-last LM iteration.
[out]GThe gain matrix.
[out]JThe jacobian as computed in the second-last LM iteration.
[in]yThe size-m input measurement vector.
[in]xaThe size-n a priori state vector.
[in]KA forward model object implementing the FowardModel class.
[in]SeinvThe inverse of the covariance matrix of the measurement error.
[in]SxinvThe inverse of the covariance matrix of the a priori error.
[in]tolThe normalized convergence criterion.
[in]max_iterThe maximum number of iterations before abortion.
[in]gamma_startThe start value of the gamma factor.
[in]gamma_scale_decThe factor to decrease gamma by if the cost function was descreased.
[in]gamma_scale_incThe factor to increase gamma by if the cost function could not be decreased.
[in]gamma_maxThe maximum gamma value. If gamma == gamma_max and the cost function can not be decreased, the iteration is aborted.
[in]gamma_thresholdThe minimum value that gamma can take on before it is set to zero.
[in]verboseIf true, log messages are printed to stdout.
Returns
Error code indicating success or failure of the method.

Definition at line 1703 of file oem.cc.

References n, ConstMatrixView::ncols(), ConstVectorView::nelem(), and ConstMatrixView::nrows().

◆ oem_linear_mform()

void oem_linear_mform ( VectorView  x,
MatrixView  G,
ConstVectorView  xa,
ConstVectorView  yf,
ConstVectorView  y,
ConstMatrixView  K,
ConstMatrixView  Se,
ConstMatrixView  Sx 
)

Linear OEM, m-form.

Computes the inverse of a linear forward model by computing the MAP solution as described in Rodgers, Inverse Methods for Atmospheric Sounding, p. 67. This function uses the m-form ( Eq. (4.6) ) which requires the solution of a linear system of equations of size m-times-m.

For the execution 1 n-times-m matrices, 2 m-times-m matrices and a vector with m elements are allocated. The given Matrix and Vector views may not overlap.

Parameters
[out]xThe optimal, estimated state vector consisting of n elements.
[in]yThe measurement vector consisting of m elements.
[in,out]yfOn input yf should contain the value of the forward model at the linearization point. On output yf should contain the fitted measurement vector.
[in]xaThe mean a priori state vector.
[in]KThe weighting function (m,n)-matrix.
[in]SeThe measurement error covariance (m,m)-matrix.
[in]SxThe a priori covariance (n,n)-matrix.
[out]GThe gain matrix.

Definition at line 1789 of file oem.cc.

References inv(), mult(), n, ConstMatrixView::ncols(), ConstVectorView::nelem(), ConstMatrixView::nrows(), and transpose().

◆ oem_linear_nform()

Index oem_linear_nform ( Vector x,
Matrix G,
Matrix J,
Vector yf,
Numeric cost_y,
Numeric cost_x,
ForwardModel &  F,
ConstVectorView  xa,
ConstVectorView  x_norm,
ConstVectorView  y,
ConstMatrixView  SeInv,
ConstMatrixView  SxInv,
const Numeric cost_start,
const bool &  verbose 
)

Linear OEM, n-form.

Computes the inverse of a linear forward model by computing the MAP solution as described in Rodgers, Inverse Methods for Atmospheric Sounding, p. 67. This function uses the n-form (eq. (4.3)) which requires the solution of a linear system of equations of size n-times-n.

Requires the inverses of the covariance matrices for the state and measurement vector to be provided as arguments.

The returned Index value indicates if the computation was successful or if an error occured. The following encoding for errors is used:

0 - Success
9 - ERROR: Evaluation of forward model failed.
Returns
Convergence status, see oem_diagnostics
Parameters
[out]xThe optimal, estimated state vector consisting of n elements.
[out]GThe gain matrix.
[in]xaThe a priori state vector
[in]yfThe value of the forward model at a priori.
[in]yThe measurement vector consisting of m elements.
[in]KThe weighting function (m,n)-matrix
[in]SeInvThe inverse of the measurement error covariance (m,m)-matrix
[in]SxInvThe inverse of the a priori covariance (n,n)-matrix
Returns
Error code indicating success or failure of the method.

Definition at line 836 of file oem.cc.

References dx, F, J, log_finalize_li(), log_init_li(), log_step_li(), n, ConstMatrixView::ncols(), ConstVectorView::nelem(), and ConstMatrixView::nrows().

◆ scale_columns()

void scale_columns ( MatrixView  A,
ConstMatrixView  B,
ConstVectorView  b 
)

Scale columns.

Scales the columns of a matrix B by the elements of a vector b.

Parameters
AThe scaled matrix
BThe matrix to scale
bThe vector containing the scaling factors

Definition at line 406 of file oem.cc.

References i, is_size(), n, ConstMatrixView::ncols(), and ConstMatrixView::nrows().

◆ scale_rows()

void scale_rows ( MatrixView  A,
ConstMatrixView  B,
ConstVectorView  b 
)

Scale rows.

Scales the rows of a matrix B by the elements of a vector b.

Parameters
AThe scaled matrix
BThe matrix to scale
bThe vector containing the scaling factors

Definition at line 436 of file oem.cc.

References lm_hitran_2017::compute(), F, i, is_size(), J, joker, lubacksub(), ludcmp(), mult(), mult_outer(), n, ConstMatrixView::ncols(), ConstMatrixView::nrows(), and transpose().

◆ separator()

void separator ( ostream &  stream,
Index  length 
)