scaling.c

/* multilarge_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_multilarge_nlinear.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>

static int init_diag_levenberg(const gsl_matrix * JTJ, gsl_vector * diag);
static int update_diag_levenberg(const gsl_matrix * JTJ, gsl_vector * diag);

static int init_diag_marquardt(const gsl_matrix * JTJ, gsl_vector * diag);
static int update_diag_marquardt (const gsl_matrix * JTJ, gsl_vector * diag);

static int init_diag_more(const gsl_matrix * JTJ, gsl_vector * diag);
static int update_diag_more(const gsl_matrix * JTJ, gsl_vector * diag);

/* Levenberg scaling, D = I */
static int
init_diag_levenberg(const gsl_matrix * JTJ, gsl_vector * diag)
{
  (void)JTJ; /* avoid unused parameter warning */
  gsl_vector_set_all(diag, 1.0);
  return GSL_SUCCESS;
}

static int
update_diag_levenberg(const gsl_matrix * JTJ, gsl_vector * diag)
{
  (void)JTJ;  /* 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 * JTJ, gsl_vector * diag)
{
  return update_diag_marquardt(JTJ, diag);
}

/* update diagonal scaling matrix D according to Marquardt method */
static int
update_diag_marquardt (const gsl_matrix * JTJ, gsl_vector * diag)
{
  const size_t p = JTJ->size2;
  size_t j;

  for (j = 0; j < p; j++)
    {
      double Jjj = gsl_matrix_get(JTJ, j, j);
      double norm;

      if (Jjj <= 0.0)
        norm = 1.0;
      else
        norm = sqrt(Jjj);

      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 * JTJ, gsl_vector * diag)
{
  int status;

  gsl_vector_set_zero(diag);
  status = update_diag_more(JTJ, diag);

  return status;
}

/* update diagonal scaling matrix D according to Eq. 6.3 of
 * More, 1978 */
static int
update_diag_more (const gsl_matrix * JTJ, gsl_vector * diag)
{
  const size_t p = JTJ->size2;
  size_t j;

  for (j = 0; j < p; j++)
    {
      double Jjj = gsl_matrix_get(JTJ, j, j);
      double *diagj = gsl_vector_ptr(diag, j);
      double norm;

      if (Jjj <= 0.0)
        norm = 1.0;
      else
        norm = sqrt(Jjj);

      *diagj = GSL_MAX(*diagj, norm);
    }

  return GSL_SUCCESS;
}

static const gsl_multilarge_nlinear_scale levenberg_type =
{
  "levenberg",
  init_diag_levenberg,
  update_diag_levenberg
};

static const gsl_multilarge_nlinear_scale marquardt_type =
{
  "marquardt",
  init_diag_marquardt,
  update_diag_marquardt
};

static const gsl_multilarge_nlinear_scale more_type =
{
  "more",
  init_diag_more,
  update_diag_more
};

const gsl_multilarge_nlinear_scale *gsl_multilarge_nlinear_scale_levenberg = &levenberg_type;
const gsl_multilarge_nlinear_scale *gsl_multilarge_nlinear_scale_marquardt = &marquardt_type;
const gsl_multilarge_nlinear_scale *gsl_multilarge_nlinear_scale_more = &more_type;