NeymanOSDE1DConvergence.hh

Go to the documentation of this file.

#ifndef NPSTAT_NEYMANOSDE1DCONVERGENCE_HH_
#define NPSTAT_NEYMANOSDE1DCONVERGENCE_HH_
 
/*!
// \file NeymanOSDE1DConvergence.hh
//
// \brief Classes for establishing convergence of the Neyman OSDE
//        optimization algorithms
//
// Author: I. Volobouev
//
// March 2023
*/
 
#include <vector>
#include <cmath>
#include <cassert>
 
#include "npstat/nm/Matrix.hh"
 
namespace npstat {
    // Convergence is established when either all
    // components of the gradient are small in
    // magnitude (useL2Distance parameter is false)
    // or the norm of the gradient is small
    // (useL2Distance parameter is true).
    //
    class GradientNeymanOSDE1DConvergenceCalc
    {
    public:
        inline explicit GradientNeymanOSDE1DConvergenceCalc(
            const double i_eps, const bool i_useL2Distance = false)
            : eps_(i_eps), useL2Distance_(i_useL2Distance)
        {
            assert(eps_ >= 0.0);
        }
 
        inline bool converged(const unsigned /* iterationNumber */,
                              const double /* prevChisq */,
                              const std::vector<double>& /* prevGradient */,
                              const std::vector<double>& /* lastStep */,
                              const double /* currentChisq */,
                              const std::vector<double>& currentGradient,
                              const Matrix<double>& /* hessian */) const
        {
            const unsigned sz = currentGradient.size();
            if (useL2Distance_)
            {
                long double sum = 0.0L;
                for (unsigned i=0; i<sz; ++i)
                    sum += currentGradient[i]*currentGradient[i];
                return sum <= sz*eps_*eps_;
            }
            else
            {
                for (unsigned i=0; i<sz; ++i)
                    if (std::abs(currentGradient[i]) > eps_)
                        return false;
                return true;
            }
        }
 
    private:
        double eps_;
        bool useL2Distance_;
    };
 
    // Convergence is established if either the
    // convergence conditions are satisfied for the
    // argument of the other convergence calculator
    // or the Hessian becomes positive-definite.
    //
    // Intended for automatic switching between
    // gradient descent and Newton's method.
    //
    template<class OtherConvergenceCalc>
    class HessianNeymanOSDE1DConvergenceCalc
    {
    public:        
        inline explicit HessianNeymanOSDE1DConvergenceCalc(
            const OtherConvergenceCalc& gcalc)
            : gcalc_(gcalc), gcalcConverged_(false) {}
 
        inline bool otherCalcConverged() const {return gcalcConverged_;}
 
        inline bool converged(const unsigned iterationNumber,
                              const double prevChisq,
                              const std::vector<double>& prevGradient,
                              const std::vector<double>& lastStep,
                              const double currentChisq,
                              const std::vector<double>& currentGradient,
                              const Matrix<double>& hessian) const
        {
            gcalcConverged_ = gcalc_.converged(iterationNumber, prevChisq,
                                               prevGradient, lastStep,
                                               currentChisq, currentGradient,
                                               hessian);
            if (gcalcConverged_)
                return true;
 
            assert(hessian.isSquare());
            const unsigned sz = hessian.nRows();
            std::vector<double> eigenvalues(sz);
            hessian.symEigen(&eigenvalues[0], sz);
            for (unsigned i=0; i<sz; ++i)
                if (eigenvalues[i] <= 0.0)
                    return false;
            return true;
        }
 
    private:
        OtherConvergenceCalc gcalc_;
        mutable bool gcalcConverged_;
    };
}
 
#endif // NPSTAT_NEYMANOSDE1DCONVERGENCE_HH_