gradientDescentForNeymanOSDE1D.hh

Go to the documentation of this file.

#ifndef NPSTAT_GRADIENTDESCENTFORNEYMANOSDE1D_HH_
#define NPSTAT_GRADIENTDESCENTFORNEYMANOSDE1D_HH_
 
/*!
// \file gradientDescentForNeymanOSDE1D.hh
//
// \brief Gradient descent method for the Neyman OSDE. Not very well
//        thought out, so it is very slow.
//
// Author: I. Volobouev
//
// March 2023
*/
 
#include <cmath>
#include <cassert>
#include <utility>
#include <algorithm>
 
#include "npstat/nm/maxAbsValue.hh"
#include "npstat/nm/scalesFromHessian.hh"
#include "npstat/nm/SimpleScalarProduct.hh"
 
#include "npstat/stat/NeymanOSDE1D.hh"
#include "npstat/stat/NeymanOSDE1DResult.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 gradientDescentForNeymanOSDE1D(
        const NeymanOSDE1D& nosde,
        const double* initialCoeffs, const unsigned nCoeffs,
        const double* shrinkages, const unsigned nShrinkages,
        const ConvergenceCalculator& conv, const unsigned maxIterations)
    {
        // Some hardwired parameters
        const double armijoC = 0.5;
        const double armijoTau = 0.5;
        const unsigned maxLineIterations = 100;
 
        std::vector<double> coeffs(initialCoeffs, initialCoeffs+nCoeffs);
        std::vector<double> lineSearchCoeffs(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);
        std::vector<double> unitStep(nCoeffs);
        std::vector<double> scales(nCoeffs);
        Matrix<double> hessian(nCoeffs, nCoeffs);
        SimpleScalarProduct<double> sp;
 
        double prevChisq = nosde.chiSquare(&coeffs[0], nCoeffs,
                                           shrinkages, nShrinkages,
                                           &polyStats[0], &prevGradient[0],
                                           &hessian);
        bool converged = false;
        unsigned iter = 0;
        for (; iter<maxIterations && !converged; ++iter)
        {
            // Try to get some idea about possible step sizes
            scalesFromHessian(hessian, 1.0, &scales[0], nCoeffs);
 
            // Figure out the "unit step" in the direction
            // opposing the gradient
            const double gradientNorm = std::sqrt(sp(&prevGradient[0], &prevGradient[0], nCoeffs));
            assert(gradientNorm > 0.0);
            for (unsigned i=0; i<nCoeffs; ++i)
                unitStep[i] = -prevGradient[i]*scales[i]/gradientNorm;
            const double unitStepLen = std::sqrt(sp(&unitStep[0], &unitStep[0], nCoeffs));
            assert(unitStepLen > 0.0);
            for (unsigned i=0; i<nCoeffs; ++i)
                unitStep[i] = -prevGradient[i]/gradientNorm;
 
            // Now perform backtracking line search
            // const double m = sp(&unitStep[0], &prevGradient[0], nCoeffs);
            // assert(m < 0.0);
            const double m = -gradientNorm;
            double alpha = unitStepLen/armijoTau;
            bool lineIterConverged = false;
            unsigned lit = 0;
            double lineChisq = 0.0;
            for (; lit<maxLineIterations && !lineIterConverged; ++lit)
            {
                alpha *= armijoTau;
                for (unsigned i=0; i<nCoeffs; ++i)
                {
                    step[i] = unitStep[i]*alpha;
                    lineSearchCoeffs[i] = coeffs[i] + step[i];
                }
                lineChisq = nosde.chiSquare(
                    &lineSearchCoeffs[0], nCoeffs, shrinkages, nShrinkages,
                    &polyStats[0], &currentGradient[0], nullptr);
                lineIterConverged = lineChisq <= prevChisq + alpha*armijoC*m;
            }
            assert(lineIterConverged);
            // std::cout << "Made " << lit << " line iterations" << std::endl;
 
            // Recalculate the Hessian
            for (unsigned i=0; i<nCoeffs; ++i)
                coeffs[i] = lineSearchCoeffs[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_GRADIENT_DESCENT, iter, converged, false);
    }
}
 
#endif // NPSTAT_GRADIENTDESCENTFORNEYMANOSDE1D_HH_