dampedNewtonForNeymanOSDE1D.hh

Go to the documentation of this file.

#ifndef NPSTAT_DAMPEDNEWTONFORNEYMANOSDE1D_HH_
#define NPSTAT_DAMPEDNEWTONFORNEYMANOSDE1D_HH_
 
/*!
// \file dampedNewtonForNeymanOSDE1D.hh
//
// \brief Damped Newton's method for the Neyman OSDE
//
// Author: I. Volobouev
//
// March 2023
*/
 
#include <cmath>
#include <cassert>
#include <utility>
#include <algorithm>
 
#include "npstat/nm/maxAbsValue.hh"
 
#include "npstat/stat/NeymanOSDE1D.hh"
#include "npstat/stat/NeymanOSDE1DResult.hh"
#include "npstat/stat/NeymanOSDE1DConvergence.hh"
#include "npstat/stat/gradientDescentForNeymanOSDE1D.hh"
 
namespace npstat {
    // The class "ConvergenceCalculator" must have a method
    // "bool converged(...) const" with the same signature as the
    // "converged" method of class "GradientNeymanOSDE1DConvergenceCalc".
    template<class ConvergenceCalculator>
    inline NeymanOSDE1DResult dampedNewtonForNeymanOSDE1D(
        const NeymanOSDE1D& nosde,
        const double* initialCoeffs, const unsigned nCoeffs,
        const double* shrinkages, const unsigned nShrinkages,
        const ConvergenceCalculator& conv, const unsigned maxIterations,
        const double dampingCoefficient = 1.0,
        const bool useGradientDescent = false)
    {
        assert(dampingCoefficient > 0.0);
 
        std::vector<double> coeffs(initialCoeffs, initialCoeffs+nCoeffs);
        std::vector<double> polyStats(nShrinkages);
        std::vector<double> prevGradient(nCoeffs);
        std::vector<double> currentGradient(nCoeffs);
        std::vector<double> hessianEigenvalues(nCoeffs);
        std::vector<double> step(nCoeffs);
        Matrix<double> hessian(nCoeffs, nCoeffs);
        Matrix<double> hessianEigenvectors(nCoeffs, nCoeffs);
 
        double prevChisq = nosde.chiSquare(&coeffs[0], nCoeffs,
                                           shrinkages, nShrinkages,
                                           &polyStats[0], &prevGradient[0],
                                           &hessian);
        bool converged = false, badHessian = false;
        unsigned iter = 0;
        for (; iter<maxIterations && !converged; ++iter)
        {
            hessian.symEigen(&hessianEigenvalues[0], nCoeffs, &hessianEigenvectors);
            for (unsigned i=0; i<nCoeffs && !badHessian; ++i)
            {
                if (hessianEigenvalues[i] <= 0.0)
                    badHessian = true;
                else
                    hessianEigenvalues[i] = 1.0/hessianEigenvalues[i];
            }
 
            if (badHessian)
            {
                // Try to switch to gradient descent.
                // Note that gradient descent by itself is very slow.
                // We are going to terminate it as soon as it brings
                // the search into a region with a positive-definite
                // Hessian.
                if (useGradientDescent)
                {
                    HessianNeymanOSDE1DConvergenceCalc<ConvergenceCalculator> descentConv(conv);
                    const NeymanOSDE1DResult& r = gradientDescentForNeymanOSDE1D(
                        nosde, &coeffs[0], nCoeffs, shrinkages, nShrinkages,
                        descentConv, maxIterations);
                    if (r.converged())
                    {
                        if (descentConv.otherCalcConverged())
                            return r;
                        else
                        {
                            // Recalculate the Hessian
                            for (unsigned i=0; i<nCoeffs; ++i)
                                coeffs[i] = r.coeffs()[i];
                            prevChisq = nosde.chiSquare(&coeffs[0], nCoeffs,
                                                        shrinkages, nShrinkages,
                                                        &polyStats[0], &prevGradient[0],
                                                        &hessian);
                            badHessian = false;
                        }
                    }
                }
 
                // Did we fix the Hessian? If yes, the new Hessian
                // has to be inverted, so we have to go to the
                // beginning of the cycle. If not, we can't do much
                // with this optimization method.
                if (badHessian)
                    break;
                else
                    continue;
            }
 
            const Matrix<double>& invEigen = diag(&hessianEigenvalues[0], nCoeffs);
            const Matrix<double>& invHess = invEigen.bilinear(hessianEigenvectors);
            invHess.timesVector(&prevGradient[0], nCoeffs, &step[0], nCoeffs);
            for (unsigned i=0; i<nCoeffs; ++i)
            {
                step[i] *= -dampingCoefficient;
                coeffs[i] += step[i];
            }
 
            const double currentChisq = nosde.chiSquare(&coeffs[0], nCoeffs,
                                                        shrinkages, nShrinkages,
                                                        &polyStats[0], &currentGradient[0],
                                                        &hessian);
            converged = conv.converged(iter, prevChisq, prevGradient, step,
                                       currentChisq, currentGradient, hessian);
 
            prevChisq = currentChisq;
            std::swap(prevGradient, currentGradient);
        }
 
        return NeymanOSDE1DResult(coeffs, polyStats, prevChisq,
                                  maxAbsValue(&prevGradient[0], prevGradient.size()),
                                  NOSDE_NEWTON, iter, converged, badHessian);
    }
}
 
#endif // NPSTAT_DAMPEDNEWTONFORNEYMANOSDE1D_HH_