|
Packit |
67cb25 |
/* multifit_nlinear/fdfvv.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2015 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 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_multifit_nlinear.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_matrix.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
fdfvv()
|
|
Packit |
67cb25 |
Compute approximate second directional derivative using
|
|
Packit |
67cb25 |
finite differences.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
See Eq. 19 of:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
M. K. Transtrum, J. P. Sethna, Improvements to the Levenberg
|
|
Packit |
67cb25 |
Marquardt algorithm for nonlinear least-squares minimization,
|
|
Packit |
67cb25 |
arXiv:1201.5885, 2012.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: h - step size for finite difference
|
|
Packit |
67cb25 |
x - parameter vector, size p
|
|
Packit |
67cb25 |
v - geodesic velocity, size p
|
|
Packit |
67cb25 |
f - vector of function values f_i(x), size n
|
|
Packit |
67cb25 |
J - Jacobian matrix J(x), n-by-p
|
|
Packit |
67cb25 |
swts - data weights
|
|
Packit |
67cb25 |
fdf - fdf struct
|
|
Packit |
67cb25 |
fvv - (output) approximate second directional derivative
|
|
Packit |
67cb25 |
vector D_v^2 f(x)
|
|
Packit |
67cb25 |
work - workspace, size p
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success or error
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
fdfvv(const double h, const gsl_vector *x, const gsl_vector *v,
|
|
Packit |
67cb25 |
const gsl_vector *f, const gsl_matrix *J, const gsl_vector *swts,
|
|
Packit |
67cb25 |
gsl_multifit_nlinear_fdf *fdf, gsl_vector *fvv, gsl_vector *work)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
const size_t n = fdf->n;
|
|
Packit |
67cb25 |
const size_t p = fdf->p;
|
|
Packit |
67cb25 |
const double hinv = 1.0 / h;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute work = x + h*v */
|
|
Packit |
67cb25 |
for (i = 0; i < p; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = gsl_vector_get(x, i);
|
|
Packit |
67cb25 |
double vi = gsl_vector_get(v, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set(work, i, xi + h * vi);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute f(x + h*v) */
|
|
Packit |
67cb25 |
status = gsl_multifit_nlinear_eval_f (fdf, work, swts, fvv);
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < n; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double fi = gsl_vector_get(f, i); /* f_i(x) */
|
|
Packit |
67cb25 |
double fip = gsl_vector_get(fvv, i); /* f_i(x + h*v) */
|
|
Packit |
67cb25 |
gsl_vector_const_view row = gsl_matrix_const_row(J, i);
|
|
Packit |
67cb25 |
double u, fvvi;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute u = sum_{ij} J_{ij} D v_j */
|
|
Packit |
67cb25 |
gsl_blas_ddot(&row.vector, v, &u);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fvvi = (2.0 * hinv) * ((fip - fi) * hinv - u);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set(fvv, i, fvvi);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_multifit_nlinear_fdfvv()
|
|
Packit |
67cb25 |
Compute approximate second directional derivative
|
|
Packit |
67cb25 |
using finite differences
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: h - step size for finite difference
|
|
Packit |
67cb25 |
x - parameter vector, size p
|
|
Packit |
67cb25 |
v - geodesic velocity, size p
|
|
Packit |
67cb25 |
f - function values f_i(x), size n
|
|
Packit |
67cb25 |
J - Jacobian matrix J(x), n-by-p
|
|
Packit |
67cb25 |
swts - sqrt data weights (set to NULL if not needed)
|
|
Packit |
67cb25 |
fdf - fdf
|
|
Packit |
67cb25 |
fvv - (output) approximate (weighted) second directional derivative
|
|
Packit |
67cb25 |
vector, size n, sqrt(W) fvv
|
|
Packit |
67cb25 |
work - workspace, size p
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success or error
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_multifit_nlinear_fdfvv(const double h, const gsl_vector *x, const gsl_vector *v,
|
|
Packit |
67cb25 |
const gsl_vector *f, const gsl_matrix *J,
|
|
Packit |
67cb25 |
const gsl_vector *swts, gsl_multifit_nlinear_fdf *fdf,
|
|
Packit |
67cb25 |
gsl_vector *fvv, gsl_vector *work)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return fdfvv(h, x, v, f, J, swts, fdf, fvv, work);
|
|
Packit |
67cb25 |
}
|