/* multifit_nlinear/scaling.c
*
* Copyright (C) 2015, 2016 Patrick Alken
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
/*
* This module handles the updating of the scaling matrix D_k in the
* trust region subproblem:
*
* min m_k (dx), || D_k dx || <= Delta_k
*
* where m_k(dx) is a model which approximates the cost function
* F(x_k + dx) near the current iteration point x_k
*
* D_k can be updated according to several different strategies.
*/
#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
static int init_diag_levenberg(const gsl_matrix * J, gsl_vector * diag);
static int update_diag_levenberg(const gsl_matrix * J,
gsl_vector * diag);
static int init_diag_marquardt(const gsl_matrix * J, gsl_vector * diag);
static int update_diag_marquardt (const gsl_matrix * J,
gsl_vector * diag);
static int init_diag_more(const gsl_matrix * J, gsl_vector * diag);
static int update_diag_more(const gsl_matrix * J, gsl_vector * diag);
/* Levenberg scaling, D = I */
static int
init_diag_levenberg(const gsl_matrix * J, gsl_vector * diag)
{
(void)J; /* avoid unused parameter warning */
gsl_vector_set_all(diag, 1.0);
return GSL_SUCCESS;
}
static int
update_diag_levenberg(const gsl_matrix * J, gsl_vector * diag)
{
(void)J; /* avoid unused parameter warning */
(void)diag; /* avoid unused parameter warning */
/* nothing to do */
return GSL_SUCCESS;
}
/* initialize diagonal scaling matrix D according to Marquardt method */
static int
init_diag_marquardt(const gsl_matrix * J, gsl_vector * diag)
{
return update_diag_marquardt(J, diag);
}
/* update diagonal scaling matrix D according to Marquardt method */
static int
update_diag_marquardt (const gsl_matrix * J, gsl_vector * diag)
{
const size_t p = J->size2;
size_t j;
for (j = 0; j < p; j++)
{
gsl_vector_const_view v = gsl_matrix_const_column(J, j);
double norm = gsl_blas_dnrm2(&v.vector);
if (norm == 0.0)
norm = 1.0;
gsl_vector_set(diag, j, norm);
}
return GSL_SUCCESS;
}
/* initialize diagonal scaling matrix D according to Eq 6.3 of
* More, 1978 */
static int
init_diag_more(const gsl_matrix * J, gsl_vector * diag)
{
int status;
gsl_vector_set_zero(diag);
status = update_diag_more(J, diag);
return status;
}
/* update diagonal scaling matrix D according to Eq. 6.3 of
* More, 1978 */
static int
update_diag_more (const gsl_matrix * J, gsl_vector * diag)
{
const size_t p = J->size2;
size_t j;
for (j = 0; j < p; j++)
{
gsl_vector_const_view v = gsl_matrix_const_column(J, j);
double norm = gsl_blas_dnrm2(&v.vector);
double *diagj = gsl_vector_ptr(diag, j);
if (norm == 0.0)
norm = 1.0;
*diagj = GSL_MAX(*diagj, norm);
}
return GSL_SUCCESS;
}
static const gsl_multifit_nlinear_scale levenberg_type =
{
"levenberg",
init_diag_levenberg,
update_diag_levenberg
};
static const gsl_multifit_nlinear_scale marquardt_type =
{
"marquardt",
init_diag_marquardt,
update_diag_marquardt
};
static const gsl_multifit_nlinear_scale more_type =
{
"more",
init_diag_more,
update_diag_more
};
const gsl_multifit_nlinear_scale *gsl_multifit_nlinear_scale_levenberg = &levenberg_type;
const gsl_multifit_nlinear_scale *gsl_multifit_nlinear_scale_marquardt = &marquardt_type;
const gsl_multifit_nlinear_scale *gsl_multifit_nlinear_scale_more = &more_type;