calculate_score_and_gradient.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef cMHN_COMMON_CALCULATE_SCORE_AND_GRADIENT_H
#define cMHN_COMMON_CALCULATE_SCORE_AND_GRADIENT_H
 
#include <map>
#include <tuple>
 
#include <cmhn/common/calculate_pTheta.hpp>
#include <cmhn/common/regulator.hpp>
#include <cmhn/common/score.hpp>
#include <cmhn/nontt/mhn_operator.hpp>
#include <cmhn/tt/mhn_operator.hpp>
 
#include <prc.hpp>
 
namespace cMHN
{
    template<class T, pRC::Size D, class S>
    std::tuple<T, pRC::Tensor<T, D, D>>
    calculateScoreAndGradient(nonTT::MHNOperator<T, D> const &op,
        std::map<S, T> const &pD, cMHN::Score<T> const &Score,
        cMHN::Regulator<T, D> const &Regulator, T const &toleranceSolverP = 1e-4,
        T const &toleranceSolverQ = 1e-4)
    {
        T score = pRC::zero();
        pRC::Tensor<T, D, D> g = pRC::zero();
 
        // use a random pInit
        pRC::SeedSequence seq(8, 16);
        pRC::RandomEngine rng(seq);
        pRC::GaussianDistribution<pRC::Float<>> dist;
        auto pInit = eval(expand(pRC::makeConstantSequence<pRC::Size, D, 2>(),
            [&](auto const... Ns)
            {
                return pRC::random<pRC::Tensor<T, Ns...>>(rng, dist);
            }));
        pInit = pInit / norm(pInit);
        
        auto const pTheta = calculatePTheta(op, pInit, toleranceSolverP);
 
        for(auto const &[k, v] : pD)
        {
            score += Score.pointwiseScore(v, pTheta(k));
        }
 
        // calculate g = dS/dtheta
        auto rhs =
            eval(expand(pRC::makeConstantSequence<pRC::Size, D, 2>(),
                [](auto const... Ns)
                {
                    return pRC::zero<pRC::Tensor<T, Ns...>>();
                }));
 
        for(auto const &[k, v] : pD)
        {
            rhs(k) = Score.pointwiseDSDP(v, pTheta(k));
        }
 
        auto const q = pRC::solve<pRC::Solver::GMRES<>,
            pRC::Operator::Transform::Transpose>(op, rhs,
            pRC::zero<decltype(rhs)>(), toleranceSolverQ);
 
        // this follows appendix C of Rudi's Thesis
        for(pRC::Index i = 0; i < D; ++i)
        {
            auto r = hadamardProduct(q,
                nonTT::applyDerivative(op, pTheta, i));
 
            pRC::range<pRC::Context::CompileTime, D>(
                [&](auto const j)
                {
                    if(i == j)
                    {
                        g(i, j) = pRC::reduce<pRC::Add>(r)();
                    }
                    else
                    {
                        g(i, j) = pRC::reduce<pRC::Add>(pRC::chip<j>(r, 1))();
                    }
                });
        }
 
        score -= Regulator.score(log(op.theta()));
        g -= Regulator.grad(log(op.theta()));
 
        return std::make_tuple(score, g);
    }
 
    template<pRC::Size RP, pRC::Size RQ, class T, pRC::Size D, class S>
    std::tuple<T, pRC::Tensor<T, D, D>>
    calculateScoreAndGradient(TT::MHNOperator<T, D> const &op,
        std::map<S, T> const &pD, cMHN::Score<T> const &Score,
        cMHN::Regulator<T, D> const &Regulator, T const &toleranceSolverP = 1e-4,
        T const &toleranceSolverQ = 1e-4)
    {
        using ModeSizes = decltype(TT::getModeSizes<D>());
        using Ranks = decltype(TT::getRanks<D, RP>());
 
        T score = pRC::zero();
        pRC::Tensor<T, D, D> g = pRC::zero();
        
        // use a random pInit
        pRC::SeedSequence seq(8, 16);
        pRC::RandomEngine rng(seq);
        pRC::GaussianDistribution<pRC::Float<>> dist;
        auto pInit = round<Ranks>(
            pRC::random<pRC::TensorTrain::Tensor<T, ModeSizes, Ranks>>
            (rng, dist));
        pInit /= scalarProduct(pInit,
            pRC::unit<pRC::TensorTrain::Tensor<T, ModeSizes>>())();
 
        auto const pTheta = 
            calculatePTheta<RP>(op, pInit, toleranceSolverP);
 
        T scoreT = pRC::zero();
        pRC::Tensor<T, D, D> gT = pRC::zero();
#if defined(_OPENMP)
    #pragma omp declare reduction(+: T, pRC::Tensor<T, D, D>: omp_out = omp_in + omp_out) \
        initializer (omp_priv(pRC::Zero()))
    #pragma omp parallel for schedule(dynamic, 10) reduction(+ : scoreT, gT)
#endif
        for(pRC::Index s = 0; s < pD.size(); ++s)
        {
            auto it = pD.cbegin();
            std::advance(it, s);
            auto const [k, v] = *it;
            auto const pThetaE = pTheta(k);
 
            scoreT += Score.pointwiseScore(v, pThetaE);
 
            auto const rhs =
                pRC::TensorTrain::Tensor<T, ModeSizes>::Single(
                    pRC::identity<T>(), k);
 
            // solve (1-Q)T * q = rhs
            auto const q = TT::ALS<RQ, pRC::Operator::Transform::Transpose>(
                op, rhs, toleranceSolverQ);
 
            pRC::Tensor<T, D, D> tmp = pRC::zero();
            for(pRC::Index i = 0; i < D; ++i)
            {
                for(pRC::Index j = 0; j < D; ++j)
                {
                    gT(i, j) = -Score.pointwiseDSDP(v, pThetaE) *
                        scalarProduct(q, op.derivative(i, j) * pTheta)();
                }
            }
        }
 
        score += scoreT - Regulator.score(log(op.theta()));
        g += gT - Regulator.grad(log(op.theta()));
 
        return std::make_tuple(score, g);
    }
}
 
#endif