Main API¶

dogleg_optimize2

This is the main call to the library for sparse Jacobians. It's declared as

 double dogleg_optimize2(double* p, unsigned int Nstate,
                         unsigned int Nmeas, unsigned int NJnnz,
                         dogleg_callback_t* f, void* cookie,
                         const dogleg_parameters2_t* parameters,
                         dogleg_solverContext_t** returnContext);

• p is the initial estimate of the state vector (and holds the final result)
• "Nstate" specifies the number of optimization variables (elements of p)
• "Nmeas" specifies the number of measurements (elements of x). "Nmeas >= Nstate" is a requirement
• "NJnnz" specifies the number of non-zero elements of the jacobian matrix df/dp. In a dense matrix "Jnnz = Nstate*Nmeas". We are dealing with sparse jacobians, so "NJnnz" should be far less. If not, libdogleg is not an appropriate routine to solve this problem.
• "f" specifies the callback function that the optimization routine calls to sample the problem being solved
• "cookie" is an arbitrary data pointer passed untouched to "f"
• "parameters" a pointer to the set of parameters to use. Set to "NULL" to use the global parameters, or call "dogleg_optimize" instead. See the Parameters section below for more details
• If not "NULL", "returnContext" can be used to retrieve the full context structure from the solver. This can be useful since this structure contains the latest operating point values. It also has an active "cholmod_common" structure, which can be reused if more CHOLMOD routines need to be called externally. You usually want "returnContext->beforeStep". If this data is requested, the user is required to free it with "dogleg_freeContext" when done.

"dogleg_optimize" returns norm2( x ) at the minimum, or a negative value if an error occurred.

dogleg_optimize

This is a flavor of "dogleg_optimize2" that implicitly uses the global parameters. It's declared as

 double dogleg_optimize(double* p, unsigned int Nstate,
                        unsigned int Nmeas, unsigned int NJnnz,
                        dogleg_callback_t* f, void* cookie,
                        dogleg_solverContext_t** returnContext);

dogleg_optimize_dense2

This is the main call to the library for dense Jacobians. Its declared as

 double dogleg_optimize_dense2(double* p, unsigned int Nstate,
                               unsigned int Nmeas,
                               dogleg_callback_dense_t* f, void* cookie,
                               const dogleg_parameters2_t* parameters,
                               dogleg_solverContext_t** returnContext);

The arguments are almost identical to those in the "dogleg_optimize" call.

• p is the initial estimate of the state vector (and holds the final result)
• "Nstate" specifies the number of optimization variables (elements of p)
• "Nmeas" specifies the number of measurements (elements of x). "Nmeas >= Nstate" is a requirement
• "f" specifies the callback function that the optimization routine calls to sample the problem being solved. Note that this callback has a different type from that in "dogleg_optimize"
• "cookie" is an arbitrary data pointer passed untouched to "f"
• "parameters" a pointer to the set of parameters to use. Set to "NULL" to use the global parameters, or call "dogleg_optimize" instead. See the Parameters section below for more details
• If not "NULL", "returnContext" can be used to retrieve the full context structure from the solver. This can be useful since this structure contains the latest operating point values. You usually want "returnContext->beforeStep". If this data is requested, the user is required to free it with "dogleg_freeContext" when done.

"dogleg_optimize" returns norm2( x ) at the minimum, or a negative value if an error occurred.

dogleg_optimize_dense

This is a flavor of "dogleg_optimize_dense2" that implicitly uses the global parameters. It's declared as

 double dogleg_optimize_dense(double* p, unsigned int Nstate,
                              unsigned int Nmeas,
                              dogleg_callback_dense_t* f, void* cookie,
                              dogleg_solverContext_t** returnContext);

dogleg_freeContext

Used to deallocate memory used for an optimization cycle. Defined as:

 void dogleg_freeContext(dogleg_solverContext_t** ctx);

If a pointer to a context is not requested (by passing "returnContext = NULL" to "dogleg_optimize"), libdogleg calls this routine automatically. If the user did retrieve this pointer, though, it must be freed with "dogleg_freeContext" when the user is finished.

dogleg_computeJtJfactorization

Computes the cholesky decomposition of JtJ. This function is only exposed if you need to touch libdogleg internals via returnContext. Sometimes after computing an optimization you want to do stuff with the factorization of JtJ, and this call ensures that the factorization is available. Most people don't need this function. If the comment wasn't clear, you don't need this function.

This is declared as

 void dogleg_computeJtJfactorization(dogleg_operatingPoint_t* point,
                                     dogleg_solverContext_t* ctx);

The arguments are

• "point" is the operating point of the system. Generally this will be "returnContext->beforeStep" where "returnContext" is from one of the "dogleg_optimize_..." functions.
• "ctx" is the dogleg context. Generally this will be "returnContext" from one of the "dogleg_optimize_..." functions

dogleg_testGradient

libdogleg requires the user to compute the jacobian matrix J. This is a performance optimization, since J could be computed by differences of x. This optimization is often worth the extra effort, but it creates a possibility that J will have a mistake and J = df/dp would not be true. To find these types of issues, the user can call

 void dogleg_testGradient(unsigned int var, const double* p0,
                          unsigned int Nstate, unsigned int Nmeas, unsigned int NJnnz,
                          dogleg_callback_t* f, void* cookie);

This function computes the jacobian with center differences and compares the result with the jacobian computed by the callback function. It is recommended to do this for every variable while developing the program that uses libdogleg.

• "var" is the index of the variable being tested
• "p0" is the state vector p where we're evaluating the jacobian
• "Nstate", "Nmeas", "NJnnz" are the number of state variables, measurements and non-zero jacobian elements, as before
• "f" is the callback function, as before
• "cookie" is the user data, as before

This function returns nothing, but prints out the test results.

dogleg_testGradient_dense

Very similar to "dogleg_testGradient", but for dense jacobians.

 void dogleg_testGradient_dense(unsigned int var, const double* p0,
                                unsigned int Nstate, unsigned int Nmeas,
                                dogleg_callback_dense_t* f, void* cookie);

This function computes the jacobian with center differences and compares the result with the jacobian computed by the callback function. It is recommended to do this for every variable while developing the program that uses libdogleg.

• "var" is the index of the variable being tested
• "p0" is the state vector p where we're evaluating the jacobian
• "Nstate", "NJnnz" are the number of state variables, measurements
• "f" is the callback function, as before
• "cookie" is the user data, as before

This function returns nothing, but prints out the test results.

dogleg_callback_t

The main user callback that specifies the sparse optimization problem has type

 typedef void (dogleg_callback_t)(const double*   p,
                                  double*         x,
                                  cholmod_sparse* Jt,
                                  void*           cookie);

• p is the current state vector
• x is the resulting f(p)
• Jt is the transpose of df/dp at p. As mentioned previously, Jt is stored column-first by CHOLMOD, which can be interpreted as storing J row-first by the user-defined callback routine
• The "cookie" is the user-defined arbitrary data passed into "dogleg_optimize".

dogleg_callback_dense_t

The main user callback that specifies the dense optimization problem has type

 typedef void (dogleg_callback_dense_t)(const double*   p,
                                        double*         x,
                                        double*         J,
                                        void*           cookie);

• p is the current state vector
• x is the resulting f(p)
• J is df/dp at p. J is stored row-first, with all the derivatives for the first measurement, then all the derivatives for the second measurement and so on.
• The "cookie" is the user-defined arbitrary data passed into "dogleg_optimize".

dogleg_solverContext_t

This is the solver context that can be retrieved through the "returnContext" parameter of the "dogleg_optimize" call. This structure contains all the internal state used by the solver. If requested, the user is responsible for calling "dogleg_freeContext" when done. This structure is defined as:

 typedef struct
 {
   cholmod_common  common;
   union
   {
     dogleg_callback_t*       f;
     dogleg_callback_dense_t* f_dense;
   };
   void*              cookie;
   // between steps, beforeStep contains the operating point of the last step.
   // afterStep is ONLY used while making the step. Externally, use beforeStep
   // unless you really know what you're doing
   dogleg_operatingPoint_t* beforeStep;
   dogleg_operatingPoint_t* afterStep;
   // The result of the last JtJ factorization performed. Note that JtJ is not
   // necessarily factorized at every step, so this is NOT guaranteed to contain
   // the factorization of the most recent JtJ
   union
   {
     cholmod_factor* factorization;
     // This is a factorization of JtJ, stored as a packed symmetric matrix
     // returned by dpptrf('L', ...)
     double*         factorization_dense;
   };
   // Have I ever seen a singular JtJ? If so, I add this constant to the diagonal
   // from that point on. This is a simple and fast way to deal with
   // singularities. This constant starts at 0, and is increased every time a
   // singular JtJ is encountered. This is suboptimal but works for me for now.
   double                   lambda;
   // Are we using sparse math (cholmod)?
   int                      is_sparse;
   int                      Nstate, Nmeasurements;
 } dogleg_solverContext_t;

Some of the members are copies of the data passed into "dogleg_optimize"; some others are internal state. Of potential interest are

• "common" is a cholmod_common structure used by all CHOLMOD calls. This can be used for any extra CHOLMOD work the user may want to do
• "beforeStep" contains the operating point of the optimum solution. The user can analyze this data without the need to re-call the callback routine.

dogleg_operatingPoint_t

An operating point of the solver. This is a part of "dogleg_solverContext_t". Various variables describing the operating point such as p, J, x, norm2(x) and Jt x are available. All of the just-mentioned variables are computed at every step and are thus always valid.

 // an operating point of the solver
 typedef struct
 {
   double*         p;
   double*         x;
   double          norm2_x;
   union
   {
     cholmod_sparse* Jt;
     double*         J_dense; // row-first: grad0, grad1, grad2, ...
   };
   double*         Jt_x;
   // the cached update vectors. It's useful to cache these so that when a step is rejected, we can
   // reuse these when we retry
   double*        updateCauchy;
   union
   {
     cholmod_dense* updateGN_cholmoddense;
     double*        updateGN_dense;
   };
   double         updateCauchy_lensq, updateGN_lensq; // update vector lengths
   // whether the current update vectors are correct or not
   int updateCauchy_valid, updateGN_valid;
   int didStepToEdgeOfTrustRegion;
 } dogleg_operatingPoint_t;

Parameters¶

The optimization is controlled by several parameters. These can be set globally for all callers of "libdogleg" in a process using the "dogleg_set....()" functions. Those global values are used by "dogleg_optimize" and "dogleg_optimize_dense". Or these can be specified independently for each invocation by passing a "parameters" argument to "dogleg_optimize2" or "dogleg_optimize_dense2". The latter is recommended because multiple instances of libdogleg in a single application would no longer conflict.

It is not required to set any of these, but it's highly recommended to set the initial trust-region size and the termination thresholds to match the problem being solved. Furthermore, it's highly recommended for the problem being solved to be scaled so that every state variable affects the objective norm2( x ) equally. This makes this method's concept of "trust region" much more well-defined and makes the termination criteria work correctly.

dogleg_setMaxIterations

To set the maximum number of solver iterations, call

 void dogleg_setMaxIterations(int n);

dogleg_setDebug

To turn on diagnostic output, call

 void dogleg_setDebug(int debug);

with a non-zero value for "debug". Two separate diagnostic streams are available: a verbose human-oriented stream, and a vnlog <http://github.com/dkogan/vnlog>.

By default, diagnostic output is disabled. The "debug" argument is a bit-mapped integer:

 if(debug == 0                 ): no diagnostic output
 if(debug &  DOGLEG_DEBUG_VNLOG): output vnlog diagnostics to stdout
 if(debug & ~DOGLEG_DEBUG_VNLOG): output human-oriented diagnostics to stderr

"DOGLEG_DEBUG_VNLOG" has a very high value, so if human diagnostics are desired, the recommended call is:

 dogleg_setDebug(1);

dogleg_setInitialTrustregion

The optimization method keeps track of a trust region size. Here, the trust region is a ball in R^Nstate. When the method takes a step p -> p + delta_p, it makes sure that

sqrt( norm2( delta_p ) ) < trust region size.

The initial value of the trust region size can be set with

 void dogleg_setInitialTrustregion(double t);

The dogleg algorithm is efficient when recomputing a rejected step for a smaller trust region, so set the initial trust region size to a value larger to a reasonable estimate; the method will quickly shrink the trust region to the correct size.

dogleg_setThresholds

The routine exits when the maximum number of iterations is exceeded, or a termination threshold is hit, whichever happens first. The termination thresholds are all designed to trigger when very slow progress is being made. If all went well, this slow progress is due to us finding the optimum. There are 3 termination thresholds:

• The function being minimized is E = norm2( x ) where x = f(p).

  dE/dp = 2*Jt*x where Jt is transpose(dx/dp).

   if( for every i  fabs(Jt_x[i]) < JT_X_THRESHOLD )
   { we are done; }
      

• The method takes discrete steps: p -> p + delta_p

   if( for every i  fabs(delta_p[i]) < UPDATE_THRESHOLD)
   { we are done; }
      

• The method dynamically controls the trust region.

   if(trustregion < TRUSTREGION_THRESHOLD)
   { we are done; }
      

To set these threholds, call

 void dogleg_setThresholds(double Jt_x, double update, double trustregion);

To leave a particular threshold alone, specify a negative value.

dogleg_setTrustregionUpdateParameters

This function sets the parameters that control when and how the trust region is updated. The default values should work well in most cases, and shouldn't need to be tweaked.

Declaration looks like

 void dogleg_setTrustregionUpdateParameters(double downFactor, double downThreshold,
                                            double upFactor,   double upThreshold);

To see what the parameters do, look at "evaluateStep_adjustTrustRegion" in the source. Again, these should just work as is.