GSL.cc

Go to the documentation of this file.

/*
 * $Id: GSL.cc 135 2013-06-29 15:09:42Z jdl3 $
 * Copyright (C) 2013 John D Lamb
 *
 * 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 2 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#ifdef HAVE_CONFIG_H
#  include<config.h>
#endif

#include"GSL.hpp"

using namespace gsl;

std::tuple<gsl::matrix,double>
ipo::detail::modifiedCholeskyDecomp( gsl::matrix const& matrix, double maxoffl ){
  // constant expressions
  static double constexpr
    eps4 { std::sqrt( std::sqrt( std::numeric_limits<double>::epsilon() ) ) };
  static double constexpr eps2 { std::sqrt( std::numeric_limits<double>::epsilon() ) };
  if( 0 == maxoffl ){
    auto D = matrix.const_diagonal();
    maxoffl = std::sqrt( std::accumulate( D.begin(), D.end(),
                      -std::numeric_limits<double>::infinity(),
                      []( double const x, double const y ){
                        return x > y ? x : y; } ) );
  }
  double const minL2 { eps2 * maxoffl };
  
  // Keep local storage: we are likely to reuse it frequently.
  static std::vector<double> D;

  size_t const N = matrix.size1();
  if( matrix.size2() != N ) throw exception( "gsl::modifiedCholeskyDecomp():"
                         "matrix must be square.", __FILE__,
                         __LINE__, exception::GSL_EBADLEN );
  // assumes maxoffl != 0
  if( maxoffl == 0 )
    throw exception( "gsl::modifiedCholeskyDecomp():"
             "maxoffl must be nonzero.", __FILE__,
             __LINE__, exception::GSL_ESANITY );

  gsl::matrix L( N, N );  // return matrix
  D.resize( N ); // This should not usually adjust the capacity of D

  double const minL { eps4 * maxoffl };
  
  double maxAdd { 0 };
  for( size_t j { 0 }; j < N; ++j ){

    double sumL2 { 0 };
    for( size_t i { 0 }; i < j; ++i ){
      double const entry { L.get( j, i ) };
      sumL2 += entry * entry;
    }
    L.set( j, j, matrix.get( j, j ) - sumL2 ); // Only diagonal set
    
    double minLjj { 0 };
    for( size_t i { j + 1 }; i < N; ++i ){
      double sumLL { 0 };
      for( size_t k { 0 }; k < j; ++k )
    sumLL += L.get( i, k ) * L.get( j, k );
      double const d { matrix.get( j, i ) - sumLL };
      L.set( i, j, d );
      minLjj = std::max( minLjj, d );
    }
    minLjj = std::max( minLjj / maxoffl, minL );

    if( L.get( j, j ) > minLjj * minLjj ){
      // Normal Cholesky iteration
      L.set( j, j, std::sqrt( L.get( j, j ) ) );
    } else {
      // Augment matrix
      minLjj = std::max( minLjj, minL2 );
      maxAdd = std::max( maxAdd, minLjj * minLjj - L.get( j, j ) );
      L.set( j, j, minLjj );
    }

    // complete decomposition
    for( size_t i { j + 1 }; i < N; ++i )
      L.set( i, j, L.get( i, j ) / L.get( j, j ) );
  }
  
  // return matrix and maxAdd
  return std::make_tuple( L, maxAdd );
}

gsl::matrix
ipo::detail::modelHessian( gsl::matrix& matrix ){
  // constant expressions
  static double constexpr eps2 { std::sqrt( std::numeric_limits<double>::epsilon() )};
  static double constexpr INF { std::numeric_limits<double>::infinity() };

  size_t const N { matrix.size1() };
  if( matrix.size2() != N ) throw exception( "gsl::modelHessian():"
                         "matrix must be square.", __FILE__,
                         __LINE__, exception::GSL_EBADLEN );
  
  double maxDiag { -INF };
  double minDiag { INF };
  for( double const& d : matrix.const_diagonal() ){
    maxDiag = std::max( maxDiag, d );
    minDiag = std::min( minDiag, d );
  }
  double const maxPosDiag { std::max( 0.0, maxDiag ) };
  
  double mu { 0 };
  if( minDiag <= eps2 * maxPosDiag ){
    mu = 2 * (maxPosDiag - minDiag) * eps2 - minDiag;
    maxDiag = mu;
  }

  double maxoff { 0 };
  for( size_t i { 0 }; i < N; ++i )
    for( size_t j { i + 1 }; j < N; ++j )
      maxoff = std::max( maxoff, matrix.get( i, j ) );
  if( maxoff * (1 + 2 * eps2) > maxDiag ){
    mu += (maxoff - maxDiag) + 2 * eps2 * maxoff;
    maxDiag = maxoff * (1 + 2 * eps2);
  }
  
  if( maxDiag == 0 ){ /* H == 0 */
    mu = 1;
    maxDiag = 1;
  }
  
  if( mu > 0 ) /* Add mu * I to H */ matrix.diagonal().add_constant( mu );
  
  double const maxoffl { std::sqrt( std::max( maxDiag, maxoff / N ) ) };
  
  auto cholDecomp = modifiedCholeskyDecomp( matrix, maxoffl );
  gsl::matrix L { std::get<0>( cholDecomp ) };
  double const maxAdd { std::get<1>( cholDecomp ) };

  if( maxAdd > 0 ){ /* H not positive definite */
    double maxev { matrix.get( 0, 0 ) };
    double minev { maxev };
    for( size_t i { 0 }; i < N; ++i ){
      double offrow { 0 };
      for( size_t j { 0 }; j + 1 < i; ++j )
    offrow += std::fabs( matrix.get( j, i ) );
      for( size_t j { i + 1 }; j < N; ++j )
    offrow += std::fabs( matrix.get( i, j ) );
      double const hii { matrix.get( i, i ) };
      maxev = std::max( maxev, hii + offrow );
      minev = std::min( minev, hii - offrow );
    }
    double const sdd { std::max( 0.0, (maxev - minev) *eps2 - minev ) };
    mu = std::min( maxAdd, sdd );
    matrix.diagonal().add_constant( mu ); /* H = H + mu * I */
  }
  cholDecomp = modifiedCholeskyDecomp( matrix, 0 );
  L = std::get<0>( cholDecomp );
  /* Fill upper triangle */
  for( size_t j { 0 }; j < N; ++j )
    for( size_t i { 0 }; i < j; ++i )
      L.set( i, j, L.get( j, i ) );

  return L;
}