devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Derivative of the incomplete Gamma function and the sample of a Gamma random variable.

Browse Source 
In addition to igamma(a, x), this code implements:
* igamma_der_a(a, x) = d igamma(a, x) / da -- derivative of igamma with respect to the parameter
* gamma_sample_der_alpha(alpha, sample) -- reparameterization derivative of a Gamma(alpha, 1) random variable sample with respect to the alpha parameter

The derivatives are computed by forward mode differentiation of the igamma(a, x) code. Although gamma_sample_der_alpha can be implemented via igamma_der_a, a separate function is more accurate and efficient due to analytical cancellation of some terms. All three functions are implemented by a method parameterized with "mode" that always computes the derivatives, but does not return them unless required by the mode. The compiler is expected to (and, based on benchmarks, does) skip the unnecessary computations depending on the mode.
This commit is contained in:
Michael Figurnov
2018-06-06 18:49:26 +01:00
parent e206f8d4a4
commit 4bd158fa37
12 changed files with 791 additions and 207 deletions
Download Patch File Download Diff File Expand all files Collapse all files

54
unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h
View File

@@ -41,6 +41,60 @@ struct functor_traits<scalar_igamma_op<Scalar> > {
};
};
/** \internal
* \brief Template functor to compute the derivative of the incomplete gamma
* function igamma_der_a(a, x)
*
* \sa class CwiseBinaryOp, Cwise::igamma_der_a
*/
template <typename Scalar>
struct scalar_igamma_der_a_op {
EIGEN_EMPTY_STRUCT_CTOR(scalar_igamma_der_a_op)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& a, const Scalar& x) const {
using numext::igamma_der_a;
return igamma_der_a(a, x);
}
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& x) const {
return internal::pigamma_der_a(a, x);
}
};
template <typename Scalar>
struct functor_traits<scalar_igamma_der_a_op<Scalar> > {
enum {
// 2x the cost of igamma
Cost = 40 * NumTraits<Scalar>::MulCost + 20 * NumTraits<Scalar>::AddCost,
PacketAccess = packet_traits<Scalar>::HasIGammaDerA
};
};
/** \internal
* \brief Template functor to compute the derivative of the sample
* of a Gamma(alpha, 1) random variable with respect to the parameter alpha
* gamma_sample_der_alpha(alpha, sample)
*
* \sa class CwiseBinaryOp, Cwise::gamma_sample_der_alpha
*/
template <typename Scalar>
struct scalar_gamma_sample_der_alpha_op {
EIGEN_EMPTY_STRUCT_CTOR(scalar_gamma_sample_der_alpha_op)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& alpha, const Scalar& sample) const {
using numext::gamma_sample_der_alpha;
return gamma_sample_der_alpha(alpha, sample);
}
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& alpha, const Packet& sample) const {
return internal::pgamma_sample_der_alpha(alpha, sample);
}
};
template <typename Scalar>
struct functor_traits<scalar_gamma_sample_der_alpha_op<Scalar> > {
enum {
// 2x the cost of igamma, minus the lgamma cost (the lgamma cancels out)
Cost = 30 * NumTraits<Scalar>::MulCost + 15 * NumTraits<Scalar>::AddCost,
PacketAccess = packet_traits<Scalar>::HasGammaSampleDerAlpha
};
};
/** \internal
* \brief Template functor to compute the complementary incomplete gamma function igammac(a, x)