TensorFlow定义支持熵的操作

＃版权所有2016 TensorFlow作者.版权所有.
＃
＃根据Apache许可证2.0版(“许可证”)许可;
＃你不能使用这个文件,除非符合许可证.
＃您可以获得许可证的副本
＃
＃http://www.apache.org/licenses/LICENSE-2.0
＃
＃除非适用法律要求或书面同意软件
根据许可证分发的＃分发在“按原样”基础上,
＃无明示或暗示的任何形式的担保或条件.
＃查看有关权限的特定语言的许可证
＃许可证下的限制.
＃ =============================================== =============================
“”“支持熵操作,参见$ {python / contrib.bayesflow.entropy}.
@@elbo_ratio
@@entropy_shannon
@@renyi_ratio
@@renyi_alpha
"""
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
import math
from tensorflow.contrib.bayesflow.python.ops import monte_carlo_impl as monte_carlo
from tensorflow.contrib.bayesflow.python.ops import variational_inference
from tensorflow.contrib.bayesflow.python.ops.monte_carlo_impl import _get_samples as get_samples
from tensorflow.python.framework import ops
from tensorflow.python.ops import check_ops
from tensorflow.python.ops import math_ops
from tensorflow.python.platform import tf_logging as logging
# Make utility functions from monte_carlo available.
# pylint: disable=protected-access
_get_samples = get_samples
_logspace_mean = monte_carlo._logspace_mean
_sample_mean = monte_carlo._sample_mean
# pylint: enable=protected-access
__all__ = [
    'elbo_ratio',
    'entropy_shannon',
    'renyi_ratio',
    'renyi_alpha',
]
ELBOForms = variational_inference.ELBOForms  # pylint: disable=invalid-name
def elbo_ratio(log_p,
               q,
               z=None,
               n=None,
               seed=None,
               form=None,
               name='elbo_ratio'):
  r"""Estimate of the ratio appearing in the `ELBO` and `KL` divergence.
  With `p(z) := exp{log_p(z)}`, this `Op` returns an approximation of
  ```
  E_q[ Log[p(Z) / q(Z)] ]
  ```
  The term `E_q[ Log[p(Z)] ]` is always computed as a sample mean.
  The term `E_q[ Log[q(z)] ]` can be computed with samples, or an exact formula
  if `q.entropy()` is defined.  This is controlled with the kwarg `form`.
  This log-ratio appears in different contexts:
  #### `KL[q || p]`
  If `log_p(z) = Log[p(z)]` for distribution `p`, this `Op` approximates
  the negative Kullback-Leibler divergence.
  ```
  elbo_ratio(log_p, q, n=100) = -1 * KL[q || p],
  KL[q || p] = E[ Log[q(Z)] - Log[p(Z)] ]
  ```
  Note that if `p` is a `Distribution`, then
  `distributions.kl_divergence(q, p)` may be defined and available as an
  exact result.
  #### ELBO
  If `log_p(z) = Log[p(z, x)]` is the log joint of a distribution `p`, this is
  the Evidence Lower BOund (ELBO):
  ```
  ELBO ~= E[ Log[p(Z, x)] - Log[q(Z)] ]
        = Log[p(x)] - KL[q || p]
       <= Log[p(x)]
  ```
  User supplies either `Tensor` of samples `z`, or number of samples to draw `n`
  Args:
    log_p:  Callable mapping samples from `q` to `Tensors` with
      shape broadcastable to `q.batch_shape`.
      For example, `log_p` works "just like" `q.log_prob`.
    q:  `tf.contrib.distributions.Distribution`.
    z:  `Tensor` of samples from `q`, produced by `q.sample(n)` for some `n`.
    n:  Integer `Tensor`.  Number of samples to generate if `z` is not provided.
    seed:  Python integer to seed the random number generator.
    form:  Either `ELBOForms.analytic_entropy` (use formula for entropy of `q`)
      or `ELBOForms.sample` (sample estimate of entropy), or `ELBOForms.default`
      (attempt analytic entropy, fallback on sample).
      Default value is `ELBOForms.default`.
    name:  A name to give this `Op`.
  Returns:
    Scalar `Tensor` holding sample mean KL divergence.  `shape` is the batch
      shape of `q`, and `dtype` is the same as `q`.
  Raises:
    ValueError:  If `form` is not handled by this function.
  """
  form = ELBOForms.default if form is None else form
  with ops.name_scope(name, values=[n, z]):
    z = _get_samples(q, z, n, seed)
    entropy = entropy_shannon(q, z=z, form=form)
    # If log_p(z) = Log[p(z)], cross entropy = -E_q[log(p(Z))]
    negative_cross_entropy = _sample_mean(log_p(z))
    return entropy + negative_cross_entropy
def entropy_shannon(p,
                    z=None,
                    n=None,
                    seed=None,
                    form=None,
                    name='entropy_shannon'):
  r"""Monte Carlo or deterministic computation of Shannon's entropy.
  Depending on the kwarg `form`, this `Op` returns either the analytic entropy
  of the distribution `p`, or the sampled entropy:
  ```
  -n^{-1} sum_{i=1}^n p.log_prob(z_i),  where z_i ~ p,
      \approx - E_p[ Log[p(Z)] ]
      = Entropy[p]
  ```
  User supplies either `Tensor` of samples `z`, or number of samples to draw `n`
  Args:
    p:  `tf.contrib.distributions.Distribution`
    z:  `Tensor` of samples from `p`, produced by `p.sample(n)` for some `n`.
    n:  Integer `Tensor`.  Number of samples to generate if `z` is not provided.
    seed:  Python integer to seed the random number generator.
    form:  Either `ELBOForms.analytic_entropy` (use formula for entropy of `q`)
      or `ELBOForms.sample` (sample estimate of entropy), or `ELBOForms.default`
      (attempt analytic entropy, fallback on sample).
      Default value is `ELBOForms.default`.
    name:  A name to give this `Op`.
  Returns:
    A `Tensor` with same `dtype` as `p`, and shape equal to `p.batch_shape`.
  Raises:
    ValueError:  If `form` not handled by this function.
    ValueError:  If `form` is `ELBOForms.analytic_entropy` and `n` was provided.
  """
  form = ELBOForms.default if form is None else form
  if n is not None and form == ELBOForms.analytic_entropy:
    raise ValueError('If form == ELBOForms.analytic_entropy, n must be None.')
  with ops.name_scope(name, values=[n, z]):
    # Entropy: -E_p[log(p(Z))].
    entropy = None
    # Try analytic path
    if form in [ELBOForms.default, ELBOForms.analytic_entropy]:
      try:
        entropy = p.entropy()
        logging.info('Using analytic entropy(p:%s)', p)
      except NotImplementedError as e:
        if form == ELBOForms.analytic_entropy:
          raise e
    elif form != ELBOForms.sample:
      raise ValueError('ELBOForm not handled by this function: %s' % form)
    # Sample path
    if entropy is None:
      logging.info('Using sampled entropy(p:%s)', p)
      entropy = -1. * monte_carlo.expectation(
          p.log_prob, p, z=z, n=n, seed=seed)
    return entropy
def renyi_ratio(log_p, q, alpha, z=None, n=None, seed=None, name='renyi_ratio'):
  r"""Monte Carlo estimate of the ratio appearing in Renyi divergence.
  This can be used to compute the Renyi (alpha) divergence, or a log evidence
  approximation based on Renyi divergence.
  #### Definition
  With `z_i` iid samples from `q`, and `exp{log_p(z)} = p(z)`, this `Op` returns
  the (biased for finite `n`) estimate:
  ```
  (1 - alpha)^{-1} Log[ n^{-1} sum_{i=1}^n ( p(z_i) / q(z_i) )^{1 - alpha},
  \approx (1 - alpha)^{-1} Log[ E_q[ (p(Z) / q(Z))^{1 - alpha} ]  ]
  ```
  This ratio appears in different contexts:
  #### Renyi divergence
  If `log_p(z) = Log[p(z)]` is the log prob of a distribution, and
  `alpha > 0`, `alpha != 1`, this `Op` approximates `-1` times Renyi divergence:
  ```
  # Choose reasonably high n to limit bias, see below.
  renyi_ratio(log_p, q, alpha, n=100)
                  \approx -1 * D_alpha[q || p],  where
  D_alpha[q || p] := (1 - alpha)^{-1} Log E_q[(p(Z) / q(Z))^{1 - alpha}]
  ```
  The Renyi (or "alpha") divergence is non-negative and equal to zero iff
  `q = p`.  Various limits of `alpha` lead to different special case results:
  ```
  alpha       D_alpha[q || p]
  -----       ---------------
  --> 0       Log[ int_{q > 0} p(z) dz ]
  = 0.5,      -2 Log[1 - Hel^2[q || p]],  (\propto squared Hellinger distance)
  --> 1       KL[q || p]
  = 2         Log[ 1 + chi^2[q || p] ],   (\propto squared Chi-2 divergence)
  --> infty   Log[ max_z{q(z) / p(z)} ],  (min description length principle).
  ```
  See "Renyi Divergence Variational Inference", by Li and Turner.
  #### Log evidence approximation
  If `log_p(z) = Log[p(z, x)]` is the log of the joint distribution `p`, this is
  an alternative to the ELBO common in variational inference.
  ```
  L_alpha(q, p) = Log[p(x)] - D_alpha[q || p]
  ```
  If `q` and `p` have the same support, and `0 < a <= b < 1`, one can show
  `ELBO <= D_b <= D_a <= Log[p(x)]`.  Thus, this `Op` allows a smooth
  interpolation between the ELBO and the true evidence.
  #### Stability notes
  Note that when `1 - alpha` is not small, the ratio `(p(z) / q(z))^{1 - alpha}`
  is subject to underflow/overflow issues.  For that reason, it is evaluated in
  log-space after centering.  Nonetheless, infinite/NaN results may occur.  For
  that reason, one may wish to shrink `alpha` gradually.  See the `Op`
  `renyi_alpha`.  Using `float64` will also help.
  #### Bias for finite sample size
  Due to nonlinearity of the logarithm, for random variables `{X_1,...,X_n}`,
  `E[ Log[sum_{i=1}^n X_i] ] != Log[ E[sum_{i=1}^n X_i] ]`.  As a result, this
  estimate is biased for finite `n`.  For `alpha < 1`, it is non-decreasing
  with `n` (in expectation).  For example, if `n = 1`, this estimator yields the
  same result as `elbo_ratio`, and as `n` increases the expected value
  of the estimator increases.
  #### Call signature
  User supplies either `Tensor` of samples `z`, or number of samples to draw `n`
  Args:
    log_p:  Callable mapping samples from `q` to `Tensors` with
      shape broadcastable to `q.batch_shape`.
      For example, `log_p` works "just like" `q.log_prob`.
    q: `tf.contrib.distributions.Distribution`.
       `float64` `dtype` recommended.
       `log_p` and `q` should be supported on the same set.
    alpha:  `Tensor` with shape `q.batch_shape` and values not equal to 1.
    z:  `Tensor` of samples from `q`, produced by `q.sample` for some `n`.
    n:  Integer `Tensor`.  The number of samples to use if `z` is not provided.
      Note that this can be highly biased for small `n`, see docstring.
    seed:  Python integer to seed the random number generator.
    name:  A name to give this `Op`.
  Returns:
    renyi_result:  The scaled log of sample mean.  `Tensor` with `shape` equal
      to batch shape of `q`, and `dtype` = `q.dtype`.
  """
  with ops.name_scope(name, values=[alpha, n, z]):
    z = _get_samples(q, z, n, seed)
    # Evaluate sample mean in logspace.  Note that _logspace_mean will compute
    # (among other things) the mean of q.log_prob(z), which could also be
    # obtained with q.entropy().  However, DON'T use analytic entropy, because
    # that increases variance, and could result in NaN/Inf values of a sensitive
    # term.
    # log_values
    # = (1 - alpha) * ( Log p - Log q )
    log_values = (1. - alpha) * (log_p(z) - q.log_prob(z))
    # log_mean_values
    # = Log[ E[ values ] ]
    # = Log[ E[ (p / q)^{1-alpha} ] ]
    log_mean_values = _logspace_mean(log_values)
    return log_mean_values / (1. - alpha)
def renyi_alpha(step,
                decay_time,
                alpha_min,
                alpha_max=0.99999,
                name='renyi_alpha'):
  r"""Exponentially decaying `Tensor` appropriate for Renyi ratios.
  When minimizing the Renyi divergence for `0 <= alpha < 1` (or maximizing the
  Renyi equivalent of elbo) in high dimensions, it is not uncommon to experience
  `NaN` and `inf` values when `alpha` is far from `1`.
  For that reason, it is often desirable to start the optimization with `alpha`
  very close to 1, and reduce it to a final `alpha_min` according to some
  schedule.  The user may even want to optimize using `elbo_ratio` for
  some fixed time before switching to Renyi based methods.
  This `Op` returns an `alpha` decaying exponentially with step:
  ```
  s(step) = (exp{step / decay_time} - 1) / (e - 1)
  t(s) = max(0, min(s, 1)),  (smooth growth from 0 to 1)
  alpha(t) = (1 - t) alpha_min + t alpha_max
  ```
  Args:
    step:  Non-negative scalar `Tensor`.  Typically the global step or an
      offset version thereof.
    decay_time:  Positive scalar `Tensor`.
    alpha_min:  `float` or `double` `Tensor`.
      The minimal, final value of `alpha`, achieved when `step >= decay_time`
    alpha_max:  `Tensor` of same `dtype` as `alpha_min`.
      The maximal, beginning value of `alpha`, achieved when `step == 0`
    name:  A name to give this `Op`.
  Returns:
    alpha:  A `Tensor` of same `dtype` as `alpha_min`.
  """
  with ops.name_scope(name, values=[step, decay_time, alpha_min, alpha_max]):
    alpha_min = ops.convert_to_tensor(alpha_min, name='alpha_min')
    dtype = alpha_min.dtype
    alpha_max = ops.convert_to_tensor(alpha_max, dtype=dtype, name='alpha_max')
    decay_time = math_ops.cast(decay_time, dtype)
    step = math_ops.cast(step, dtype)
    check_scalars = [
        check_ops.assert_rank(step, 0, message='step must be scalar'),
        check_ops.assert_rank(
            decay_time, 0, message='decay_time must be scalar'),
        check_ops.assert_rank(alpha_min, 0, message='alpha_min must be scalar'),
        check_ops.assert_rank(alpha_max, 0, message='alpha_max must be scalar'),
    ]
    check_sign = [
        check_ops.assert_non_negative(
            step, message='step must be non-negative'),
        check_ops.assert_positive(
            decay_time, message='decay_time must be positive'),
    ]
    with ops.control_dependencies(check_scalars + check_sign):
      theta = (math_ops.exp(step / decay_time) - 1.) / (math.e - 1.)
      theta = math_ops.minimum(math_ops.maximum(theta, 0.), 1.)
      return alpha_max * (1. - theta) + alpha_min * theta