Source code for neural_tangents.predict

# Copyright 2019 Google LLC
#
# 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
#
#     https://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.

"""Functions to make predictions on the train/test set using NTK/NNGP.

Most functions in this module accept training data as inputs and return a new
function `predict_fn` that computes predictions on the train set / given test
set / timesteps.

WARNING: `trace_axes` parameter supplied to prediction functions must match the
respective parameter supplied to the function used to compute the kernel.
Namely, this is the same `trace_axes` used to compute the empirical kernel
(`utils/empirical.py`; `diagonal_axes` must be `()`), or `channel_axis` in the
output of the top layer used to compute the closed-form kernel (`stax.py`; note
that closed-form kernels currently only support a single `channel_axis`).
"""


import collections
from functools import lru_cache
from typing import Union, Tuple, Callable, Iterable, Optional, Dict, NamedTuple, Sequence, Generator

import jax
from jax import grad
from jax.experimental import ode
from jax.lib import xla_bridge
import jax.numpy as np
import jax.scipy as sp
from jax.tree_util import tree_map, tree_all
from neural_tangents.utils import utils, dataclasses
from neural_tangents.utils.typing import KernelFn, Axes, Get
import numpy as onp
import scipy as osp


"""Alias for optional arrays or scalars."""
ArrayOrScalar = Union[None, int, float, np.ndarray]


[docs]def gradient_descent_mse( k_train_train: np.ndarray, y_train: np.ndarray, learning_rate: float = 1., diag_reg: float = 0., diag_reg_absolute_scale: bool = False, trace_axes: Axes = (-1,) ) -> Callable[ [ArrayOrScalar, ArrayOrScalar, ArrayOrScalar, Optional[np.ndarray]], Union[np.ndarray, Tuple[np.ndarray, np.ndarray]] ]: r"""Predicts the outcome of function space gradient descent training on MSE. Solves in closed form for the continuous-time version of gradient descent. Uses the closed-form solution for gradient descent on an MSE loss in function space detailed in [*,**] given a Neural Tangent or Neural Network Gaussian Process Kernel over the dataset. Given NNGP or NTK, this function will return a function that predicts the time evolution for function space points at arbitrary time[s] (training step[s]) `t`. Note that these time[s] (step[s]) are continuous and are interpreted in units of the `learning_rate` so `absolute_time = learning_rate * t`, and the scales of `learning_rate` and `t` are interchangeable. Note that first invocation of the returned `predict_fn` will be slow and allocate a lot of memory for its whole lifetime, as either eigendecomposition (`t` is a scalar or an array) or Cholesky factorization (`t=None`) of `k_train_train` is performed and cached for future invocations (or both, if the function is called on both finite and infinite (`t=None`) times). [*] https://arxiv.org/abs/1806.07572 [**] https://arxiv.org/abs/1902.06720 Example: >>> from neural_tangents import empirical_ntk_fn >>> from neural_tangents import predict >>> >>> t = 1e-7 >>> kernel_fn = empirical_ntk_fn(f) >>> k_train_train = kernel_fn(x_train, None, params) >>> k_test_train = kernel_fn(x_test, x_train, params) >>> >>> predict_fn = predict.gradient_descent_mse(k_train_train, y_train) >>> >>> fx_train_0 = f(params, x_train) >>> fx_test_0 = f(params, x_test) >>> >>> fx_train_t, fx_test_t = predict_fn(t, fx_train_0, fx_test_0, >>> k_test_train) Args: k_train_train: kernel on the training data. Must have the shape of `zip(y_train.shape, y_train.shape)` with `trace_axes` absent. y_train: targets for the training data. learning_rate: learning rate, step size. diag_reg: a scalar representing the strength of the diagonal regularization for `k_train_train`, i.e. computing `k_train_train + diag_reg * I` during Cholesky factorization or eigendecomposition. diag_reg_absolute_scale: `True` for `diag_reg` to represent regularization in absolute units, `False` to be `diag_reg * np.mean(np.trace(k_train_train))`. trace_axes: `f(x_train)` axes such that `k_train_train` lacks these pairs of dimensions and is to be interpreted as :math:`\Theta \otimes I`, i.e. block-diagonal along `trace_axes`. These can can be specified either to save space and compute, or to even improve approximation accuracy of the infinite-width or infinite-samples limit, since in in these limits the covariance along channel / feature / logit axes indeed converges to a constant-diagonal matrix. However, if you target linearized dynamics of a specific finite-width network, `trace_axes=()` will yield most accurate result. Returns: A function of signature `predict_fn(t, fx_train_0, fx_test_0, k_test_train)` that returns output train [and test] set[s] predictions at time[s] `t`. """ _, odd, first, _ = _get_axes(k_train_train) trace_axes = utils.canonicalize_axis(trace_axes, y_train) trace_axes = tuple(-y_train.ndim + a for a in trace_axes) n_t_axes, n_non_t_axes = len(trace_axes), y_train.ndim - len(trace_axes) last_t_axes = tuple(range(-n_t_axes, 0)) non_t_axes = tuple(range(-y_train.ndim, -n_t_axes)) @lru_cache(1) def get_predict_fn_inf(): with jax.core.eval_context(): solve = _get_cho_solve(k_train_train, diag_reg, diag_reg_absolute_scale) def predict_fn_inf(fx_train_0, fx_test_0, k_test_train): fx_train_t = y_train.astype(k_train_train.dtype) if fx_test_0 is None: return fx_train_t rhs = y_train if fx_train_0 is None else y_train - fx_train_0 dfx_test = np.tensordot(k_test_train, solve(rhs, trace_axes), (odd, first)) dfx_test = np.moveaxis(dfx_test, last_t_axes, trace_axes) fx_test_t = fx_test_0 + dfx_test if fx_train_0 is None: return fx_test_t return fx_train_t, fx_test_t return predict_fn_inf @lru_cache(1) def get_predict_fn_finite(): with jax.core.eval_context(): expm1_fn, inv_expm1_fn = _get_fns_in_eigenbasis( k_train_train, diag_reg, diag_reg_absolute_scale, (_make_expm1_fn(y_train.size), _make_inv_expm1_fn(y_train.size)) ) rhs_shape = tuple(y_train.shape[a] for a in trace_axes) def predict_fn_finite(t, fx_train_0, fx_test_0, k_test_train): t = np.array(t) * learning_rate t_shape, t_ndim = t.shape, t.ndim first_t_axes = tuple(range(t_ndim)) t = t.reshape((-1, 1)) rhs = -y_train if fx_train_0 is None else fx_train_0 - y_train rhs = np.moveaxis(rhs, trace_axes, last_t_axes).reshape( (-1,) + rhs_shape) shape = t_shape + k_train_train.shape[1::2] + rhs_shape if fx_train_0 is not None: dfx_train = expm1_fn(rhs, t).reshape(shape) dfx_train = np.moveaxis(dfx_train, last_t_axes, trace_axes) fx_train_t = np.expand_dims(fx_train_0, first_t_axes) + dfx_train if fx_test_0 is not None: dfx_test = inv_expm1_fn(rhs, t).reshape(shape) dfx_test = np.tensordot(k_test_train, dfx_test, (odd, non_t_axes)) dfx_test = np.moveaxis( dfx_test, tuple(range(n_non_t_axes, n_non_t_axes + t_ndim)) + last_t_axes, tuple(range(t_ndim)) + trace_axes) fx_test_t = np.expand_dims(fx_test_0, first_t_axes) + dfx_test if fx_train_0 is not None and fx_test_0 is not None: return fx_train_t, fx_test_t if fx_test_0 is None: return fx_train_t return fx_test_t return predict_fn_finite def predict_fn( t: Optional[ArrayOrScalar] = None, fx_train_0: ArrayOrScalar = 0., fx_test_0: Optional[ArrayOrScalar] = None, k_test_train: Optional[np.ndarray] = None ) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]: """Return output predictions on train [and test] set[s] at time[s] `t`. Args: t: a scalar of array of scalars of any shape. `t=None` is treated as infinity and returns the same result as `t=np.inf`, but is computed using identity or linear solve for train and test predictions respectively instead of eigendecomposition, saving time and precision. Equivalent of training steps (but can be fractional). fx_train_0: output of the network at `t == 0` on the training set. `fx_train_0=None` means to not compute predictions on the training set. fx_test_0: output of the network at `t == 0` on the test set. `fx_test_0=None` means to not compute predictions on the test set. k_test_train: kernel relating test data with training data. Must have the shape of `zip(y_test.shape, y_train.shape)` with `trace_axes` absent. Pass `k_test_train=None` if you only need non-regularized (`diag_reg=0`) predictions on the training set. For regularized train-set predictions, pass `k_test_train=k_train_train`. Returns: `fx_train_t` or `(fx_train_t, fx_test_t)` if `fx_test_0 != None` with potentially additional leading time dimensions matching `t.shape`. Raises: ValueError: if `fx_test_0` is not `None`, but `k_test_train` is `None`. """ _check_inputs(fx_train_0, fx_test_0, k_test_train) # Infinite time if t is None: return get_predict_fn_inf()(fx_train_0, fx_test_0, k_test_train) # Finite time return get_predict_fn_finite()(t, fx_train_0, fx_test_0, k_test_train) return predict_fn
[docs]@dataclasses.dataclass class ODEState: """ODE state dataclass holding outputs and auxiliary variables.""" fx_train: Optional[np.ndarray] = None fx_test: Optional[np.ndarray] = None qx_train: Optional[np.ndarray] = None qx_test: Optional[np.ndarray] = None
[docs]def gradient_descent( loss: Callable[[np.ndarray, np.ndarray], float], k_train_train: np.ndarray, y_train: np.ndarray, learning_rate: float = 1., momentum: Optional[float] = None, trace_axes: Axes = (-1,) ) -> Callable[ [ArrayOrScalar, Union[ArrayOrScalar, ODEState], ArrayOrScalar, Optional[np.ndarray]], Union[np.ndarray, Tuple[np.ndarray, np.ndarray], ODEState] ]: r"""Predicts the outcome of function space training using gradient descent. Uses an ODE solver. If `momentum != None`, solves a continuous-time version of gradient descent with momentum (note: this case uses standard momentum as opposed to Nesterov momentum). Solves the function space ODE for [momentum] gradient descent with a given `loss` (detailed in [*]) given a Neural Tangent Kernel[s] over the dataset[s] at arbitrary time[s] (step[s]) `t`. Note that for gradient descent `absolute_time = learning_rate * t` and the scales of the learning rate and query step[s] `t` are interchangeable. However, the momentum gradient descent ODE is solved in the units of `learning_rate**0.5`, and therefore `absolute_time = learning_rate**0.5 * t`, hence the `learning_rate` and training time[s] (step[s]) `t` scales are not interchangeable. [*] https://arxiv.org/abs/1902.06720 Example: >>> from neural_tangents import empirical_ntk_fn >>> from neural_tangents import predict >>> >>> t = 1e-7 >>> learning_rate = 1e-2 >>> momentum = 0.9 >>> >>> kernel_fn = empirical_ntk_fn(f) >>> k_test_train = kernel_fn(x_test, x_train, params) >>> >>> from jax.experimental import stax >>> cross_entropy = lambda fx, y_hat: -np.mean(stax.logsoftmax(fx) * y_hat) >>> predict_fn = predict.gradient_descent(cross_entropy, k_train_train, >>> y_train, learning_rate, momentum) >>> >>> fx_train_0 = f(params, x_train) >>> fx_test_0 = f(params, x_test) >>> >>> fx_train_t, fx_test_t = predict_fn(t, fx_train_0, fx_test_0, >>> k_test_train) Args: loss: a loss function whose signature is `loss(f(x_train), y_train)`. Note: the loss function should treat the batch and output dimensions symmetrically. k_train_train: kernel on the training data. Must have the shape of `zip(y_train.shape, y_train.shape)` with `trace_axes` absent. y_train: targets for the training data. learning_rate: learning rate, step size. momentum: momentum scalar. trace_axes: `f(x_train)` axes such that `k_train_train` lacks these pairs of dimensions and is to be interpreted as :math:`\Theta \otimes I`, i.e. block-diagonal along `trace_axes`. These can can be specified either to save space and compute, or to even improve approximation accuracy of the infinite-width or infinite-samples limit, since in in these limits the covariance along channel / feature / logit axes indeed converges to a constant-diagonal matrix. However, if you target linearized dynamics of a specific finite-width network, `trace_axes=()` will yield most accurate result. Returns: A function that returns output train [and test] set[s] predictions at time[s] `t`. """ _, odd, _, _ = _get_axes(k_train_train) trace_axes = utils.canonicalize_axis(trace_axes, y_train) non_t_axes = tuple(a for a in range(y_train.ndim) if a not in trace_axes) last_t_axes = range(-len(trace_axes), 0) dtype = k_train_train.dtype grad_loss = grad(lambda fx: loss(fx, y_train)) if momentum is not None: learning_rate **= 0.5 momentum = (momentum - 1.0) / learning_rate def get_state_0(fx_train_or_state_0, fx_test_0, fx_test_shape): if isinstance(fx_train_or_state_0, ODEState): fx_train_0 = fx_train_or_state_0.fx_train fx_test_0 = fx_train_or_state_0.fx_test qx_train_0 = fx_train_or_state_0.qx_train qx_test_0 = fx_train_or_state_0.qx_test else: fx_train_0 = fx_train_or_state_0 qx_train_0 = qx_test_0 = None if fx_train_0 is None: fx_train_0 = np.zeros_like(y_train, dtype) else: fx_train_0 = np.broadcast_to(fx_train_0, y_train.shape) if fx_test_0 is not None: fx_test_0 = np.broadcast_to(fx_test_0, fx_test_shape) if momentum is None: if qx_train_0 is not None or qx_test_0 is not None: raise ValueError('Got passed momentum state variables, while ' '`momentum is None`.') else: qx_train_0 = (np.zeros_like(y_train, dtype) if qx_train_0 is None else np.broadcast_to(qx_train_0, y_train.shape)) qx_test_0 = (None if fx_test_0 is None else (np.zeros(fx_test_shape, dtype) if qx_test_0 is None else np.broadcast_to(qx_test_0, fx_test_shape))) return ODEState(fx_train_0, fx_test_0, qx_train_0, qx_test_0) # pytype: disable=wrong-arg-count def get_dstate_dt(k_test_train): def dstate_dt(state_t: ODEState, unused_t) -> ODEState: fx_train_t, fx_test_t, qx_train_t, qx_test_t = ( state_t.fx_train, state_t.fx_test, state_t.qx_train, state_t.qx_test) dy_df_t = grad_loss(fx_train_t) fx_train_t = -np.moveaxis( np.tensordot(k_train_train, dy_df_t, (odd, non_t_axes)), last_t_axes, trace_axes ) if fx_test_t is not None: fx_test_t = -np.moveaxis( np.tensordot(k_test_train, dy_df_t, (odd, non_t_axes)), last_t_axes, trace_axes ) if momentum is None: return ODEState(fx_train_t, fx_test_t) # pytype: disable=wrong-arg-count fx_train_t += momentum * qx_train_t if qx_test_t is not None: fx_test_t += momentum * qx_test_t return ODEState(qx_train_t, qx_test_t, fx_train_t, fx_test_t) # pytype: disable=wrong-arg-count return dstate_dt def predict_fn( t: Optional[ArrayOrScalar] = None, fx_train_or_state_0: Union[ArrayOrScalar, ODEState] = 0., fx_test_0: Optional[ArrayOrScalar] = None, k_test_train: Optional[np.ndarray] = None ) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray], ODEState]: """Return output predictions on train [and test] set[s] at time[s] `t`. Args: t: a scalar or array of scalars of any shape in strictly increasing order. `t=None` is equivalent to `t=np.inf` and may not converge. Equivalent of training steps (but can be fractional). fx_train_or_state_0: either (a) output of the network at `t == 0` on the training set or (b) complete ODE state (`predict.ODEState`). Pass an ODE state if you want to operate on the full ODE state instead of output variables only (useful for inspecting auxiliary variables or resuming an optimizer with auxiliary variables from a specific state. Note that only `momentum != None` optimizer currently has auxiliary variables. To initialize an ODE state from scratch, call `predict.ODEState(fx_train_0, fx_test_0)`. If an ODE state is passed, an ODE state is returned. `fx_train_0=None` means to not compute predictions on the training set. fx_test_0: output of the network at `t == 0` on the test set. `fx_test_0=None` means to not compute predictions on the test set. k_test_train: kernel relating test data with training data. Must have the shape of `zip(y_test.shape, y_train.shape)` with `trace_axes` absent. Pass `k_test_train=None` if you only need predictions on the training set. Returns: `fx_train_t` or `(fx_train_t, fx_test_t)` if `fx_test_0 != None` with potentially additional leading time dimensions matching `t.shape`. Alternatively can return an `ODEState` at time[s] `t`. Raises: ValueError: if `fx_test_0` is not `None`, but `k_test_train` is `None`. """ _check_inputs(fx_train_or_state_0, fx_test_0, k_test_train) t = np.array(t if t is not None else np.inf, dtype) * learning_rate t_shape = t.shape t = t.reshape((-1,)) # ODE solver requires `t[0]` to be the time where `fx_train_0` [and # `fx_test_0`] are evaluated, but also a strictly increasing sequence of # timesteps, so we always temporarily append an [almost] `0` at the start. t0 = np.where(t[0] == 0, np.full((1,), -1e-24, t.dtype), np.zeros((1,), t.dtype)) t = np.concatenate([t0, t]) # Solve the ODE. fx_test_shape = _get_fx_test_shape(y_train, k_test_train, trace_axes) state_0 = get_state_0(fx_train_or_state_0, fx_test_0, fx_test_shape) state_t = ode.odeint(get_dstate_dt(k_test_train), state_0, t) # Remove the added `t0`. trim = lambda x: x[1:].reshape(t_shape + x.shape[1:]) trim_tree = lambda tree: tree_map(trim, tree) state_t = trim_tree(state_t) # `ODEState` -> `ODEState` if isinstance(fx_train_or_state_0, ODEState): return state_t # `np.ndarray` -> `np.ndarray` fx_train_t, fx_test_t = state_t.fx_train, state_t.fx_test if fx_train_or_state_0 is not None and fx_test_0 is None: return fx_train_t if fx_test_0 is not None and fx_train_or_state_0 is None: return fx_test_t return fx_train_t, fx_test_t return predict_fn
[docs]class Gaussian(NamedTuple): """A `(mean, covariance)` convenience namedtuple.""" mean: np.ndarray covariance: np.ndarray
[docs]def gp_inference( k_train_train, y_train: np.ndarray, diag_reg: float = 0., diag_reg_absolute_scale: bool = False, trace_axes: Axes = (-1,)): r"""Compute the mean and variance of the 'posterior' of NNGP/NTK/NTKGP. NNGP - the exact posterior of an infinitely wide Bayesian NN. NTK - exact distribution of an infinite ensemble of infinitely wide NNs trained with gradient flow for infinite time. NTKGP - posterior of a GP (Gaussian process) with the NTK covariance (see https://arxiv.org/abs/2007.05864 for how this can correspond to infinite ensembles of infinitely wide NNs as well). Note that first invocation of the returned `predict_fn` will be slow and allocate a lot of memory for its whole lifetime, as a Cholesky factorization of `k_train_train.nngp` or `k_train_train.ntk` (or both) is performed and cached for future invocations. Args: k_train_train: train-train kernel. Can be (a) `np.ndarray`, (b) `Kernel` namedtuple, (c) `Kernel` object. Must contain the necessary `nngp` and/or `ntk` kernels for arguments provided to the returned `predict_fn` function. For example, if you request to compute posterior test [only] NTK covariance in future `predict_fn` invocations, `k_train_train` must contain both `ntk` and `nngp` kernels. y_train: train targets. diag_reg: a scalar representing the strength of the diagonal regularization for `k_train_train`, i.e. computing `k_train_train + diag_reg * I` during Cholesky factorization. diag_reg_absolute_scale: `True` for `diag_reg` to represent regularization in absolute units, `False` to be `diag_reg * np.mean(np.trace(k_train_train))`. trace_axes: `f(x_train)` axes such that `k_train_train`, `k_test_train`[, and `k_test_test`] lack these pairs of dimensions and are to be interpreted as :math:`\Theta \otimes I`, i.e. block-diagonal along `trace_axes`. These can can be specified either to save space and compute, or to even improve approximation accuracy of the infinite-width or infinite-samples limit, since in in these limits the covariance along channel / feature / logit axes indeed converges to a constant-diagonal matrix. However, if you target linearized dynamics of a specific finite-width network, `trace_axes=()` will yield most accurate result. Returns: A function of signature `predict_fn(get, k_test_train, k_test_test)` computing 'posterior' Gaussian distribution (mean or mean and covariance) on a given test set. """ even, odd, first, last = _get_axes(_get_first(k_train_train)) trace_axes = utils.canonicalize_axis(trace_axes, y_train) @lru_cache(2) def solve(g: str): k_dd = _get_attr(k_train_train, g) return _get_cho_solve(k_dd, diag_reg, diag_reg_absolute_scale) @lru_cache(2) def k_inv_y(g: str): return solve(g)(y_train, trace_axes) @utils.get_namedtuple('Gaussians') def predict_fn(get: Optional[Get] = None, k_test_train=None, k_test_test=None ) -> Dict[str, Union[np.ndarray, Gaussian]]: """`test`-set posterior given respective covariance matrices. Args: get: string, the mode of the Gaussian process, either "nngp", "ntk", "ntkgp", (see https://arxiv.org/abs/2007.05864) or a tuple, or `None`. If `None` then both `nngp` and `ntk` predictions are returned. k_test_train: test-train kernel. Can be (a) `np.ndarray`, (b) `Kernel` namedtuple, (c) `Kernel` object. Must contain the necessary `nngp` and/or `ntk` kernels for arguments provided to the returned `predict_fn` function. For example, if you request to compute posterior test [only] NTK covariance, `k_test_train` must contain both `ntk` and `nngp` kernels. If `None`, returns predictions on the training set. Note that train-set outputs are always `N(y_train, 0)` and mostly returned for API consistency. k_test_test: test-test kernel. Can be (a) `np.ndarray`, (b) `Kernel` namedtuple, (c) `Kernel` object. Must contain the necessary `nngp` and/or `ntk` kernels for arguments provided to the returned `predict_fn` function. Provide if you want to compute test-test posterior covariance. `k_test_test=None` means to not compute it. If `k_test_train is None`, pass any non-`None` value (e.g. `True`) if you want to get non-regularized (`diag_reg=0`) train-train posterior covariance. Note that non-regularized train-set outputs will always be the zero-variance Gaussian `N(y_train, 0)` and mostly returned for API consistency. For regularized train-set posterior outputs according to a positive `diag_reg`, pass `k_test_train=k_train_train`, and, optionally, `k_test_test=nngp_train_train`. Returns: Either a `Gaussian('mean', 'variance')` namedtuple or `mean` of the GP posterior on the `test` set. """ if get is None: get = ('nngp', 'ntk') out = {} for g in get: k = g if g != 'ntkgp' else 'ntk' k_dd = _get_attr(k_train_train, k) k_td = None if k_test_train is None else _get_attr(k_test_train, k) if k_td is None: # Train set predictions. y = y_train.astype(k_dd.dtype) else: # Test set predictions. y = np.tensordot(k_td, k_inv_y(k), (odd, first)) y = np.moveaxis(y, range(-len(trace_axes), 0), trace_axes) if k_test_test is not None: if k_td is None: out[g] = Gaussian(y, np.zeros_like(k_dd, k_dd.dtype)) else: if (g == 'ntk' and (not hasattr(k_train_train, 'nngp') or not hasattr(k_test_train, 'nngp'))): raise ValueError( 'If `"ntk" in get`, and `k_test_test is not None`, ' 'and `k_test_train is not None`, i.e. you request the ' 'NTK posterior covariance on the test set, you need ' 'both NTK and NNGP train-train and test-train matrices ' 'contained in `k_test_train` and `k_train_train`. ' 'Hence they must be `namedtuple`s with `nngp` and ' '`ntk` attributes.') # kernel of wide NN at initialization g_init = 'nngp' if g != 'ntkgp' else 'ntk' k_td_g_inv_y = solve(k)(_get_attr(k_test_train, g_init), even) k_tt = _get_attr(k_test_test, g_init) if g == 'nngp' or g == 'ntkgp': cov = np.tensordot(k_td, k_td_g_inv_y, (odd, first)) cov = k_tt - utils.zip_axes(cov) out[g] = Gaussian(y, cov) elif g == 'ntk': term_1 = solve(g)(k_td, even) cov = np.tensordot(_get_attr(k_train_train, 'nngp'), term_1, (odd, first)) cov = np.tensordot(term_1, cov, (first, first)) term_2 = np.tensordot(k_td, k_td_g_inv_y, (odd, first)) term_2 += np.moveaxis(term_2, first, last) cov = utils.zip_axes(cov - term_2) + k_tt out[g] = Gaussian(y, cov) else: raise ValueError(g) else: out[g] = y return out return predict_fn
"""Helper type to fit cache dictionaries to `get` API.""" _Kernel = collections.namedtuple('Kernel', 'nngp ntk') _Kernel.__new__.__defaults__ = (None,) * len(_Kernel._fields)
[docs]def gradient_descent_mse_ensemble( kernel_fn: KernelFn, x_train: np.ndarray, y_train: np.ndarray, learning_rate: float = 1., diag_reg: float = 0.0, diag_reg_absolute_scale: bool = False, trace_axes: Axes = (-1,), **kernel_fn_train_train_kwargs): r"""Predicts the gaussian embedding induced by gradient descent on MSE loss. This is equivalent to an infinite ensemble of infinite-width networks after marginalizing out the initialization, if `kernel_fn` is the kernel function of the infinite-width network. Note that `kernel_fn` can in principle also be an empirical / Monte Carlo finite-width kernel function, but in this case the returned output will not have a simple interpretation (unless these functions are used to approximate the infinite-width kernel). Note that first invocation of the returned `predict_fn` will be slow and allocate a lot of memory for its whole lifetime, as the kernel computation, and either eigendecomposition (`t` is a scalar or an array) or Cholesky factorization (`t=None`) of `kernel_fn(x_train, None, get)` is performed and cached for future invocations (or both, if the function is called on both finite and infinite (`t=None`) times). Args: kernel_fn: A kernel function that computes NNGP and/or NTK. Must have a signature `kernel_fn(x1, x2, get, **kernel_fn_kwargs)` and return a `Kernel` object or a `namedtuple` with `nngp` and/or `ntk` attributes. Therefore, it can be an `AnalyticKernelFn`, but also a `MonteCarloKernelFn`, or an `EmpiricalKernelFn` (but only `nt.empirical_kernel_fn` and not `nt.empirical_ntk_fn` or `ntk.empirical_nngp_fn`, since the latter two do not accept a `get` argument). Note that for meaningful outputs, the kernel function must represent or at least approximate the infinite-width kernel. x_train: training inputs. y_train: training targets. learning_rate: learning rate, step size. diag_reg: a scalar representing the strength of the diagonal regularization for `kernel_fn(x_train, None, get)`, i.e. computing `kernel_fn(x_train, None, get) + diag_reg * I` during Cholesky factorization or eigendecomposition. diag_reg_absolute_scale: `True` for `diag_reg` to represent regularization in absolute units, `False` to be `diag_reg * np.mean(np.trace(kernel_fn(x_train, None, get)))`. trace_axes: `f(x_train)` axes such that `kernel_fn(x_train, None, get)`, `kernel_fn(x_test, x_train, get)`[, and `kernel_fn(x_test, None, get)`] lack these pairs of dimensions and are to be interpreted as :math:`\Theta \otimes I`, i.e. block-diagonal along `trace_axes`. These can can be specified either to save space and compute, or to even improve approximation accuracy of the infinite-width or infinite-samples limit, since in in these limits the covariance along channel / feature / logit axes indeed converges to a constant-diagonal matrix. However, if you target linearized dynamics of a specific finite-width network, `trace_axes=()` will yield most accurate result. **kernel_fn_train_train_kwargs: optional keyword arguments passed to `kernel_fn`. For train-train kernel, these are passed to `kernel_fn` without changes. For test-test kernel, they are passed to `kernel_fn`, unless overwritten by a similar `**kernel_fn_test_test_kwargs` arguments passed to the `predict_fn` function call. Finally, for test-train kernel, values that are tuples of arrays (destined for calls of the finite-width network on training and testing data) will be tuples of values combined from `**kernel_fn_train_train_kwargs` and `**kernel_fn_test_test_kwargs`, and all other values must match. Returns: A function with signature `predict_fn(t, x_test, get, compute_cov)` returning either mean or mean and covariance of the infinite ensemble of infinite-width networks outputs on `x_test` at time[s] `t`, in the `get` regime (`"nngp"`, `"ntk"`, or `("nngp", "ntk")`). """ expm1 = _make_expm1_fn(y_train.size) inv_expm1 = _make_inv_expm1_fn(y_train.size) trace_axes = utils.canonicalize_axis(trace_axes, y_train) trace_axes = tuple(-y_train.ndim + a for a in trace_axes) n_trace_axes = len(trace_axes) last_t_axes = range(-n_trace_axes, 0) trace_shape = tuple(y_train.shape[a] for a in trace_axes) y_train_flat = np.moveaxis(y_train, trace_axes, last_t_axes).reshape( (-1,) + trace_shape) k_dd_cache = {} def get_k_train_train(get: Sequence[str]) -> _Kernel: if len(get) == 1: get = get[0] if get not in k_dd_cache: k_dd_cache[get] = kernel_fn(x_train, None, get, **kernel_fn_train_train_kwargs) elif len(get) == 2: if not any(g in k_dd_cache for g in get): k_dd_cache.update( kernel_fn(x_train, None, get, **kernel_fn_train_train_kwargs)._asdict()) else: for g in get: if g not in k_dd_cache: k_dd_cache[g] = kernel_fn(x_train, None, g, **kernel_fn_train_train_kwargs) else: raise ValueError(get) return _Kernel(**k_dd_cache) @lru_cache(2) def eigenspace(get: str): k_dd = getattr(get_k_train_train((get,)), get) k_dd = _add_diagonal_regularizer(utils.make_2d(k_dd), diag_reg, diag_reg_absolute_scale) evals, evecs = np.linalg.eigh(k_dd) evals = np.expand_dims(evals, 0) return evals, evecs @lru_cache(4) def predict_inf(get: Get): _, get = utils.canonicalize_get(get) k_dd = get_k_train_train(get) return gp_inference(k_dd, y_train, diag_reg, diag_reg_absolute_scale, trace_axes) def get_kernels(get: Get, x_test: Optional[np.ndarray], compute_cov: bool, **kernel_fn_test_test_kwargs): get = _get_dependency(get, compute_cov) k_dd = get_k_train_train(get) if x_test is None: k_td = None nngp_tt = compute_cov or None else: args_train, _ = utils.split_kwargs(kernel_fn_train_train_kwargs, x_train) args_test, _ = utils.split_kwargs(kernel_fn_test_test_kwargs, x_test) def is_array(x): return tree_all(tree_map( lambda x: isinstance(x, (onp.ndarray, np.ndarray)), x)) kwargs_td = dict(kernel_fn_train_train_kwargs) kwargs_tt = dict(kernel_fn_train_train_kwargs) for k in kernel_fn_test_test_kwargs: v_tt = kernel_fn_test_test_kwargs[k] v_dd = kernel_fn_train_train_kwargs[k] if is_array(v_dd) and is_array(v_tt): if (isinstance(v_dd, tuple) and len(v_dd) == 2 and isinstance(v_tt, tuple) and len(v_tt) == 2): v_td = (args_test[k], args_train[k]) else: v_td = v_tt elif v_dd != v_tt: raise ValueError(f'Same keyword argument {k} of `kernel_fn` is set to' f'different values {v_dd} != {v_tt} when computing ' f'the train-train and test-train/test-test kernels. ' f'If this is your intention, please submit a feature' f' request at ' f'https://github.com/google/neural-tangents/issues') else: v_td = v_tt kwargs_td[k] = v_td kwargs_tt[k] = v_tt k_td = kernel_fn(x_test, x_train, get, **kwargs_td) if compute_cov: nngp_tt = kernel_fn(x_test, None, 'nngp', **kwargs_tt) else: nngp_tt = None return k_dd, k_td, nngp_tt @utils.get_namedtuple('Gaussians') def predict_fn(t: Optional[ArrayOrScalar] = None, x_test: Optional[np.ndarray] = None, get: Optional[Get] = None, compute_cov: bool = False, **kernel_fn_test_test_kwargs) -> Dict[str, Gaussian]: """Return output mean and covariance on the test set at time[s] `t`. Args: t: a scalar of array of scalars of any shape. `t=None` is treated as infinity and returns the same result as `t=np.inf`, but is computed using linear solve for test predictions instead of eigendecomposition, saving time and precision. x_test: test inputs. `None` means to return non-regularized (`diag_reg=0`) predictions on the train-set inputs. For regularized predictions, pass `x_test=x_train`. get: string, the mode of the Gaussian process, either "nngp" or "ntk", or a tuple. `get=None` is equivalent to `get=("nngp", "ntk")`. compute_cov: if `True` computing both `mean` and `variance` and only `mean` otherwise. **kernel_fn_test_test_kwargs: optional keyword arguments passed to `kernel_fn`. See also `kernel_fn_train_train_kwargs` argument of the parent function. Returns: `fx_test_mean_t` or `(fx_test_mean_t, fx_test_cov_t)` if `compute_cov == True` with potentially additional leading time dimensions. """ if get is None: get = ('nngp', 'ntk') # train-train, test-train, test-test. k_dd, k_td, nngp_tt = get_kernels(get, x_test, compute_cov, **kernel_fn_test_test_kwargs) # Infinite time. if t is None: return predict_inf(get)(get=get, k_test_train=k_td, k_test_test=nngp_tt) # Finite time. t = np.array(t) * learning_rate t_shape = t.shape t = t.reshape((-1, 1)) def reshape_mean(mean): k = _get_first(k_dd if k_td is None else k_td) mean = mean.reshape(t_shape + k.shape[::2] + trace_shape) mean = np.moveaxis(mean, last_t_axes, trace_axes) return mean def reshape_cov(cov): k = _get_first(k_dd if k_td is None else k_td) cov_shape_t = t_shape + k.shape[::2] * 2 return utils.zip_axes(cov.reshape(cov_shape_t), len(t_shape)) out = {} for g in get: evals, evecs = eigenspace(g) # Training set. if k_td is None: mean = np.einsum( 'ji,ti,ki,k...->tj...', evecs, -expm1(evals, t), evecs, y_train_flat, optimize=_optimize()) # Test set. else: neg_inv_expm1 = -inv_expm1(evals, t) ktd_g = utils.make_2d(getattr(k_td, g)) mean = np.einsum( 'lj,ji,ti,ki,k...->tl...', ktd_g, evecs, neg_inv_expm1, evecs, y_train_flat, optimize=_optimize()) mean = reshape_mean(mean) if nngp_tt is not None: nngp_dd = utils.make_2d(k_dd.nngp) # Training set. if k_td is None: if g == 'nngp': cov = np.einsum( 'ji,ti,ki->tjk', evecs, (np.maximum(evals, 0.) * np.exp(- 2 * np.maximum(evals, 0.) * t / y_train.size)), evecs, optimize=_optimize()) elif g == 'ntk': exp = np.einsum( 'mi,ti,ki->tmk', evecs, np.exp(-np.maximum(evals, 0.) * t / y_train.size), evecs, optimize=_optimize()) cov = np.einsum( 'tmk,kl,tnl->tmn', exp, nngp_dd, exp, optimize=_optimize()) else: raise ValueError(g) # Test set. else: _nngp_tt = np.expand_dims(utils.make_2d(nngp_tt), 0) if g == 'nngp': cov = _nngp_tt - np.einsum( 'mj,ji,ti,ki,lk->tml', ktd_g, evecs, -inv_expm1(evals, 2 * t), evecs, ktd_g, optimize=_optimize()) elif g == 'ntk': term_1 = np.einsum( 'mi,ti,ki,lk->tml', evecs, neg_inv_expm1, evecs, ktd_g, optimize=_optimize()) term_2 = np.einsum( 'mj,ji,ti,ki,lk->tml', ktd_g, evecs, neg_inv_expm1, evecs, utils.make_2d(k_td.nngp), # pytype:disable=attribute-error optimize=_optimize()) term_2 += np.moveaxis(term_2, 1, 2) cov = np.einsum( 'tji,jk,tkl->til', term_1, nngp_dd, term_1, optimize=_optimize()) cov += -term_2 + _nngp_tt else: raise ValueError(g) out[g] = Gaussian(mean, reshape_cov(cov)) else: out[g] = mean return out return predict_fn
[docs]def max_learning_rate( ntk_train_train: np.ndarray, y_train_size: Optional[int] = None, momentum=0., eps: float = 1e-12) -> float: r"""Computes the maximal feasible learning rate for infinite width NNs. The network is assumed to be trained using mini-/full-batch GD + momentum with mean squared loss. The loss is assumed to have the form `1/(2 * batch_size * output_size) \|f(train_x) - train_y\|^2`. For vanilla SGD (i.e. `momentum = 0`) the maximal feasible learning rate is the largest `\eta` such that the operator `(I - \eta / (batch_size * output_size) * NTK)` is a contraction, which is `2 * batch_size * output_size * lambda_max(NTK)`. When `momentum > 0`, we use `2 * (1 + momentum) * batch_size * output_size * lambda_max(NTK)` (see `The Dynamics of Momentum` section in https://distill.pub/2017/momentum/). Args: ntk_train_train: analytic or empirical NTK on the training data. y_train_size: total training set output size, i.e. `f(x_train).size == y_train.size`. If `output_size=None` it is inferred from `ntk_train_train.shape` assuming `trace_axes=()`. momentum: The `momentum` for momentum optimizers. eps: a float to avoid zero divisor. Returns: The maximal feasible learning rate for infinite width NNs. """ ntk_train_train = utils.make_2d(ntk_train_train) factor = ntk_train_train.shape[0] if y_train_size is None else y_train_size if utils.is_on_cpu(ntk_train_train): max_eva = osp.linalg.eigvalsh(ntk_train_train, eigvals=(ntk_train_train.shape[0] - 1, ntk_train_train.shape[0] - 1))[-1] else: max_eva = np.linalg.eigvalsh(ntk_train_train)[-1] lr = 2 * (1 + momentum) * factor / (max_eva + eps) return lr
# INTERNAL UTILITIES def _optimize() -> str: """Return contraction order for `np.einsum` based on platform. Introduced after https://github.com/google/jax/pull/7512 since TPU seems to be more precise in `greeedy` mode. """ return 'greedy' if xla_bridge.get_backend().platform == 'tpu' else 'optimal' def _get_dependency(get: Get, compute_cov: bool) -> Tuple[str, ...]: """Figure out dependency for get.""" _, get = utils.canonicalize_get(get) for g in get: if g not in ['nngp', 'ntk']: raise NotImplementedError( 'Can only get either "nngp" or "ntk" predictions, got %s.' % g) get_dependency = () if 'nngp' in get or ('ntk' in get and compute_cov): get_dependency += ('nngp',) if 'ntk' in get: get_dependency += ('ntk',) return get_dependency def _get_fns_in_eigenbasis( k_train_train: np.ndarray, diag_reg: float, diag_reg_absolute_scale: bool, fns: Iterable[Callable[[np.ndarray, np.ndarray], np.ndarray]] ) -> Generator[Callable[[np.ndarray, np.ndarray], np.ndarray], None, None]: """Build functions of a matrix in its eigenbasis. Args: k_train_train: an n x n matrix. diag_reg: diagonal regularizer strength. diag_reg_absolute_scale: `True` to use absolute (vs relative to mean trace) regulatization. fns: a sequence of functions that add on the eigenvalues (evals, dt) -> modified_evals. Returns: A tuple of functions that act as functions of the matrix mat acting on vectors: `transform(vec, dt) = fn(mat, dt) @ vec` """ k_train_train = utils.make_2d(k_train_train) k_train_train = _add_diagonal_regularizer(k_train_train, diag_reg, diag_reg_absolute_scale) evals, evecs = np.linalg.eigh(k_train_train) evals = np.expand_dims(evals, 0) def to_eigenbasis(fn): """Generates a transform given a function on the eigenvalues.""" def new_fn(y_train, t): return np.einsum('ji,ti,ki,k...->tj...', evecs, fn(evals, t), evecs, y_train, optimize=_optimize()) return new_fn return (to_eigenbasis(fn) for fn in fns) def _add_diagonal_regularizer(A: np.ndarray, diag_reg: float, diag_reg_absolute_scale: bool) -> np.ndarray: dimension = A.shape[0] if not diag_reg_absolute_scale: diag_reg *= np.trace(A) / dimension return A + diag_reg * np.eye(dimension) def _get_cho_solve(A: np.ndarray, diag_reg: float, diag_reg_absolute_scale: bool, lower: bool = False) -> Callable[[np.ndarray, Axes], np.ndarray]: x_non_channel_shape = A.shape[1::2] A = utils.make_2d(A) A = _add_diagonal_regularizer(A, diag_reg, diag_reg_absolute_scale) C = sp.linalg.cho_factor(A, lower) def cho_solve(b: np.ndarray, b_axes: Axes) -> np.ndarray: b_axes = utils.canonicalize_axis(b_axes, b) last_b_axes = range(-len(b_axes), 0) x_shape = x_non_channel_shape + tuple(b.shape[a] for a in b_axes) b = np.moveaxis(b, b_axes, last_b_axes) b = b.reshape((A.shape[1], -1)) x = sp.linalg.cho_solve(C, b) x = x.reshape(x_shape) return x return cho_solve def _get_fx_test_shape(y_train: np.ndarray, k_test_train: np.ndarray, y_axes: Axes) -> Tuple[int, ...]: if k_test_train is None: return y_train.shape shape = list(k_test_train.shape[::2]) y_axes = utils.canonicalize_axis(y_axes, y_train) for i, c in enumerate(y_train.shape): if i in y_axes: shape.insert(i, c) return tuple(shape) def _make_expm1_fn(normalization: float): def expm1_fn(evals: np.ndarray, t: np.ndarray): # Since our matrix really should be positive semidefinite, # we can threshold the eigenvalues to squash ones that are negative # for numerical reasons. return np.expm1(-np.maximum(evals, 0.) * t / normalization) return expm1_fn def _make_inv_expm1_fn(normalization: float): expm1_fn = _make_expm1_fn(normalization) def _inv_expm1_fn(evals: np.ndarray, t: np.ndarray): return expm1_fn(evals, t) / np.abs(evals) return _inv_expm1_fn def _check_inputs(fx_train_or_state_0: Union[ArrayOrScalar, ODEState], fx_test_0: ArrayOrScalar, k_test_train: Optional[np.ndarray]): if isinstance(fx_train_or_state_0, ODEState): if fx_test_0 is not None: raise ValueError('`fx_test_0` is included in `ODEState` and must be set ' 'to `None`.') fx_train_0 = fx_train_or_state_0.fx_train fx_test_0 = fx_train_or_state_0.fx_test else: fx_train_0 = fx_train_or_state_0 if fx_train_0 is None and fx_test_0 is None: raise ValueError('Both `fx_train_0` and `fx_test_0` are `None`, i.e. no ' 'predictions will be computed.') if fx_test_0 is not None and k_test_train is None: raise ValueError('To get predictions on the test set, please provide ' '`k_test_train` kernel to the parent function.') def _get_axes(x: np.ndarray): n = x.ndim return ( tuple(range(0, n, 2)), tuple(range(1, n, 2)), tuple(range(0, n // 2)), tuple(range(n // 2, n)) ) def _get_first(k) -> np.ndarray: if isinstance(k, (onp.ndarray, np.ndarray)): return k for g in ('nngp', 'ntk'): if hasattr(k, g): v = getattr(k, g) if v is not None: return v raise ValueError(k) def _get_attr(k, g: str) -> np.ndarray: if isinstance(k, (onp.ndarray, np.ndarray)): return k return getattr(k, g)