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