# Copyright 2021 the GPflow authors.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
from collections import namedtuple
from typing import NamedTuple, Union
import tensorflow as tf
from ..base import InputData, MeanAndVariance, Parameter, RegressionData
from ..config import default_float, default_int
from ..covariances import Kuf
from ..utilities import add_noise_cov, to_default_float
from .sgpr import SGPR
[docs]class CGLB(SGPR):
"""
Conjugate Gradient Lower Bound. The key reference is
::
@InProceedings{pmlr-v139-artemev21a,
title = {Tighter Bounds on the Log Marginal Likelihood of
Gaussian Process Regression Using Conjugate Gradients},
author = {Artemev, Artem and Burt, David R. and van der Wilk, Mark},
booktitle = {Proceedings of the 38th International Conference on Machine Learning},
pages = {362--372},
year = {2021}
}
:param cg_tolerance: Determines accuracy to which conjugate
gradient is run when evaluating the elbo. Running more
iterations of CG would increase the ELBO by at most
`cg_tolerance`.
:param max_cg_iters: Maximum number of iterations of CG to run
per evaluation of the ELBO (or mean prediction).
:param restart_cg_iters: How frequently to restart the CG iteration.
Can be useful to avoid build up of numerical errors when
many steps of CG are run.
:param v_grad_optimization: If False, in every evaluation of the
ELBO, CG is run to select a new auxilary vector `v`. If
False, no CG is run when evaluating the ELBO but
gradients with respect to `v` are tracked so that it can
be optimized jointly with other parameters.
"""
def __init__(
self,
data: RegressionData,
*args,
cg_tolerance: float = 1.0,
max_cg_iters: int = 100,
restart_cg_iters: int = 40,
v_grad_optimization: bool = False,
**kwargs,
):
super(CGLB, self).__init__(data, *args, **kwargs)
n, b = self.data[1].shape
self._v = Parameter(tf.zeros((b, n), dtype=default_float()), trainable=v_grad_optimization)
self._cg_tolerance = cg_tolerance
self._max_cg_iters = max_cg_iters
self._restart_cg_iters = restart_cg_iters
@property
def aux_vec(self):
return self._v
[docs] def logdet_term(self, common: NamedTuple) -> tf.Tensor:
"""
Compute a lower bound on -0.5 * log |K + σ²I| based on a
low-rank approximation to K.
.. math::
log |K + σ²I| <= log |Q + σ²I| + n * log(1 + tr(K - Q)/(σ²n)).
This bound is at least as tight as
.. math::
log |K + σ²I| <= log |Q + σ²I| + tr(K - Q)/σ²,
which appears in SGPR.
"""
LB = common.LB
AAT = common.AAT
x, y = self.data
num_data = to_default_float(tf.shape(y)[0])
output_dim = to_default_float(tf.shape(y)[1])
sigma_sq = self.likelihood.variance
kdiag = self.kernel(x, full_cov=False)
# t / σ²
trace = tf.reduce_sum(kdiag) / sigma_sq - tf.reduce_sum(tf.linalg.diag_part(AAT))
logdet_b = tf.reduce_sum(tf.math.log(tf.linalg.diag_part(LB)))
logsigma_sq = num_data * tf.math.log(sigma_sq)
# Correction term from Jensen's inequality
logtrace = num_data * tf.math.log(1 + trace / num_data)
return -output_dim * (logdet_b + 0.5 * logsigma_sq + 0.5 * logtrace)
[docs] def quad_term(self, common: NamedTuple) -> tf.Tensor:
"""
Computes a lower bound on the quadratic term in the log
marginal likelihood of conjugate GPR. The bound is based on
an auxiliary vector, v. For :math:`Q ≺ K` and :math:`r=y - Kv`
.. math::
-0.5 * (rᵀQ⁻¹r + 2yᵀv - vᵀ K v ) <= -0.5 * yᵀK⁻¹y <= -0.5 * (2yᵀv - vᵀKv).
Equality holds if :math:`r=0`, i.e. :math:`v = K⁻¹y`.
If `self.aux_vec` is trainable, gradients are computed with
respect to :math:`v` as well and :math:`v` can be optimized
using gradient based methods.
Otherwise, :math:`v` is updated with the method of conjugate
gradients (CG). CG is run until :math:`0.5 * rᵀQ⁻¹r <= ϵ`,
which ensures that the maximum bias due to this term is not
more than :math:`ϵ`. The :math:`ϵ` is the CG tolerance.
"""
x, y = self.data
err = y - self.mean_function(x)
sigma_sq = self.likelihood.variance
K = add_noise_cov(self.kernel.K(x), sigma_sq)
A = common.A
LB = common.LB
preconditioner = NystromPreconditioner(A, LB, sigma_sq)
err_t = tf.transpose(err)
v_init = self.aux_vec
if not v_init.trainable:
v = cglb_conjugate_gradient(
K,
err_t,
v_init,
preconditioner,
self._cg_tolerance,
self._max_cg_iters,
self._restart_cg_iters,
)
else:
v = v_init
Kv = v @ K
r = err_t - Kv
_, error_bound = preconditioner(r)
lb = tf.reduce_sum(v * (r + 0.5 * Kv))
ub = lb + 0.5 * error_bound
if not v_init.trainable:
v_init.assign(v)
return -ub
[docs] def predict_f(
self,
xnew: InputData,
full_cov=False,
full_output_cov=False,
cg_tolerance: Union[float, None] = 1e-3,
) -> MeanAndVariance:
"""
The posterior mean for CGLB model is given by
.. :math::
m(xs) = K_{sf}v + Q_{ff}Q⁻¹r
where :math:`r = y - K v` is the residual from CG.
Note that when :math:`v=0`, this agree with the SGPR mean,
while if :math:`v = K⁻¹ y`, then :math:`r=0`, and the exact
GP mean is recovered.
:param cg_tolerance: float or None: If None, the cached value of
:math:`v` is used. If float, conjugate gradient is run until :math:`rᵀQ⁻¹r < ϵ`.
"""
x, y = self.data
err = y - self.mean_function(x)
kxx = self.kernel(x, x)
ksf = self.kernel(xnew, x)
sigma_sq = self.likelihood.variance
sigma = tf.sqrt(sigma_sq)
iv = self.inducing_variable
kernel = self.kernel
matmul = tf.linalg.matmul
trisolve = tf.linalg.triangular_solve
kmat = add_noise_cov(kxx, sigma_sq)
common = self._common_calculation()
A, LB, L = common.A, common.LB, common.L
v = self.aux_vec
if cg_tolerance is not None:
preconditioner = NystromPreconditioner(A, LB, sigma_sq)
err_t = tf.transpose(err)
v = cglb_conjugate_gradient(
kmat,
err_t,
v,
preconditioner,
cg_tolerance,
self._max_cg_iters,
self._restart_cg_iters,
)
cg_mean = matmul(ksf, v, transpose_b=True)
res = err - matmul(kmat, v, transpose_b=True)
Kus = Kuf(iv, kernel, xnew)
Ares = matmul(A, res) # The god of war!
c = trisolve(LB, Ares, lower=True) / sigma
tmp1 = trisolve(L, Kus, lower=True)
tmp2 = trisolve(LB, tmp1, lower=True)
sgpr_mean = matmul(tmp2, c, transpose_a=True)
if full_cov:
var = (
kernel(xnew)
+ matmul(tmp2, tmp2, transpose_a=True)
- matmul(tmp1, tmp1, transpose_a=True)
)
var = tf.tile(var[None, ...], [self.num_latent_gps, 1, 1]) # [P, N, N]
else:
var = (
kernel(xnew, full_cov=False)
+ tf.reduce_sum(tf.square(tmp2), 0)
- tf.reduce_sum(tf.square(tmp1), 0)
)
var = tf.tile(var[:, None], [1, self.num_latent_gps])
mean = sgpr_mean + cg_mean + self.mean_function(xnew)
return mean, var
[docs] def predict_y(
self,
xnew: InputData,
full_cov: bool = False,
full_output_cov: bool = False,
cg_tolerance: Union[float, None] = 1e-3,
) -> MeanAndVariance:
"""
Compute the mean and variance of the held-out data at the
input points.
"""
if full_cov or full_output_cov:
# See https://github.com/GPflow/GPflow/issues/1461
raise NotImplementedError(
"The predict_y method currently supports only the argument values full_cov=False and full_output_cov=False"
)
f_mean, f_var = self.predict_f(
xnew, full_cov=full_cov, full_output_cov=full_output_cov, cg_tolerance=cg_tolerance
)
return self.likelihood.predict_mean_and_var(f_mean, f_var)
[docs] def predict_log_density(
self,
data: RegressionData,
full_cov: bool = False,
full_output_cov: bool = False,
cg_tolerance: Union[float, None] = 1e-3,
) -> tf.Tensor:
"""
Compute the log density of the data at the new data points.
"""
if full_cov or full_output_cov:
raise NotImplementedError(
"The predict_log_density method currently supports only the argument values full_cov=False and full_output_cov=False"
)
x, y = data
f_mean, f_var = self.predict_f(
x, full_cov=full_cov, full_output_cov=full_output_cov, cg_tolerance=cg_tolerance
)
return self.likelihood.predict_log_density(f_mean, f_var, y)
[docs]class NystromPreconditioner:
"""
Preconditioner of the form :math:`Q=(Q_ff + σ²I)⁻¹`,
where L is lower triangular with :math: `LLᵀ = Kᵤᵤ`
:math:`A = σ⁻²L⁻¹Kᵤₓ` and :math:`B = AAᵀ + I = LᵦLᵦᵀ`
"""
def __init__(self, A: tf.Tensor, LB: tf.Tensor, sigma_sq: float):
self.A = A
self.LB = LB
self.sigma_sq = sigma_sq
def __call__(self, v):
"""
Computes :math:`vᵀQ^{-1}` and `vᵀQ^{-1}v`. Note that this is
implemented as multipication of a row vector on the right.
:param v: Vector we want to backsolve. Shape [B, N].
"""
sigma_sq = self.sigma_sq
A = self.A
LB = self.LB
trans = tf.transpose
trisolve = tf.linalg.triangular_solve
matmul = tf.linalg.matmul
sum = tf.reduce_sum
v = trans(v)
Av = matmul(A, v)
LBinvAv = trisolve(LB, Av)
LBinvtLBinvAv = trisolve(trans(LB), LBinvAv, lower=False)
rv = v - matmul(A, LBinvtLBinvAv, transpose_a=True)
vtrv = sum(rv * v)
return trans(rv) / sigma_sq, vtrv / sigma_sq
[docs]def cglb_conjugate_gradient(
K, b, initial, preconditioner, cg_tolerance, max_steps, restart_cg_step
):
"""
Conjugate gradient algorithm used in CGLB model. The method of
conjugate gradient (Hestenes and Stiefel, 1952) produces a
sequence of vectors :math:`v_0, v_1, v_2, ..., v_N` such that
:math:`v_0` = initial, and (in exact arithmetic)
:math:`Kv_n = b`. In practice, the v_i often converge quickly to
approximate :math:`K^{-1}b`, and the algorithm can be stopped
without running N iterations.
We assume the preconditioner, :math:`Q`, satisfies :math:`Q ≺ K`,
and stop the algorithm when :math:`r_i = b - Kv_i` satisfies
:math:`||rᵢᵀ||_{Q⁻¹r}^2 = rᵢᵀQ⁻¹rᵢ <= ϵ`.
:param K: Matrix we want to backsolve from. Must be PSD. Shape [N, N].
:param b: Vector we want to backsolve. Shape [B, N].
:param initial: Initial vector solution. Shape [N].
:param preconditioner: Preconditioner function.
:param cg_tolerance: Expected maximum error. This value is used
as a decision boundary against stopping criteria.
:param max_steps: Maximum number of CG iterations.
:param restart_cg_step: Restart step at which the CG resets the
internal state to the initial position using the currect
solution vector :math:`v`. Can help avoid build up of
numerical errors.
:return: `v` where `v` approximately satisfies :math:`Kv = b`.
"""
CGState = namedtuple("CGState", ["i", "v", "r", "p", "rz"])
def stopping_criterion(state):
return (0.5 * state.rz > cg_tolerance) and (state.i < max_steps)
def cg_step(state: CGState):
Ap = state.p @ K
denom = tf.reduce_sum(state.p * Ap, axis=-1)
gamma = state.rz / denom
v = state.v + gamma * state.p
i = state.i + 1
r = tf.cond(
state.i % restart_cg_step == restart_cg_step - 1,
lambda: b - v @ K,
lambda: state.r - gamma * Ap,
)
z, new_rz = preconditioner(r)
p = tf.cond(
state.i % restart_cg_step == restart_cg_step - 1,
lambda: z,
lambda: z + state.p * new_rz / state.rz,
)
return [CGState(i, v, r, p, new_rz)]
Kv = initial @ K
r = b - Kv
z, rz = preconditioner(r)
p = z
i = tf.constant(0, dtype=default_int())
initial_state = CGState(i, initial, r, p, rz)
final_state = tf.while_loop(stopping_criterion, cg_step, [initial_state])
final_state = tf.nest.map_structure(tf.stop_gradient, final_state)
return final_state[0].v