pysolvegn.solve#

solve(terms, x0, *, max_iteration=None, max_time=None, ftol=None, xtol=None, gtol=None, atol=None, ptol=None, callback_func=None, update_func=None, verbosity=0, naninf=True, history=False, history_details=None)[source]#

Function to solve the least squares problem using the Gauss-Newton method with robust cost functions.

The function accepts multiple terms in the least squares problem, each with its own residual function, Jacobian function, weight, and loss function solving the following optimization problem :

\[\min_{\mathbf{p}} \frac{1}{2} \sum_{i} w_i \sum_j \rho_i\left(\| \mathbf{r}_{i,j}(\mathbf{p}) \|^2\right)\]

where \(w_i\) is a weight for each sub least squares problem, and \(\rho_i\) is a robust cost function for each sub least squares problem.

See also

For more details on notations for the optimization problem and the algorithm, please refer to the mathematical section of the documentation.

The cost of each term defined by the residuals is estimated as :

\[\frac{1}{2} \sum_j \rho_i\left(\| \mathbf{R}_{i,j}(\mathbf{p}) \|^2\right)\]

Important

For terms defined directly by a Hessian and gradient, only a local quadratic model of the cost is available. Since this model is defined up to an arbitrary constant, an absolute cost value cannot be uniquely determined. For reporting consistency, the cost contribution of these terms is set to 0.

Parameters:

x0 (ArrayLike) – The initial guess for the parameters of the optimization with shape (n_parameters,). The array will not be modified by this function, a copy of the parameters will be used for the optimization process.
max_iteration (Optional[Integral], optional (default=None)) – Stop criterion by the number of iterations. The optimization process is stopped when the number of iterations exceeds max_iteration. If None, no limit on the number of iterations is considered.
max_time (Optional[Real], optional (default=None)) – Stop criterion by the elapsed time of optimization. The optimization process is stopped when the time elapsed since the beginning of the optimization exceeds max_time seconds. If None, no limit on the computation time is considered.
ftol (Optional[Real], optional (default=None)) – Stop criterion by the change of the cost function value. The optimization process is stopped when dF < ftol * F where F is the cost function value and dF is the change of the cost function value between two iterations. If None, this criterion is not considered.
xtol (Optional[Real], optional (default=None)) – Stop criterion by the change of the parameters. The optimization process is stopped when ||dp|| < xtol * (xtol + ||p||) where p is the parameters and dp is the change of the parameters between two iterations. If None, this criterion is not considered.
gtol (Optional[Real], optional (default=None)) – Stop criterion by the optimality value. The optimization process is stopped when the optimality verifies norm(g, ord=numpy.inf) < gtol where g is the scaled second term. If None, this criterion is not considered.
atol (Optional[Real], optional (default=None)) – Stop criterion by the absolute cost function value. The optimization process is stopped when F < atol where F is the cost function value. If None, this criterion is not considered.
ptol (Optional[Real], optional (default=None)) – Stop criterion by the change of the parameters value. The optimization process is stopped when norm(dp, ord=numpy.inf) < ptol where dp is the change of the parameters. If None, this criterion is not considered.
callback_func (Optional[Callable[[Dict], bool]], optional (default=None)) – A callback function that is called at the end of each iteration of the optimization process. The callback function should take a dictionary similar to the one returned in the history of the optimization process as input. If the return value of the callback function is True (or None), the optimization process will continue. If the return value is False, the optimization process will be stopped. If None, no callback function is used.
update_func (Optional[Callable[[numpy.ndarray, numpy.ndarray], ArrayLike]], optional (default=None)) – A function that is called at the end of each iteration of the optimization process to compute the parameters for the next iteration. The function should take the current parameters and update (p, dp) as inputs and return the parameters with shape (n_parameters,) to use for the next iteration. If None, the parameters for the next iteration will be computed as p + dp where p is the current parameters and dp is the update computed by solving the linear system.
verbosity (Integral, optional (default=0)) – The level of verbosity for logging the optimization process. 0: No logging 1: Log only the final results of the optimization process. 2: Log the results at each iteration of the optimization process. 3: Details logging for debugging purposes.
naninf (bool, optional (default=True)) – If True, the optimization process will be stopped if any NaN or Inf value is encountered in the cost function value or parameters during the optimization process. If False, the optimization process will continue.
history (bool (optional, default=False)) – If True, the function will also return a tuple containing the history of the optimization process.
history_details (Optional[Union[str, Sequence[str]]], optional (default=None)) – The details to include in the history of the optimization process. Can be a single string or a sequence of strings. See the Notes section for the available details. By default, the history will include : “iteration”, “elapsed_time”, “parameters”, “delta_params”, “delta_cost”, “cost”, and “optimality” details. For a specified list of details, only the details in the list below will be included in the history.
terms (Sequence[Term])

Returns:

parameters (numpy.ndarray) – The parameters that minimize the least squares problem with shape (n_parameters,).
history (List[Dict], optional) – Only returned if history is True. A list containing the history of the optimization process. Each element of the list is a dictionary containing the keys described below.

Return type:

Tuple[ndarray, List[Dict] | None]

Notes

The history contains the following keys (if requested in history_details):

“iteration”: Integer representing the iteration number.
“elapsed_time”: Float representing the time elapsed since the beginning of the optimization process in seconds.
“parameters”: Numpy array representing the parameters at the current iteration.
“delta_params”: Numpy array representing the change of parameters between the current and previous iteration (Only for iterations > 0).
“costs”: List of floats representing the cost function value for each term in the least squares problem at the current iteration compute as \(\frac{1}{2} \sum_j \rho_i\left(\| \mathbf{R}_{i,j}(\mathbf{p}) \|^2\right)\) for each term \(i\).
“cost”: Float representing the cost function value at the current iteration computed as \(\frac{1}{2} \sum_i w_i \sum_j \rho_i\left(\| \mathbf{R}_{i,j}(\mathbf{p}) \|^2\right)\).
“delta_cost”: Float representing the change of the cost function value between the current and previous iteration (Only for iterations > 0).
“optimality”: Float representing the optimality value at the current iteration computed as norm(b, ord=numpy.inf) where b is the scaled second term.
“residuals”: A list of numpy arrays representing the residuals for each term in the least squares problem at the current iteration (Only for rJ terms).
“second_term”: The scaled second term \(\mathbf{g}\) of the linear system at the current iteration.
“hessian”: The Hessian approximation \(\mathbf{H}\) of the linear system at the current iteration.

The optimization process is performed by iteratively solving a linearized version of the least squares problem at each iteration for different values of the regularization factor, and then selecting the optimal regularization factor based on the L-curve analysis. The step of each loop iteration is as follows:

While NOT converged:
Build the system ``H Δp = -g``.
Call the 'callback' function and STOP if it returns False.
Check the stopping criteria and STOP if any of them is satisfied.
If CONTINUE solve the linear system ``H Δp = -g`` for Δp.
Update the parameters as ``p = p + Δp`` (or using 'update_func' if provided).

Version#

0.1.0: Initial version for the solver method.