Online Natural Gradient as a Kalman Filter

TL;DR

This paper demonstrates that natural gradient descent in statistical learning can be exactly interpreted as a Kalman filter, providing new insights into hyperparameters and extending to recurrent models.

Contribution

It establishes a rigorous equivalence between natural gradient methods and Kalman filtering, including for recurrent models, offering new theoretical insights.

Findings

Natural gradient corresponds to Kalman filtering in i.i.d. cases.

Joint Kalman filter over states and parameters aligns with natural gradient on RTRL.

Provides interpretations for natural gradient hyperparameters.

Abstract

We cast Amari's natural gradient in statistical learning as a specific case of Kalman filtering. Namely, applying an extended Kalman filter to estimate a fixed unknown parameter of a probabilistic model from a series of observations, is rigorously equivalent to estimating this parameter via an online stochastic natural gradient descent on the log-likelihood of the observations. In the i.i.d. case, this relation is a consequence of the "information filter" phrasing of the extended Kalman filter. In the recurrent (state space, non-i.i.d.) case, we prove that the joint Kalman filter over states and parameters is a natural gradient on top of real-time recurrent learning (RTRL), a classical algorithm to train recurrent models. This exact algebraic correspondence provides relevant interpretations for natural gradient hyperparameters such as learning rates or initialization and regularization of the Fisher information matrix.

Tables

Table 1. Table 1: Kalman filter objects vs natural gradient objects. The inputs are u t subscript 𝑢 𝑡 u_{t} , the predicted values are y ^ t subscript ^ 𝑦 𝑡 \hat{y}_{t} , and the model parameters are θ 𝜃 \theta .

iid (static, non-recurrent) model ${\hat{y}}_t=h(\theta ,u_t)$
Extended Kalman filter on static	Online natural gradient on $\theta$ with
parameter $\theta$	learning rate ${\eta }_t=1/\left(t+1\right)$
Covariance matrix $P_t$	Fisher information matrix $J_t={\eta }_t{P}_t^{-1}$
Bayesian prior $P_0$	Fisher matrix initialization $J_0=P_0^{-1}$
Fading memory	Larger or constant learning rate
Fading memory+constant prior	Fisher matrix regularization
Recurrent (state space) model ${\hat{y}}_t=\Phi ({\hat{y}}_{t-1},\theta ,u_t)$
Extended Kalman filter on $(\theta ,\hat{y})$	RTRL+natural gradient+state correction
Covariance of $\theta$ alone, $P^{\theta }$	Fisher matrix $J_t={\eta }_t{\left(P^{\theta }\right)}^{-1}$
Correlation between $\theta$ and ${\hat{y}}_t$	RTRL gradient estimate $\partial {\hat{y}}_t/\partial \theta$

Equations

$$y_t={\hat{y}}_t+\mathcal{N}(0,R_t),{\hat{y}}_t=h(\theta ,u_t)$$

$${\ell }_t=\frac{1}{2}({\hat{y}}_t-y_t{)}^{\top }{R}_t^{-1}({\hat{y}}_t-y_t)=-\ln p(y_t|{\hat{y}}_t)$$

$$p_{\beta }(y):=\frac{1}{Z(\beta )}{\mathrm{e}}^{{\sum }_k{\beta }_k{T}_k(y)}\lambda (dy)$$

$${\theta }_t\leftarrow {\theta }_{t-1}-{\eta }_t{\frac{\partial {\ell }_t(y_t)}{\partial \theta }}^{\top }$$

$$J(\theta ):={\mathbb{E}}_{y\sim p(y|\theta )}\left[{\frac{\partial \ln p(y|\theta )}{\partial \theta }}^{\otimes 2}\right]=-{\mathbb{E}}_{y\sim p(y|\theta )}\left[\frac{{\partial }^2\ln p(y|\theta )}{\partial {\theta }^2}\right]$$

$$\sum _t{\ell }_t(y_t),{\ell }_t(y):=-\ln p(y|{\hat{y}}_t)$$

$$J_t$$

$${\theta }_t$$

$$h_t$$

$${\theta }_t^{\ast }$$

$$J_0\leftarrow \mathrm{diag}(\text{fan-in})$$

$$L_t(\theta )=\ln p_{\theta }(y_t)+(1-{\lambda }_t)L_{t-1}(\theta ),L_0(\theta ):=-\frac{1}{2}(\theta -{\theta }_0{)}^{\top }{P}_0^{-1}(\theta -{\theta }_0)$$

$${\eta }_t={\gamma }_t,P_0={\eta }_0{J}_0^{-1},$$

$$1-{\lambda }_t={\eta }_{t-1}/{\eta }_t-{\eta }_{t-1}\text{for t⩾1}$$

$$L_t(\theta )=\ln p_{\theta }(y_t)+(1-{\lambda }_t)L_{t-1}(\theta ),L_0:=0$$

$${\theta }_t$$

$${\eta }_t={\gamma }_t,1-{\lambda }_t={\eta }_{t-1}/{\eta }_t-{\eta }_{t-1},{\eta }_0:={\eta }_1$$

$$s_t=f(s_{t-1},u_t)+\mathcal{N}(0,Q_t),{\hat{y}}_t=h(s_t,u_t),y_t\sim p(y|{\hat{y}}_t)$$

$$s_{t|t-1}$$

$$F_{t-1}$$

$$P_{t|t-1}$$

$${\hat{y}}_t$$

$$E_t$$

$$R_t$$

$$H_t$$

$$K_t$$

$$P_t$$

$$s_t$$

$$s_t=s_{t|t-1}-P_t{\left(\frac{\partial {\ell }_t(y_t)}{\partial s_{t|t-1}}\right)}^{\top }$$

$$E_t=R_t{\left(\frac{\partial \ln p(y_t|{\hat{y}}_t)}{\partial {\hat{y}}_t}\right)}^{\top }$$

$$\frac{\partial \ln p(y|\beta )}{\partial {\beta }_i}=T_i(y)-\mathbb{E}{T}_i=T_i(y)-{\hat{y}}_i$$

$$E_t=T(y_t)-{\hat{y}}_t={\left(\frac{\partial \ln p(y_t|\beta )}{\partial \beta }\right)}^{\top }$$

$$\frac{\partial \ln p(y_t|\beta )}{\partial \beta }=\frac{\partial \ln p(y_t|\hat{y})}{\partial \hat{y}}\frac{\partial \hat{y}}{\partial \beta }$$

$$P_t^{-1}\leftarrow P_{t|t-1}^{-1}+H_t^{\top }{R}_t^{-1}{H}_t$$

$$P_t^{-1}$$

$$s_t$$

Full text

In the i.i.d. case, this relation is a consequence of the “information filter” phrasing of the extended Kalman filter. In the recurrent (state space, non-i.i.d.) case, we prove that the joint Kalman filter over states and parameters is a natural gradient on top of real-time recurrent learning (RTRL), a classical algorithm to train recurrent models.

This exact algebraic correspondence provides relevant interpretations for natural gradient hyperparameters such as learning rates or initialization and regularization of the Fisher information matrix.

In statistical learning, stochastic gradient descent is a widely used tool to estimate the parameters of a model from empirical data, especially when the parameter dimension and the amount of data are large [BL03] (such as is typically the case with neural networks, for instance). The natural gradient [Ama98] is a tool from information geometry, which aims at correcting several shortcomings of the widely ordinary stochastic gradient descent, such as its sensitivity to rescalings or simple changes of variables in parameter space [Oll15]. The natural gradient modifies the ordinary gradient by using the information geometry of the statistical model, via the Fisher information matrix (see formal definition in Section 1.2; see also [Mar14]). The natural gradient comes with a theoretical guarantee of asymptotic optimality [Ama98] that the ordinary gradient lacks, and with the theoretical knowledge and various connections from information geometry, e.g., [AN00, OAAH17]. In large dimension, its computational complexity makes approximations necessary, e.g., [LMB07, Oll15, MCO16, GS15, MG15]; this has limited its adoption despite many desirable theoretical properties.

The extended Kalman filter (see e.g., the textbooks [Sim06, Sä13, Jaz70]) is a generic and effective tool to estimate in real time the state of a nonlinear dynamical system, from noisy measurements of some part or some function of the system. (The ordinary Kalman filter deals with linear systems.) Its use in navigation systems (GPS, vehicle control, spacecraft…), time series analysis, econometrics, etc. [Sä13], is extensive to the point it can been described as “one of the great discoveries of mathematical engineering” [GA15].

The goal of this text is to show that the natural gradient, when applied online, is a particular case of the extended Kalman filter. Indeed, the extended Kalman filter can be used to estimate the parameters of a statistical model (probability distribution), by viewing the parameters as the hidden state of a “static” dynamical system, and viewing i.i.d. samples as noisy observations depending on the parameters 111For this we slightly extend the definition of the Kalman filter to include discrete observations, by defining (Def. 5) the measurement error as $T(y)-\hat{y}$ instead of $y-\hat{y}$, where $T$ is the sufficient statistics of an exponential family model for output noise with mean $\hat{y}$. This reduces to the standard filter for Gaussian output noise, and naturally covers categorical outputs as often used in statistical learning (with $\hat{y}$ the class probabilities in a softmax classifier and $T$ a “one-hot” encoding of $y$).. We show that doing so is exactly equivalent to performing an online stochastic natural gradient descent (Theorem 2).

This results in a rigorous dictionary between the natural gradient objects from statistical learning, and the objects appearing in Kalman filtering; for instance, a larger learning rate for the natural gradient descent exactly corresponds to a fading memory in the Kalman filter (Proposition 3).

Table 1 lists a few correspondences between objects from the Kalman filter side and from the natural gradient side, as results from the theorems and propositions below. Note that the correspondence is one-sided: the online natural gradient is exactly an extended Kalman filter, but only corresponds to a particular use of the Kalman filter for parameter estimation problems (i.e., with static dynamics on the parameter part of the system).

Beyond the static case, we also consider the learning of the parameters of a general dynamical system, where subsequent observations exhibit temporal patterns instead of being i.i.d.; in statistical learning this is called a recurrent model, for instance, a recurrent neural network. We refer to [Jae02] for an introduction to recurrent models in statistical learning (recurrent neural networks) and the afferent techniques (including Kalman filters), and to [Hay01] for a clear, in-depth treatment of Kalman filtering for recurrent models. We prove (Theorem 12) that the extended Kalman filter applied jointly to the state and parameter, amounts to a natural gradient on top of real-time recurrent learning (RTRL), a classical (and costly) online algorithm for recurrent network training [Jae02].

Thus, we provide a bridge between techniques from large-scale statistical learning (natural gradient, RTRL) and a central object from mathematical engineering, signal processing, and estimation theory. Casting the natural gradient as a specific case of the extended Kalman filter is an instance of the provocative statement from [LS83] that “there is only one recursive identification method” that is optimal on quadratic functions. Indeed, the online natural gradient descent fits into the framework of [LS83, §3.4.5]. Arguably, this statement is limited to linear models, and for non-linear models one would expect different algorithms to coincide only at a certain order, or asymptotically; however, all the correspondences presented below are exact.

Related work.

In the i.i.d. (static) case, the natural gradient/Kalman filter correspondence follows from the information filter phrasing of Kalman filtering [Sim06, §6.2] by relatively direct manipulations. Nevertheless, we could find no reference in the literature explicitly identifying the two. [SW88] is an early example of the use of Kalman filtering for training feedforward neural networks in statistical learning, but does not mention the natural gradient. [RRK*+*92] argue that for neural networks, backpropagation, i.e., ordinary gradient descent, “is a degenerate form of the extended Kalman filter”. [Ber96] identifies the extended Kalman filter with a Gauss–Newton gradient descent for the specific case of nonlinear regression. [dFNG00] interprets process noise in the static Kalman filter as an adaptive, per-parameter learning rate, thus akin to a preconditioning matrix. [ŠKT01] uses the Fisher information matrix to study the variance of parameter estimation in Kalman-like filters, without using a natural gradient; [BL03] comment on the similarity between Kalman filtering and a version of Amari’s natural gradient for the specific case of least squares regression; [Mar14] and [Oll15] mention the relationship between natural gradient and the Gauss–Newton Hessian approximation; [Pat16] exploits the relationship between second-order gradient descent and Kalman filtering in specific cases including linear regression; [LCL*+*17] use a natural gradient descent over Gaussian distributions for an auxiliary problem arising in Kalman-like Bayesian filtering, a problem independent from the one treated here.

For the recurrent (non-i.i.d.) case, our result is that joint Kalman filtering is essentially a natural gradient on top of the classical RTRL algorithm for recurrent models [Jae02]. [Wil92] already observed that starting with the Kalman filter and introducing drastic simplifications (doing away with the covariance matrix) results in RTRL, while [Hay01, §5] contains statements that can be interpreted as relating Kalman filtering and preconditioned RTRL-like gradient descent for recurrent models (Section 3.2).

Perspectives.

In this text our goal is to derive the precise correspondence between natural gradient and Kalman filtering for parameter estimation (Thm. 2, Prop. 3, Prop. 4, Thm. 12), and to work out an exact dictionary between the mathematical objects on both sides. This correspondence suggests several possible venues for research, which nevertheless are not explored here.

First, the correspondence with the Kalman filter brings new interpretations and suggestions for several natural gradient hyperparameters, such as Fisher matrix initialization, equality between Fisher matrix decay rate and learning rate, or amount of regularization to the Fisher matrix (Section 2.2). The natural gradient can be quite sensitive to these hyperparameters. A first step would be to test the matrix decay rate and regularization values suggested by the Bayesian interpretation (Prop. 4) and see if they help with the natural gradient, or if these suggestions are overriden by the various approximations needed to apply the natural gradient in practice. These empirical tests are beyond the scope of the present study.

Next, since statistical learning deals with either continuous or categorical data, we had to extend the usual Kalman filter to such a setting. Traditionally, non-Gaussian output models have been treated by applying a nonlinearity to a standard Gaussian noise (Section 2.3). Instead, modeling the measurement noise as an exponential family (Appendix and Def. 5) allows for a unified treatment of the standard case (Gaussian output noise with known variance), of discrete categorical observations, or other exponential noise models (e.g., Gaussian noise with unknown variance). We did not test the empirical consequences of this choice, but it certainly makes the mathematical treatment flow smoothly, in particular the view of the Kalman filter as preconditioned gradient descent (Prop. 6).

Neither the natural gradient nor the extended Kalman filter scale well to large-dimensional models as currently used in machine learning, so that approximations are required. The correspondence raises the possibility that various methods developed for Kalman filtering (e.g., particle or unscented filters) or for natural gradient approximations (e.g., matrix factorizations such as the Kronecker product [MG15] or quasi-diagonal reductions [Oll15, MCO16]) could be transferred from one viewpoint to the other.

In statistical learning, other means have been developed to attain the same asymptotic efficiency as the natural gradient, notably trajectory averaging (e.g. [PJ92], or [Mar14] for the relationship to natural gradient) at little algorithmic cost. One may wonder if these can be generalized to filtering problems.

Proof techniques could be transferred as well: for instance, Amari [Ama98] gave a strong but sometimes informal argument that the natural gradient is Fisher-efficient, i.e., the resulting parameter estimate is asymptotically optimal for the Cramér–Rao bound; alternate proofs could be obtained by transferring related statements for the extended Kalman filter, e.g., combining techniques from [ŠKT01, BRD97, LS83].

Organization of the text.

In Section 1 we set the notation, recall the definition of the natural gradient (Def. 1), and explain how Kalman filtering can be used for parameter estimation in statistical learning (Section 1.3); the definition of the Kalman filter is included in Def. 5. Section 2 gives the main statements for viewing the natural gradient as an instance of an extended Kalman filter for i.i.d. observations (static systems), first intuitively via a heuristic asymptotic argument (Section 2.1), then rigorously (Thm. 2, Prop. 3, Prop. 4). The proof of these results appears in Section 2.3 and sheds some light on the geometry of Kalman filtering. Finally, the case of non-i.i.d. observations (recurrent or state space model) is treated in Section 3.

Acknowledgments.

Many thanks to Silvère Bonnabel, Gaétan Marceau-Caron, and the anonymous reviewers for their careful reading of the manuscript, corrections, and suggestions for the presentation and organization of the text. I would also like to thank Shun-ichi Amari, Frédéric Barbaresco, and Nando de Freitas for additional comments and for pointing out relevant references.

1 Problem setting, natural gradient, Kalman filter

1.1 Problem setting

In statistical learning, we have a series of observation pairs $(u_{1},y_{1}),\ldots,(u_{t},y_{t}),\ldots$ and want to predict $y_{t}$ from $u_{t}$ using a probabilistic model $p_{\theta}$. Assume for now that $y_{t}$ is real-valued (regression problem) and that the model for $y_{t}$ is a Gaussian centered on a predicted value $\hat{y}_{t}$, with known covariance matrix $R_{t}$, namely

[TABLE]

The function $h$ may represent any computation, for instance, a feedforward neural network with input $u$, parameters $\theta$, and output $\hat{y}$. The goal is to find the parameters $\theta$ such that the prediction $\hat{y}_{t}=h(\theta,u_{t})$ is as close as possible to $y_{t}$: the loss function is

[TABLE]

up to an additive constant.

For non-Gaussian outputs, we assume that the noise model on $y_{t}$ given $\hat{y}_{t}$ belongs to an exponential family, namely, that $\hat{y}_{t}$ is the mean parameter of an exponential family of distributions 222 The Appendix contains a reminder on exponential families. An exponential family of probability distributions on $y$, with sufficient statistics $T_{1}(y),\ldots,T_{K}(y)$, and with parameter $\beta\in{\mathbb{R}}^{K}$, is given by

$p_{\beta}(y)\mathrel{\mathop{:}}=\frac{1}{Z(\beta)}\,\mathrm{e}^{\sum_{k}\beta_{k}T_{k}(y)}\,\lambda(\operatorname{d}\!{}y)$

(1.3)

where $Z(\beta)$ is a normalizing constant, and $\lambda(\operatorname{d}\!{}y)$ is any reference measure on $y$. For instance, if $y\in{\mathbb{R}}^{K}$, $T_{k}(y)=y_{k}$ and $\lambda(\operatorname{d}\!{}y)$ is a Gaussian measure centered on [math], by varying $\beta$ one gets all Gaussian measures with the same covariance matrix and another mean. $y$ may be discrete, e.g., Bernoulli distributions correspond to $\lambda$ the uniform measure on $y\in\{0,1\}$ and a single sufficient statistic $T(0)=0$, $T(1)=1$. Often, the mean parameter $\bar{T}\mathrel{\mathop{:}}=\mathbb{E}_{y\sim p_{\beta}}T(y)$ is a more convenient parameterization than $\beta$. Exponential families maximize entropy (minimize information divergence from $\lambda$) for a given mean of $T$. over $y_{t}$; we again define the loss function as $\ell_{t}\mathrel{\mathop{:}}=-\ln p(y_{t}|\hat{y}_{t})$, and the output noise $R_{t}$ can be defined as the covariance matrix of the sufficient statistics of $y_{t}$ given this mean (Def. 5). For a Gaussian output noise this works as expected. For instance, for a classification problem, the output is categorical, $y_{t}\in\{1,\ldots,K\}$, and $\hat{y}_{t}$ will be the set of probabilities $\hat{y}_{t}=(p_{1},\ldots,p_{K-1})$ to have $y_{t}=1,\ldots,{K-1}$. In that case $R_{t}$ is the $(K-1)\times(K-1)$ matrix $(R_{t})_{kk^{\prime}}=\operatorname{diag}(p_{k})-p_{k}p_{k^{\prime}}$. (The last probability $p_{K}$ is determined by the others via $\sum p_{k}=1$ and has to be excluded to obtain a non-degenerate parameterization and an invertible covariance matrix $R_{t}$.)

This convention allows us to extend the definition of the Kalman filter to such a setting (Def. 5) in a natural way, just by replacing the measurement error $y_{t}-\hat{y}_{t}$ with $T(y_{t})-\hat{y}_{t}$, with $T$ the sufficient statistics for the exponential family. (For Gaussian noise this is the same, as $T(y)$ is $y$.)

In neural network terms, this means that the output layer of the network is fed to a loss function that is the log-loss of an exponential family, but places no restriction on the rest of the model.

General notation.

In statistical learning, the external inputs or regressor variables are often denoted $x$. In Kalman filtering, $x$ often denotes the state of the system, while the external inputs are often $u$. Thus we will avoid $x$ altogether and denote by $u$ the inputs and by $s$ the state of the system.

The variable to be predicted at time $t$ will be $y_{t}$, and $\hat{y}_{t}$ is the corresponding prediction. In general $\hat{y}_{t}$ and $y_{t}$ may be different objects in that $\hat{y}_{t}$ encodes a full probabilistic prediction for $y_{t}$. For Gaussians with known variance, $\hat{y}_{t}$ is just the predicted mean of $y_{t}$, so in this case $y_{t}$ and $\hat{y}_{t}$ are the same type of object. For Gaussians with unknown variance, $\hat{y}$ encodes both the mean and second moment of $y$. For discrete categorical data, $\hat{y}$ encodes the probability of each possible outcome $y$.

Thus, the formal setting for this text is as follows: we are given a sequence of finite-dimensional observations $(y_{t})$ with each $y_{t}\in{\mathbb{R}}^{\dim(y)}$, a sequence of inputs $(u_{t})$ with each $u_{t}\in{\mathbb{R}}^{\dim(u)}$, a parametric model $\hat{y}=h(\theta,u_{t})$ with parameter $\theta\in{\mathbb{R}}^{\dim(\theta)}$ and $h$ some fixed smooth function from ${\mathbb{R}}^{\dim(\theta)}\times{\mathbb{R}}^{\dim(u)}$ to ${\mathbb{R}}^{\dim(\hat{y})}$. We are given an exponential family (output noise model) $p(y|\hat{y})$ on $y$ with mean parameter $\hat{y}$ and sufficient statistics $T(y)$ (see the Appendix), and we define the loss function $\ell_{t}\mathrel{\mathop{:}}=-\ln p(y_{t}|\hat{y}_{t})$.

The natural gradient descent on parameter $\theta_{t}$ will use the Fisher matrix $J_{t}$. The Kalman filter will have posterior covariance matrix $P_{t}$.

For multidimensional quantities $x$ and $y=f(x)$, we denote by $\frac{\partial y}{\partial x}$ the Jacobian matrix of $y$ w.r.t. $x$, whose $(i,j)$ entry is $\frac{\partial f_{i}(x)}{\partial x_{j}}$. This satisfies the chain rule $\frac{\partial z}{\partial y}\frac{\partial y}{\partial x}=\frac{\partial z}{\partial x}$. With this convention, gradients of real-valued functions are row vectors, so that a gradient descent takes the form $x\leftarrow x-\eta\,(\partial f/\partial x)^{\!\top}\!$.

For a column vector $u$, $u^{\otimes 2}$ is synonymous with $uu^{\!\top}\!$, and with $u^{\!\top}\!u$ for a row vector.

1.2 Natural gradient descent

A standard approach to optimize the parameter $\theta$ of a probabilistic model, given a sequence of observations $(y_{t})$, is an online gradient descent

[TABLE]

with learning rate $\eta_{t}$. This simple gradient descent is particularly suitable for large datasets and large-dimensional models [BL03], but has several practical and theoretical shortcomings. For instance, it uses the same non-adaptive learning rate for all parameter components. Moreover, simple changes in parameter encoding or in data presentation (e.g., encoding black and white in images by 0/1 or 1/0) can result in different learning performance.

This motivated the introduction of the natural gradient [Ama98]. It is built to achieve invariance with respect to parameter re-encoding; in particular, learning become insensitive to the characteristic scale of each parameter direction, so that different directions naturally get suitable learning rates. The natural gradient is the only general way to achieve such invariance [AN00, §2.4].

The natural gradient preconditions the gradient descent with $J(\theta)^{-1}$ where $J$ is the Fisher information matrix [Kul97] with respect to the parameter $\theta$. For a smooth probabilistic model $p(y|\theta)$ over a random variable $y$ with parameter $\theta$, the latter is defined as

[TABLE]

Definition 1 below formally introduces the online natural gradient. If the model for $y$ involves an input $u$, then an expectation or empirical average over the input is introduced in the definition of $J$ [AN00, §8.2] [Mar14, §5].

However, this comes at a large computational cost for large-dimensional models: just storing the Fisher matrix already costs $O((\dim\theta)^{2})$. Various strategies are available to approximate the natural gradient for complex models such as neural networks, using diagonal or block-diagonal approximation schemes for the Fisher matrix, e.g., [LMB07, Oll15, MCO16, GS15, MG15].

Definition 1 (Online natural gradient).

Consider a statistical model with parameter $\theta$ that predicts an output $y$ given an input $u$. Suppose that the prediction takes the form $y\sim p(y|\hat{y})$ where $\hat{y}=h(\theta,u)$ depends on the input via a model $h$ with parameter $\theta$. Given observation pairs $(u_{t},y_{t})$, the goal is to minimize, online, the loss function

[TABLE]

as a function of $\theta$.

The online natural gradient maintains a current estimate $\theta_{t}$ of the parameter $\theta$, and a current approximation $J_{t}$ of the Fisher matrix. The parameter is estimated by a gradient descent with preconditioning matrix $J_{t}^{-1}$, namely

[TABLE]

with learning rate $\eta_{t}$ and Fisher matrix decay rate $\gamma_{t}$.

In the Fisher matrix update, the expectation over all possible values $y\sim p(y|\hat{y})$ can often be computed algebraically, but this is sometimes computationally bothersome (for instance, in neural networks, it requires $\dim(\hat{y}_{t})$ distinct backpropagation steps [Oll15]). A common solution [APF00, LMB07, Oll15, PB13] is to just use the value $y=y_{t}$ (outer product approximation) instead of the expectation over $y$. Another is to use a Monte Carlo approximation with a single sample of $y\sim p(y|\hat{y}_{t})$ [Oll15, MCO16], namely, using the gradient of a synthetic sample instead of the actual observation $y_{t}$ in the Fisher matrix. These latter two solutions are often confused; only the latter provides an unbiased estimate, see discussion in [Oll15, PB13].

The online “smoothed” update of the Fisher matrix in (1.7) mixes past and present estimates (this or similar updates are used in [LMB07, MCO16]). The reason is at least twofold. First, the “genuine” Fisher matrix involves an expectation over the inputs $u_{t}$ [AN00, §8.2]: this can be approximated online only via a moving average over inputs (e.g., $\gamma_{t}=1/t$ realizes an equal-weight average over all inputs seen so far). Second, the expectation over $y\sim p(y|\hat{y}_{t})$ in (1.7) is often replaced with a Monte Carlo estimation with only one value of $y$, and averaging over time compensates for this Monte Carlo sampling.

As a consequence, since $\theta_{t}$ changes over time, this means that the estimate $J_{t}$ mixes values obtained at different values of $\theta$, and converges to the Fisher matrix only if $\theta_{t}$ changes slowly, i.e., if $\eta_{t}\to 0$. The correspondence below with Kalman filtering suggests using $\gamma_{t}=\eta_{t}$.

1.3 Kalman filtering for parameter estimation

One possible definition of the extended Kalman filter is as follows [Sim06, §15.1]. We are trying to estimate the current state of a dynamical system $s_{t}$ whose evolution equation is known but whose precise value is unknown; at each time step, we have access to a noisy measurement $y_{t}$ of a quantity $\hat{y}_{t}=h(s_{t})$ which depends on this state.

The Kalman filter maintains an approximation of a Bayesian posterior on $s_{t}$ given the observations $y_{1},\ldots,y_{t}$. The posterior distribution after $t$ observations is approximated by a Gaussian with mean $s_{t}$ and covariance matrix $P_{t}$. (Indeed, Bayesian posteriors always tend to Gaussians asymptotically under mild conditions, by the Bernstein–von Mises theorem [vdV00].) The Kalman filter prescribes a way to update $s_{t}$ and $P_{t}$ when new observations become available.

The Kalman filter update is summarized in Definition 5 below. It is built to provide the exact value of the Bayesian posterior in the case of linear dynamical systems with Gaussian measurements and a Gaussian prior. In that sense, it is exact at first order.

The Kalman filtering viewpoint on a statistical learning problem is that we are facing a system with hidden variable $\theta$, with an unknown value that does not evolve in time, and that the observations $y_{t}$ bring more and more information on $\theta$. Thus, a statistical learning problem can be tackled by applying the extended Kalman filter to the unknown variable $s_{t}=\theta$, whose underlying dynamics from time $t$ to time $t+1$ is just to remain unchanged ($f=\operatorname{Id}$ and noise on $s$ is [math] in Definition 5). In such a setting, the posterior covariance matrix $P_{t}$ will generally tend to [math] as observations accumulate and the parameter is identified better333But $P_{t}$ must still be maintained even if it tends to [math], since it is used to update the parameter at the correct rate. (this occurs at rate $1/t$ for the basic filter, which estimates from all $t$ past observations at time $t$, or at other rates if fading memory is included, see below). The initialization $\theta_{0}$ and its covariance $P_{0}$ can be interpreted as Bayesian priors on $\theta$ [SW88, LS83].

We will refer to this as a static Kalman filter. In the static case and without fading memory, the posterior covariance $P_{t}$ after $t$ observations will decrease like $O(1/t)$, so that the parameter gets updated by $O(1/t)$ after each new observation. Introducing fading memory for past observations (equivalent to adding noise on $\theta$ at each step, $Q_{t}\propto P_{t|t-1}$ in Def. 5) leads to a larger covariance and faster updates.

An example: Feedforward neural networks.

The Kalman approach above can be applied to any parametric statistical model. For instance [SW88] treat the case of a feedforward neural network. In our setting this is described as follows. Let $u$ be the input of the model and $y$ the true (desired) output. A feedforward neural network can be described as a function $\hat{y}=h(\theta,u)$ where $\theta$ is the set of all parameters of the network, where $h$ represents all computations performed by the network on input $u$, and $\hat{y}$ encodes the network prediction for the value of the output $y$ on input $u$. For categorical observations $y$, $\hat{y}$ is usually a set of predicted probabilities for all possible classes; while for regression problems, $\hat{y}$ is directly the predicted value. In both cases, the error function to be minimized can be defined as $\ell(y)\mathrel{\mathop{:}}=-\ln p(y|\hat{y})$: in the regression case, $\hat{y}$ is interpreted as a mean of a Gaussian model on $y$, so that $-\ln p(y|\hat{y})$ is the square error up to a constant.

Training the neural network amounts to estimating the network parameter $\theta$ from the observations. Applying a static Kalman filter for this problem [SW88] amounts to using Def. 5 with $s=\theta$, $f=\operatorname{Id}$ and $Q=0$. At first glance this looks quite different from the common gradient descent (backpropagation) approach for neural networks. The backpropagation operation is represented in the Kalman filter by the computation of $H=\frac{\partial h(s,u)}{\partial s}$ (2.17) where $s$ is the parameter. We show that the additional operations of the Kalman filter correspond to using a natural gradient instead of a vanilla gradient.

Unfortunately, for models with high-dimensional parameters such as neural networks, the Kalman filter is computationally costly and requires block-diagonal approximations for $P_{t}$ (which is a square matrix of size $\dim\theta$); moreover, computing $H_{t}=\partial\hat{y}_{t}/\partial\theta$ is needed in the filter, and requires doing one separate backpropagation for each component of the output $\hat{y}_{t}$.

2 Natural gradient as a Kalman filter: the static (i.i.d.) case

We now write the explicit correspondence between an online natural gradient to estimate the parameter of a statistical model from i.i.d. observations, and a static extended Kalman filter. We first give a heuristic argument that outlines the main ideas from the proof (Section 2.1).

Then we state the formal correspondences. First, the static Kalman filter corresponds to an online natural gradient with learning rate $1/t$ (Thm. 2). The rate $1/t$ arises because such a filter takes into account all previous evidence without decay factors (and with process noise $Q=0$ in the Kalman filter), thus the posterior covariance matrix decreases like $O(1/t)$. Asymptotically, this is the optimal rate in statistical learning [Ama98]. (Note, however, that the online natural gradient and extended Kalman filter are identical at every time step, not only asymptotically.)

The $1/t$ rate is often too slow in practical applications, especially when starting far away from an optimal parameter value. The natural gradient/Kalman filter correspondence is not specific to the $O(1/t)$ rate. Larger learning rates in the natural gradient correspond to a fading memory Kalman filter (adding process noise $Q$ proportional to the posterior covariance at each step, corresponding to a decay factor for the weight of previous observations); this is Proposition 3. In such a setting, the posterior covariance matrix in the Kalman filter does not decrease like $O(1/t)$; for instance, a fixed decay factor for the fading memory corresponds to a constant learning rate.

Finally, a fading memory in the Kalman filter may erase prior Bayesian information $(\theta_{0},P_{0})$ too fast; maintaining the weight of the prior in a fading memory Kalman filter is treated in Proposition 4 and corresponds, on the natural gradient side, to a so-called weight decay [Bis06] towards $\theta_{0}$ together with a regularization of the Fisher matrix, at specific rates.

2.1 Natural gradient as a Kalman filter: heuristics

As a first ingredient in the correspondence, we interpret Kalman filters as gradient descents: the extended Kalman filter actually performs a gradient descent on the log-likelihood of each new observation, with preconditioning matrix equal to the posterior covariance matrix. This is Proposition 6 below. This relies on having an exponential family as the output noise model.

Meanwhile, the natural gradient uses the Fisher matrix as a preconditioning matrix. The Fisher matrix is the average Hessian of log-likelihood, thanks to the classical double definition of the Fisher matrix as square gradient or Hessian, $J(\theta)=\mathbb{E}_{y\sim p(y|\theta)}\left[\frac{\partial\ln p(y)}{\partial\theta}^{\otimes 2}\right]=-\mathbb{E}_{y\sim p(y|\theta)}\left[\frac{\partial^{2}\ln p(y)}{\partial\theta^{2}}\right]$ for any probabilistic model $p(y|\theta)$ [Kul97].

Assume that the probability of the data given the parameter $\theta$ is approximately Gaussian, $p(y_{1},\ldots,y_{t}|\theta)\propto\exp(-(\theta-\theta^{\ast})^{\!\top}\!\Sigma^{-1}(\theta-\theta^{\ast}))$ with covariance $\Sigma$. This often holds asymptotically thanks to the Bernstein–von Mises theorem; moreover, the posterior covariance $\Sigma$ typically decreases like $1/t$. Then the Hessian (w.r.t. $\theta$) of the total log-likelihood of $(y_{1},\ldots,y_{t})$ is $\Sigma^{-1}$, the inverse covariance of $\theta$. So the average Hessian per data point, the Fisher matrix $J$, is approximately $J\approx\Sigma^{-1}/t$. Since a Kalman filter to estimate $\theta$ is essentially a gradient descent preconditioned with $\Sigma$, it will be the same as using a natural gradient with learning rate $1/t$. Using a fading memory Kalman filter will estimate $\Sigma$ from fewer past observations and provide larger learning rates.

Another way to understand the link between natural gradient and Kalman filter is as a second-order Taylor expansion of data log-likelihood. Assume that the total data log-likelihood at time $t$, $L_{t}(\theta)\mathrel{\mathop{:}}=-\sum_{s=1}^{t}\ln p(y_{s}|\theta)$, is approximately quadratic as a function of $\theta$, with a minimum at $\theta^{\ast}_{t}$ and a Hessian $h_{t}$, namely, $L_{t}(\theta)\approx\frac{1}{2}(\theta-\theta^{\ast}_{t})^{\!\top}\!h_{t}(\theta-\theta^{\ast}_{t})$. Then when new data points become available, this quadratic approximation would be updated as follows (online Newton method):

[TABLE]

and indeed these are equalities for a quadratic log-likelihood. Namely, the update of $\theta^{\ast}_{t}$ is a gradient ascent on log-likelihood, preconditioned by the inverse Hessian (Newton method). Note that $h_{t}$ grows like $t$ (each data point adds its own contribution). Thus, $h_{t}$ is $t$ times the empirical average of the Hessian, i.e., approximately $t$ times the Fisher matrix of the model ($h_{t}\approx tJ$). So this update is approximately a natural gradient descent with learning rate $1/t$.

Meanwhile, the Bayesian posterior on $\theta$ (with uniform prior) after observations $y_{1},\ldots,y_{t}$ is proportional to $e^{-L_{t}}$ by definition of $L_{t}$. If $L_{t}\approx\frac{1}{2}(\theta-\theta^{\ast}_{t})^{\!\top}\!h_{t}(\theta-\theta^{\ast}_{t})$, this is a Gaussian distribution centered at $\theta^{\ast}_{t}$ with covariance matrix $h_{t}^{-1}$. The Kalman filter is built to maintain an approximation $P_{t}$ of this covariance matrix $h_{t}^{-1}$, and then performs a gradient step preconditioned on $P_{t}$ similar to (2.2).

The simplest situation corresponds to an asymptotic rate $O(1/t)$, i.e., estimating the parameter based on all past evidence; the update (2.1) of the Hessian is additive, so that $h_{t}$ grows like $t$ and $h_{t}^{-1}$ in (2.2) produces an effective learning rate $O(1/t)$. Introducing a decay factor for older observations, multiplying the term $h_{t-1}$ in (2.1), produces a fading memory effect and results in larger learning rates.

These heuristics justify the statement from [LS83] that “there is only one recursive identification method”. Close to an optimum (so that the Hessian is positive), all second-order algorithms are essentially an online Newton step (2.1)-(2.2) approximated in various ways.

But even though this heuristic argument appears to be approximate or asymptotic, the correspondence between online natural gradient and Kalman filter presented below is exact at every time step.

2.2 Statement of the correspondence, static (i.i.d.) case

For the statement of the correspondence, we assume that the output noise on $y$ given $\hat{y}$ is modelled by an exponential family with mean parameter $\hat{y}$. This covers the traditional Gaussian case $y=\mathcal{N}(\hat{y},\Sigma)$ with fixed $\Sigma$ often used in Kalman filters. The Appendix contains necessary background on exponential families.

Theorem 2 (Natural gradient as a static Kalman filter).

These two algorithms are identical under the correspondence $(\theta_{t},J_{t})\leftrightarrow(s_{t},P_{t}^{-1}/(t+1))$:

The online natural gradient (Def. 1) with learning rates $\eta_{t}=\gamma_{t}=1/(t+1)$, applied to learn the parameter $\theta$ of a model that predicts observations $(y_{t})$ with inputs $(u_{t})$, using a probabilistic model $y\sim p(y|\hat{y})$ with $\hat{y}=h(\theta,u)$, where $h$ is any model and $p(y|\hat{y})$ is an exponential family with mean parameter $\hat{y}$. 2. 2.

The extended Kalman filter (Def. 5) to estimate the state $s$ from observations $(y_{t})$ and inputs $(u_{t})$, using a probabilistic model $y\sim p(y|\hat{y})$ with $\hat{y}=h(s,u)$ and $p(y|\hat{y})$ an exponential family with mean parameter $\hat{y}$, with static dynamics and no added noise on $s$ ($f(s,u)=s$ and $Q=0$ in Def. 5).

Namely, if at startup $(\theta_{0},J_{0})=(s_{0},P_{0}^{-1})$, then $(\theta_{t},J_{t})=(s_{t},P_{t}^{-1}/(t+1))$ for all $t\geqslant 0$.

The correspondence is exact only if the Fisher metric is updated before the parameter in the natural gradient descent (as in Definition 1).

The correspondence with a Kalman filter provides an interpretation for various hyper-parameters of online natural gradient descent. In particular, $J_{0}=P_{0}^{-1}$ can be interpreted as the inverse covariance of a Bayesian prior on $\theta$ [SW88]. This relates the initialization $J_{0}$ of the Fisher matrix to the initialization of $\theta$: for instance, in neural networks it is recommended to initialize the weights according to a Gaussian of covariance $\operatorname{diag}(1/\text{fan-in})$ (number of incoming weights) for each neuron; interpreting this as a Bayesian prior on weights, one may recommend to initialize the Fisher matrix to the inverse of this covariance, namely,

[TABLE]

Indeed this seemed to perform quite well in small-scale experiments.

Learning rates, fading memory, and metric decay rate.

Theorem 2 exhibits a $1/(t+1)$ learning rate for the online natural gradient. This is because the static Kalman filter for i.i.d. observations approximates the maximum a posteriori (MAP) of the parameter $\theta$ based on all past observations; MAP and maximum likelihood estimators change by $O(1/t)$ when a new data point is observed.

However, for nonlinear systems, optimality of the $1/t$ rate only occurs asymptotically, close enough to the optimum. In general, a $1/(t+1)$ learning rate is far from optimal if optimization does not start close to the optimum or if one is not using the exact Fisher matrix $J_{t}$ or covariance matrix $P_{t}$.

Larger effective learning rates are achieved thanks to so-called “fading memory” variants of the Kalman filter, which put less weight on older observations. For instance, one may multiply the log-likelihood of previous points by a forgetting factor $(1-\lambda_{t})$ before each new observation. This is equivalent to an additional step $P_{t-1}\leftarrow P_{t-1}/(1-\lambda_{t})$ in the Kalman filter, or to the addition of an artificial process noise $Q_{t}$ proportional to $P_{t-1}$ in the model. Such strategies are reported to often improve performance, especially when the data do not truly follow the model [Sim06, §5.5, §7.4], [Hay01, §5.2.2]. See for instance [Ber96] for the relationship between Kalman fading memory and gradient descent learning rates (in a particular case).

Proposition 3 (Natural gradient rates and fading memory).

Under the same model and assumptions as in Theorem 2, the following two algorithms are identical via the correspondence $(\theta_{t},J_{t})\leftrightarrow(s_{t},\eta_{t}P_{t}^{-1})$:

• •

An online natural gradient step with learning rate $\eta_{t}$ and metric decay rate $\gamma_{t}$

• •

A fading memory Kalman filter with an additional step $P_{t-1}\leftarrow P_{t-1}/(1-\lambda_{t})$ before the transition step; such a filter iteratively optimizes a weighted log-likelihood function $L_{t}$ of recent observations, with decay $(1-\lambda_{t})$ at each step, namely:

[TABLE]

provided the following relations are satisfied:

[TABLE]

For example, taking $\eta_{t}=1/(t+\mathrm{cst})$ corresponds to $\lambda_{t}=0$, no decay for older observations, and an initial covariance $P_{0}=J_{0}^{-1}/\mathrm{cst}$. Taking a constant learning rate $\eta_{t}=\eta_{0}$ corresponds to a constant decay factor $\lambda=\eta_{0}$.

The proposition above computes the fading memory decay factors $1-\lambda_{t}$ from the natural gradient learning rates $\eta_{t}$ via (2.6). In the other direction, one can start with the decay factors $\lambda_{t}$ and obtain the learning rates $\eta_{t}$ via the cumulated sum of weights $S_{t}$: $S_{0}\mathrel{\mathop{:}}=1/\eta_{0}$ then $S_{t}\mathrel{\mathop{:}}=(1-\lambda_{t})S_{t-1}+1$, then $\eta_{t}\mathrel{\mathop{:}}=1/S_{t}$. This clarifies how $\lambda_{t}=0$ corresponds to $\eta_{t}=1/(t+\mathrm{cst})$ where the constant is $S_{0}$.

The learning rates also control the weight given to the Bayesian prior and to the starting point $\theta_{0}$. For instance, with $\eta_{t}=1/(t+t_{0})$ and large $t_{0}$, the gradient descent will move away slowly from $\theta_{0}$; in the Kalman interpretation this corresponds to $\lambda_{t}=0$ and a small initial covariance $P_{0}=J_{0}^{-1}/t_{0}$ around $\theta_{0}$, so that the prior weighs as much as $t_{0}$ observations.

This result suggests to set $\gamma_{t}=\eta_{t}$ in the online natural gradient descent of Definition 1. The intuitive explanation for this setting is as follows: Both the Kalman filter and the natural gradient build a second-order approximation of the log-likelihood of past observations as a function of the parameter $\theta$, as explained in Section 2.1. Using a fading memory corresponds to putting smaller weights on past observations; these weights affect the first-order and the second-order parts of the approximation in the same way. In the gradient viewpoint, the learning rate $\eta_{t}$ corresponds to the first-order term (comparing (1.8) and (2.2)) while the Fisher matrix decay rate corresponds to the rate at which the second-order information is updated. Thus, the setting $\eta_{t}=\gamma_{t}$ in the natural gradient corresponds to using the same decay weights for the first-order and second-order expansion of the log-likelihood of past observations.

Still, one should keep in mind that the extended Kalman filter is itself only an approximation for nonlinear systems. Moreover, from a statistical point of view, the second-order object $J_{t}$ is higher-dimensional than the first-order information, so that estimating $J_{t}$ based on more past observations may be more stable. Finally, for large-dimensional problems the Fisher matrix is always approximated, which affects optimality of the learning rates. So in practice, considering $\gamma_{t}$ and $\eta_{t}$ as hyperparameters to be tuned independently may still be beneficial, though $\gamma_{t}=\eta_{t}$ seems a good place to start.

Regularization of the Fisher matrix and Bayesian priors.

A potential downside of fading memory in the Kalman filter is that the Bayesian interpretation is partially lost, because the Bayesian prior is forgotten too quickly. For instance, with a constant learning rate, the weight of the Bayesian prior decreases exponentially; likewise, with $\eta_{t}=O(1/\sqrt{t})$, the filter essentially works with the $O(\sqrt{t})$ most recent observations, while the weight of the prior decreases like $\approx e^{-\sqrt{t}}$ (as does the weight of the earliest observations; this is the product $\prod(1-\lambda_{t})$). But precisely, when working with fewer data points one may wish the prior to play a greater role.

The Bayesian interpretation can be restored by explicitly optimizing a combination of the log-likelihood of recent points, and the log-likelihood of the prior. This is implemented in Proposition 4.

From the natural gradient viewpoint, this translates both as a regularization of the Fisher matrix (often useful in practice to numerically stabilize its inversion) and of the gradient step. With a Gaussian prior $\mathcal{N}(\theta_{\mathrm{prior}},\operatorname{Id})$, this manifests as an additional step towards $\theta_{\mathrm{prior}}$ and adding $\varepsilon.\operatorname{Id}$ to the Fisher matrix, known respectively as weight decay and Tikhonov regularization [Bis06, §3.3, §5.5] in statistical learning.

Proposition 4 (Bayesian regularization of the Fisher matrix).

Let $\pi=\mathcal{N}(\theta_{\mathrm{prior}},\Sigma_{0})$ be a Gaussian prior on $\theta$. Under the same model and assumptions as in Theorem 2, the following two algorithms are equivalent:

• •

A modified fading memory Kalman filter that iteratively optimizes $L_{t}(\theta)+n_{\mathrm{prior}}\ln\pi(\theta)$ where $L_{t}$ is a weighted log-likelihood function of recent observations with decay $(1-\lambda_{t})$:

[TABLE]

initialized with $P_{0}=\frac{\eta_{1}}{1+n_{\mathrm{prior}}\eta_{1}}\Sigma_{0}$.

• •

A regularized online natural gradient step with learning rate $\eta_{t}$ and metric decay rate $\gamma_{t}$, initialized with $J_{0}=\Sigma_{0}^{-1}$,

[TABLE]

provided the following relations are satisfied:

[TABLE]

Thus, the regularization terms are fully determined by choosing the learning rates $\eta_{t}$, a prior such as $\mathcal{N}(0,1/\text{fan-in})$ (for neural networks), and a value of $n_{\mathrm{prior}}$ such as $n_{\mathrm{prior}}=1$ (the prior weighs as much as $n_{\mathrm{prior}}$ data points). This holds both for regularization of the Fisher matrix $J_{t}+\eta_{t}n_{\mathrm{prior}}\Sigma_{0}^{-1}$, and for regularization of the parameter via the extra gradient step $\lambda_{t}n_{\mathrm{prior}}\Sigma_{0}^{-1}(\theta-\theta_{\text{prior}})$.

The relative strength of regularization in the Fisher matrix decreases like $\eta_{t}$. In particular, a constant learning rate results in a constant regularization.

The added gradient step $\lambda_{t}n_{\mathrm{prior}}\Sigma_{0}^{-1}(\theta-\theta_{\text{prior}})$ is modulated by $\lambda_{t}$ which depends on $\eta_{t}$; this extra term pulls towards the prior $\theta_{\mathrm{prior}}$. The Bayesian viewpoint guarantees that this extra term will not ultimately prevent convergence of the gradient descent (as the influence of the prior vanishes when the number of observations increases).

It is not clear how much these recommendations for natural gradient descent coming from its Bayesian interpretation are sensitive to using only an approximation of the Fisher matrix.

2.3 Proofs for the static case

The proof of Theorem 2 starts with the interpretation of the Kalman filter as a gradient descent (Proposition 6).

We first recall the exact definition and the notation we use for the extended Kalman filter.

Definition 5 (Extended Kalman filter).

Consider a dynamical system with state $s_{t}$, inputs $u_{t}$ and outputs $y_{t}$,

[TABLE]

where $p(\cdot|\hat{y})$ denotes an exponential family with mean parameter $\hat{y}$ (e.g., $y=\mathcal{N}(\hat{y},R)$ with fixed covariance matrix $R$).

The extended Kalman filter for this dynamical system estimates the current state $s_{t}$ given observations $y_{1},\ldots,y_{t}$ in a Bayesian fashion. At each time, the Bayesian posterior distribution of the state given $y_{1},\ldots,y_{t}$ is approximated by a Gaussian $\mathcal{N}(s_{t},P_{t})$ so that $s_{t}$ is the approximate maximum a posteriori, and $P_{t}$ is the approximate posterior covariance matrix. (The prior is $\mathcal{N}(s_{0},P_{0})$ at time [math].) Each time a new observation $y_{t}$ is available, these estimates are updated as follows.

The transition step (before observing $y_{t}$) is

[TABLE]

and the observation step after observing $y_{t}$ is

[TABLE]

For non-Gaussian output noise, the definition of $E_{t}$ and $R_{t}$ above via the mean parameter $\hat{y}$ of an exponential family, differs from the practice of modelling non-Gaussian noise via a nonlinear function applied to Gaussian noise. This allows for a straightforward treatment of various output models, such as discrete outputs or Gaussians with unknown variance. In the Gaussian case with known variance our definition is fully standard. 444Non-Gaussian output noise is often modelled in Kalman filtering via a continuous nonlinear function applied to a Gaussian noise [Sim06, 13.1]; this cannot easily represent discrete random variables. Moreover, since the filter linearizes the function around the [math] value of the noise [Sim06, 13.1], the noise is still implicitly Gaussian, though with a state-dependent variance.

The proof starts with the interpretation of the Kalman filter as a gradient descent preconditioned by $P_{t}$. Compare this result and Lemma 9 to [Hay01, (5.68)–(5.73)].

Proposition 6 (Kalman filter as preconditioned gradient descent).

The update of the state $s$ in a Kalman filter can be seen as an online gradient descent on data log-likelihood, with preconditioning matrix $P_{t}$. More precisely, denoting $\ell_{t}(y)\mathrel{\mathop{:}}=-\ln p(y|\hat{y}_{t})$, the update (2.20) is equivalent to

[TABLE]

where in the derivative, $\ell_{t}$ depends on $s_{t|t-1}$ via $\hat{y}_{t}=h(s_{t|t-1},u_{t})$.

Lemma 7 (Errors and gradients).

When the output model is an exponential family with mean parameter $\hat{y}_{t}$, the error $E_{t}$ is related to the gradient of the log-likelihood of the observation $y_{t}$ with respect to the prediction $\hat{y}_{t}$ by

[TABLE]

Proof of the lemma.

.

For a Gaussian $y_{t}=\mathcal{N}(\hat{y}_{t},R)$, this is just a direct computation. For a general exponential family, consider the natural parameter $\beta$ of the exponential family which defines the law of $y$, namely, $p(y|\beta)=\exp(\sum_{i}\beta_{i}T_{i}(y))/Z(\beta)$ with sufficient statistics $T_{i}$ and normalizing constant $Z$. An elementary computation (Appendix, (A.3)) shows that

[TABLE]

by definition of the mean parameter $\hat{y}$. Thus,

[TABLE]

where the derivative is with respect to the natural parameter $\beta$. To express the derivative with respect to $\hat{y}$, we apply the chain rule

[TABLE]

and use the fact that, for exponential families, the Jacobian matrix of the mean parameter $\frac{\partial\hat{y}}{\partial\beta}$ is equal to the covariance matrix $R_{t}$ of the sufficient statistics (Appendix, (A.11) and (A.6)). ∎

Lemma 8.

The extended Kalman filter satisfies $K_{t}R_{t}=P_{t}H_{t}^{\!\top}\!$.

Proof of the lemma.

.

This relation is known, e.g., [Sim06, (6.34)]. Indeed, using the definition of $K_{t}$, we have $K_{t}R_{t}=K_{t}(R_{t}+H_{t}P_{t|t-1}H_{t}^{\!\top}\!)-K_{t}H_{t}P_{t|t-1}H_{t}^{\!\top}\!=P_{t|t-1}H_{t}^{\!\top}\!-K_{t}H_{t}P_{t|t-1}H_{t}^{\!\top}\!=(\operatorname{Id}-K_{t}H_{t})P_{t|t-1}H_{t}^{\!\top}\!=P_{t}H_{t}^{\!\top}\!$. ∎

Proof of Proposition 6.

.

By definition of the Kalman filter we have $s_{t}=s_{t|{t-1}}+K_{t}E_{t}$. By Lemma 7, $\displaystyle E_{t}=R_{t}\left(\frac{\partial\ell_{t}}{\partial\hat{y}_{t}}\right)^{\!\top}\!$. Thanks to Lemma 8 we find $\displaystyle s_{t}=s_{t|{t-1}}+K_{t}R_{t}\left(\frac{\partial\ell_{t}}{\partial\hat{y}_{t}}\right)^{\!\top}\!=s_{t|{t-1}}+P_{t}H_{t}^{\!\top}\!\left(\frac{\partial\ell_{t}}{\partial\hat{y}_{t}}\right)^{\!\top}\!=s_{t|{t-1}}+P_{t}\left(\frac{\partial\ell_{t}}{\partial\hat{y}_{t}}H_{t}\right)^{\!\top}\!$. But by the definition of $H$, $H_{t}$ is $\displaystyle\frac{\partial\hat{y}_{t}}{\partial s_{t|{t-1}}}$ so that $\displaystyle\frac{\partial\ell_{t}}{\partial\hat{y}_{t}}H_{t}$ is $\displaystyle\frac{\partial\ell_{t}}{\partial s_{t|{t-1}}}$. ∎

The first part of the next lemma is known as the information filter in the Kalman filter literature, and states that the observation step for $P$ is additive when considered on $P^{-1}$ [Sim06, §6.2]: after each observation, the Fisher information matrix of the latest observation is added to $P^{-1}$.

Lemma 9 (Information filter).

The update (2.18)–(2.19) of $P_{t}$ in the extended Kalman filter is equivalent to

[TABLE]

(assuming $P_{t|t-1}$ and $R_{t}$ are invertible).

In particular, for static dynamical systems ($f(s,u)=s$ and $Q_{t}=0$), the whole extended Kalman filter (2.12)-(2.20) is equivalent to

[TABLE]

Proof.

.

The first statement is well-known for Kalman filters [Sim06, (6.33)]. Indeed, expanding the definition of $K_{t}$ in the update (2.19) of $P_{t}$, we have

[TABLE]

but this is equal to $(P_{t|t-1}^{-1}+H_{t}^{\!\top}\!R_{t}^{-1}H_{t})^{-1}$ thanks to the Woodbury matrix identity.

The second statement follows from Proposition 6 and the fact that for $f(s,u)=s$, the transition step of the Kalman filter is just $s_{t|t-1}=s_{t-1}$ and $P_{t|t-1}=P_{t-1}$. ∎

Lemma 10.

For exponential families $p(y|\hat{y})$, the term $H_{t}^{\!\top}\!R_{t}^{-1}H_{t}$ appearing in Lemma 9 is equal to the Fisher information matrix of $y$ with respect to the state $s$,

[TABLE]

where $\ell_{t}(y)=-\ln p(y|\hat{y}_{t})$ depends on $s$ via $\hat{y}=h(s,u)$.

Proof.

.

Let us omit time indices for brevity. We have $\displaystyle\frac{\partial\ell(y)}{\partial s}=\frac{\partial\ell(y)}{\partial\hat{y}}\frac{\partial\hat{y}}{\partial s}=\frac{\partial\ell(y)}{\partial\hat{y}}H$. Consequently, $\displaystyle\mathbb{E}_{y}\left[\frac{\partial\ell(y)}{\partial s}^{\otimes 2}\right]=H^{\!\top}\!\,\mathbb{E}_{y}\left[\frac{\partial\ell(y)}{\partial\hat{y}}^{\otimes 2}\right]H$. The middle term $\displaystyle\mathbb{E}_{y}\left[\frac{\partial\ell(y)}{\partial\hat{y}}^{\otimes 2}\right]$ is the Fisher matrix of the random variable $y$ with respect to $\hat{y}$.

Now, for an exponential family $y\sim p(y|\hat{y})$ in mean parameterization $\hat{y}$, the Fisher matrix with respect to $\hat{y}$ is equal to the inverse covariance matrix of the sufficient statistics of $y$ (Appendix, (A.16)), that is, $R_{t}^{-1}$. ∎

Proof of Theorem 2.

.

By induction on $t$. By the combination of Lemmas 9 and 10, the update of the Kalman filter with static dynamics ($s_{t|t-1}=s_{t-1}$) is

[TABLE]

Defining $J_{t}=P_{t}^{-1}/(t+1)$, this update is equivalent to

[TABLE]

Under the identification $s_{t-1}\leftrightarrow\theta_{t-1}$, this is the online natural gradient update with learning rate $\eta_{t}=1/(t+1)$ and metric update rate $\gamma_{t}=1/(t+1)$. ∎

The proof of Proposition 3 is similar, with additional factors $(1-\lambda_{t})$. Proposition 4 is proved by applying a fading memory Kalman filter to a modified log-likelihood $\bar{L}_{0}\mathrel{\mathop{:}}=n_{\mathrm{prior}}\ln\pi(\theta)$, $\bar{L}_{t}\mathrel{\mathop{:}}=\ln p_{\theta}(y_{t})+(1-\lambda_{t})\bar{L}_{t-1}+\lambda_{t}n_{\mathrm{prior}}\ln\pi(\theta)$ so that the prior is kept constant in $\bar{L}_{t}$.

3 Natural gradient as a Kalman filter: the state space (recurrent)

case

3.1 Recurrent models, RTRL

Let us now consider non-memoryless models, i.e., models defined by a recurrent or state space equation

[TABLE]

with $u_{t}$ the observations at time $t$. To save notation, here we dump into $\hat{y}_{t}$ the whole state of the model, including both the part that contains the prediction about $y_{t}$ and all state or internal variables (e.g., all internal and output layers of a recurrent neural network, not only the output layer). The state $\hat{y}_{t}$, or a part thereof, defines a loss function $\ell_{t}(y_{t})\mathrel{\mathop{:}}=-\ln p(y_{t}|\hat{y}_{t})$ for each observation $y_{t}$.

The current state $\hat{y}_{t}$ can be seen as a function which depends on $\theta$ via the whole trajectory. The derivative of the current state with respect to $\theta$ can be computed inductively just by differentiating the recurrent equation (3.1) defining $\hat{y}_{t}$:

[TABLE]

Real-time recurrent learning [Jae02] uses this equation to keep an estimate $G_{t}$ of $\frac{\partial\hat{y}_{t}}{\partial\theta}$. RTRL then uses $G_{t}$ to estimate the gradient of the loss function $\ell_{t}$ with respect to $\theta$ via the chain rule, $\partial\ell_{t}/\partial\theta=(\partial\ell_{t}/\partial\hat{y}_{t})(\partial\hat{y}_{t}/\partial\theta)=(\partial\ell_{t}/\partial\hat{y}_{t})G_{t}$.

Definition 11 (Real-time recurrent learning).

Given a recurrent model $\hat{y}_{t}=\Phi(\hat{y}_{t-1},\theta_{t-1},u_{t})$, real-time recurrent learning (RTRL) learns the parameter $\theta$ via

[TABLE]

Since $\theta$ changes at each step, the actual estimate $G_{t}$ in RTRL is only an approximation of the gradient $\frac{\partial\hat{y}_{t}}{\partial\theta}$ at $\theta=\theta_{t}$, valid in the limit of small learning rates $\eta_{t}$.

In practice, RTRL has a high computational cost due to the necessary storage of $G_{t}$, a matrix of size $\dim\theta\times\dim\hat{y}$. For large-dimensional models, backpropagation through time is usually preferred, truncated to a certain length in the past [Jae02]; [OTC15, TO17] introduce a low-rank, unbiased approximation of $G_{t}$.

3.2 Statement of the correspondence, recurrent case

There are several ways in which a Kalman filter can be used to estimate $\theta$ for such recurrent models.

A first possibility is to view each $\hat{y}_{t}$ as a function of $\theta$ via the whole trajectory, and to apply a Kalman filter on $\theta$. This would require, in principle, recomputing the whole trajectory from time [math] to time $t$ using the new value of $\theta$ at each step, and using RTRL to compute $\partial\hat{y}_{t}/\partial\theta$, which is needed in the filter. In practice, the past trajectory is not updated, and truncated backpropagation through time is used to approximate the derivatice $\partial\hat{y}_{t}/\partial\theta$ [Jae02, Hay01]. 2. 2.

A second possibility is the joint Kalman filter, namely, a Kalman filter on the pair $(\theta,\hat{y}_{t})$ [Hay01, §5], [Sim06, §13.4]. This does not require going back in time, as $\hat{y}_{t}$ is a function of $\hat{y}_{t-1}$ and $\theta$. This is the version appearing in Theorem 12 below. 3. 3.

A third possibility is the dual Kalman filter [WN96]: a Kalman filter for $\theta$ given $\hat{y}$, and another one for $\hat{y}$ given $\theta$. This requires to explicitly couple the two Kalman filters by manually adding RTRL-like terms to account for the (linearized) dependency of $\hat{y}$ on $\theta$ [Hay01, §5].

Intuitively, the joint Kalman filter maintains a covariance matrix on $(\theta,\hat{y}_{t})$, whose off-diagonal term is the covariance between $\hat{y}_{t}$ and $\theta$. This term captures how the current state would change if another value of the parameter had been used. The decomposition (3.13) in the theorem makes this intuition precise in relation to RTRL: the Kalman covariance between $\hat{y}_{t}$ and $\theta$ is directly given by the RTRL gradient $G_{t}$.

Theorem 12 (Kalman filter on $(\theta,\hat{y})$ as RTRL+natural

gradient+state correction).

Consider a recurrent model $\hat{y}_{t}=\Phi(\hat{y}_{t-1},\theta_{t-1},u_{t})$. Assume that the observations $y_{t}$ are predicted with a probabilistic model $p(y|\hat{y}_{t})$ that is an exponential family with mean parameter a subset of $\hat{y}_{t}$.

Given an estimate $G_{t}$ of $\partial\hat{y}_{t}/\partial\theta$, and an observation $y$, denote

[TABLE]

the corresponding estimate of $\partial\ell_{t}(y)/\partial\theta$.

Then these two algorithms are equivalent:

• •

The extended Kalman filter on the pair $(\theta,\hat{y})$ with transition function $(\operatorname{Id},\Phi)$, initialized with covariance matrix $P_{0}^{(\theta,\hat{y})}=\begin{pmatrix}P^{\theta}_{0}&0\\ 0&0\end{pmatrix}$, and with no process noise ($Q=0$).

• •

A natural gradient RTRL algorithm with learning rate $\eta_{t}=1/(t+1)$, defined as follows. The state, RTRL gradient and Fisher matrix have a transition step

[TABLE]

and after observing $y_{t}$, the state and parameter are updated as

[TABLE]

Moreover, at each time $t$, the covariance matrix of the extended Kalman filter over $(\theta,\hat{y})$ is related to $G_{t}$ and $J_{t}$ via

[TABLE]

This result may explain an observation from [Wil92, §4.2] that RTRL can be obtained by introducing some drastic simplifications in the Kalman filter equations (changing the formula of the Kalman optimal gain and neglecting the covariance matrix update).

Again, the expectation for the Fisher matrix in (3.9) may be estimated by a Monte Carlo sample $y\sim p(y|\hat{y}_{t})$, or by just using the current observation $y=y_{t}$, as discussed after Definition 1.

As before, learning rates $\eta_{t}$ different from $1/(t+1)$ can be obtained by introducing a fading memory (i.e., process noise $Q$ proportional to $P$) in the joint Kalman filter. We omit the statement for simplicity, but it is analogous to Propositions 3 and 4.

The algorithm above features a state update (3.12) together with the parameter update; this is not commonly used in online recurrent neural network algorithms. In small-scale experiments, we have not found any clear effect of this; besides, such state updates must be applied cautiously if the range of possible values for the state $\hat{y}$ is somehow constrained.

In the result above, the Kalman filter is initialized with a covariance matrix in which every uncertainty comes from uncertainty on $\theta$ rather than the initial state $\hat{y}_{0}$. This has the advantage of making the correspondence algebraically simple, but is not a fundamental restriction. If modelling an initial uncertainty on $\hat{y}_{0}$ is important, one can always apply the theorem by incorporating the initial condition as an additional component of the parameter $\theta$, with its own variance; in this case, $G_{0}$ must be initialized to $\operatorname{Id}$ on the corresponding component of $\theta$, namely

[TABLE]

and then Theorem 12 can be applied to $\theta^{+\mathrm{init}}$.

Actually this operation is often not needed at all: indeed, if the dynamical system is such that the initial condition is forgotten reasonably quickly, then the initial covariance of $\hat{y}_{0}$ decreases (terms $W$ in the proof below) and the Kalman covariance tends to the type (3.13) above exponentially fast, even without using $\theta^{+\mathrm{init}}$. This is the case, for instance, for any stable linear dynamical system, as a consequence of Lemmas 13-14, and more generally for any system with geometric memory in the sense that $\frac{\partial\hat{y}_{t}}{\partial\hat{y}_{t-1}}$ is contracting for a fixed parameter and a given input.

The filtering literature contains updates similar to the above for $G_{t}$, but more complex [LS83, Hay01]; this is, first, because they are expressed over the variable $\operatorname{Cov}(\hat{y}_{t},\theta)=G_{t}J_{t}^{-1}$ instead of $G_{t}$ alone, second, because we have initialized the uncertainty on $\hat{y}_{0}$ to [math], and, third, because in dual rather than joint filter approaches, higher-order terms depending on second derivatives of $F$ are sometimes included. Interestingly, there is some debate in the literature about whether to add some second-order corrections to the joint Kalman filter (especially [LS83, §2.3.3], see discussion in [Hay01, §5.3.4]). The interpretation in Theorem 12 makes it clear which terms are neglected: in particular, in RTRL $G_{t}$ is not recomputed after the update of $\theta$ and $\hat{y}_{t}$, so that $G_{t}$ contains a mixture of derivatives at different values of $\theta$ over time. Correcting for this would involve second derivatives of $F$ (as in [Hay01, §5, Appendix A]), thus amounting to a partial implementation of a second-order extended Kalman filter (EKF2, [Sim06, §13.3]).

In terms of computational cost, for recurrent neural networks (RNNs), RTRL alone is already as costly as the joint Kalman filter [Wil92]. Indeed, RTRL requires $(\dim\theta)$ forward tangent propagations at each step, each of which costs $O(\dim\theta)$ for a standard RNN model [Jae02], thus for a total cost of $O((\dim\theta)^{2})$ per time step. The Fisher matrix is of size $(\dim\theta)^{2}$; if a single Monte Carlo sample $y\sim p(y|\hat{y}_{t})$ is used, then the Fisher matrix update is rank-one and costs $O((\dim\theta)^{2})$; the update of the inverse Fisher matrix can be maintained at the same cost thanks to the Woodbury matrix identity (as done, e.g., in [LMB07]). Thus, if RTRL is computationally affordable, there is little point in not using the Fisher matrix on top.

3.3 Proofs for the recurrent case

We now turn to the proof for the recurrent case, involving a joint Kalman filter on $(\theta,\hat{y})$. The key is to decompose the Kalman covariance matrix of the pair $(\theta,\hat{y})$ into three variables (3.17): the covariance of $\theta$, the correlation between $\theta$ and $\hat{y}$, and the part of the covariance of $\hat{y}$ that does not come from its correlation with $\theta$ (its so-called Schur complement). This provides a nice expression for the transition step of the Kalman filter (Lemma 13).

Then we find that the correlation between $\theta$ and $\hat{y}$ is exactly the gradient $G=\frac{\partial\hat{y}}{\partial\theta}$ maintained by RTRL (Corollary 15); meanwhile, we find $\theta$ and its covariance essentially follow a standalone Kalman filter related to the observations via $G$, which is a natural gradient for the same reasons as in the static case.

In the recurrent case, we are applying an extended Kalman filter to the state $s=\begin{pmatrix}\theta\\ \hat{y}\end{pmatrix}$ with transition function $f=\begin{pmatrix}\operatorname{Id}\\ \Phi\end{pmatrix}$. Let us decompose the covariance matrix $P_{t}$ of this system as

[TABLE]

From now on, for simplicity we omit the time indices when no ambiguity is present.

By the theory of Schur complement for positive-semidefinite matrices [BV04, Appendix A.5.5], letting $P^{+}$ be any generalized inverse of $P^{\theta}$, we know that $P^{\hat{y}}-P^{\theta\hat{y}}P^{+}{P^{\theta\hat{y}}}^{\!\top}\!$ is positive-semidefinite and that $P^{\theta\hat{y}}(\operatorname{Id}-P^{+}P^{\theta})=0$. The latter rewrites as $P^{\theta\hat{y}}=P^{\theta\hat{y}}P^{+}P^{\theta}$. Let us set

[TABLE]

Then $P^{\theta\hat{y}}=GP^{\theta}$ and $W=P^{\hat{y}}-GP^{\theta}G^{\!\top}\!$. Thus, at each time $t$ we can decompose $P_{t}$ as

[TABLE]

without loss of generality, where $W$ is positive-semidefinite. This decomposition tells us which part of the covariance of the current state $\hat{y}$ comes from the covariance of the parameter $\theta$ via the dynamics of the system.

First, we will show that if $W_{0}=0$, then $W_{t}=0$ for all $t$, and that in this case $G_{t}$ satisfies the RTRL equation.

Lemma 13.

Consider the extended Kalman filter on the pair $s=(\theta,\hat{y})^{\!\top}\!$ with transition function $f=(\theta,\Phi(\hat{y},\theta,u))^{\!\top}\!$ and no added noise ($Q_{t}=0$). Then the Kalman transition step (2.13) on $P$, expressed in the decomposition (3.17), is equivalent to

[TABLE]

This equation for $G$ is the RTRL update.

Proof of the lemma.

.

This is a direct computation using the Kalman transition step (2.13) for $P$. Indeed, the decomposition (3.17) of $P$ rewrites as

[TABLE]

Now, the Kalman transition step (2.13) for $P$ is $P_{t|t-1}=\frac{\partial f}{\partial s}P_{t-1}\frac{\partial f}{\partial s}^{\!\top}\!$ when $Q=0$. So the update is equivalent to replacing $\begin{pmatrix}\operatorname{Id}&0\\ G&\operatorname{Id}\end{pmatrix}$ with $\frac{\partial f}{\partial s}\begin{pmatrix}\operatorname{Id}&0\\ G&\operatorname{Id}\end{pmatrix}$ on both sides in (3.21). Here we have $\frac{\partial f}{\partial s}=\begin{pmatrix}\operatorname{Id}&0\\ \frac{\partial\Phi}{\partial\theta}&\frac{\partial\Phi}{\partial\hat{y}}\end{pmatrix}$. This yields the result. ∎

Lemma 14.

Consider the extended Kalman filter on the pair $s=(\theta,\hat{y})^{\!\top}\!$ with transition function $f=(\theta,\Phi(\hat{y},\theta,u))^{\!\top}\!$. Then the observation update (2.18)–(2.19) of $P_{t}$, expressed in the variables $P^{\theta}$, $W$, and $G$, is given by

[TABLE]

in that order, where $R$ is given by (2.16). Moreover, if $P^{\theta}$ or $W$ are invertible then their respective updates are equivalent to

[TABLE]

and

[TABLE]

Thus, the updates for $W$, (3.19) and (3.23), are just the updates of an extended Kalman filter on $\hat{y}$ alone, with covariance matrix $W$ and noise measurement $R$. The update for $P^{\theta}$ is identical to an extended Kalman filter on $\theta$ where measurements are made on $\hat{y}$, with $\hat{y}$ seen as a function of $\theta$ with derivative $\partial\hat{y}/\partial\theta=G$, and where the measurement noise on $\hat{y}$ is $R+W$ (the measurement noise on $y$ plus the covariance of $\hat{y}$). Thus, these two lemmas relate the joint Kalman filter on $(\theta,\hat{y})$ to the dual Kalman filter that filters separately $\theta$ given $\hat{y}$ and $\hat{y}$ given $\theta$, together with an estimate of $\partial\hat{y}/\partial\theta$. As far as we could check, this decomposition is specific to a situation in which one component (the parameter) is supposed to have static underlying dynamics, $\theta_{t+1}=\theta_{t}$.

Proof of the lemma.

.

In our case, the function $h$ of the extended Kalman filter is the function that sends $(\theta,\hat{y})$ to $\hat{y}$. In particular, $H_{t}=(0,\operatorname{Id})$.

First, if $P^{\theta}$ and $W$ are invertible, then the updates (3.22), (3.23) for $P^{\theta}$ and $W$ follow from the updates (3.25), (3.26) on their inverses, thanks to the Woodbury matrix identity. Since working on the inverses is simpler, we shall prove only the latter. Since (3.22), (3.23) are continuous in $P^{\theta}$ and $W$, the non-invertible case follows by continuity.

Starting again with the decomposition of $P_{t}$ as a product (3.21), the inverse of $P_{t}$ is

[TABLE]

From Lemma 9, the Kalman observation udpate for $P_{t}$ amounts to adding $H^{\!\top}\!R^{-1}H$ to $P_{t}^{-1}$. Here $H=(0,\operatorname{Id})$ so that $H^{\!\top}\!R^{-1}H$ is $\begin{pmatrix}0&0\\ 0&R^{-1}\end{pmatrix}$. So the update for $P_{t}$ amounts to

[TABLE]

To interpret this as an update on $P^{\theta}$, $W$ and $G$, we have to introduce new variables $\tilde{W}$, $\tilde{G}$, and $\tilde{P}^{\theta}$ such that (3.29) takes the original form (3.28) in these new variables.

Introducing $\tilde{W}^{-1}\mathrel{\mathop{:}}=W^{-1}+R^{-1}$, the update rewrites as

[TABLE]

Introducing $\tilde{G}\mathrel{\mathop{:}}=(\operatorname{Id}-\tilde{W}R^{-1})G$ and $(\tilde{P}^{\theta})^{-1}\mathrel{\mathop{:}}=(P^{\theta})^{-1}+G^{\!\top}\!W^{-1}G-\tilde{G}^{\!\top}\!\tilde{W}^{-1}\tilde{G}$, we get back the original form (3.28). This provides the updates $\tilde{W}$ and $\tilde{G}$ for $W$ and $G$. We still have to find a more explicit expression for $(\tilde{P}^{\theta})^{-1}$.

Since we have defined $\tilde{W}$ and $\tilde{G}$ by identifying (3.29) with the original form (3.28), we have $\tilde{W}\tilde{G}=WG$ by construction. Thus

[TABLE]

Thanks to the identity $A^{-1}-A^{-1}BA^{-1}=(A+(B^{-1}-A^{-1})^{-1})^{-1}$ for any matrices $A$ and $B$ (this follows from the matrix inversion formula $(A+C)^{-1}=A^{-1}-A^{-1}(A^{-1}+C^{-1})^{-1}A^{-1}$ applied to $C=(B^{-1}-A^{-1})^{-1}$), we find

[TABLE]

but by definition, $\tilde{W}^{-1}=W^{-1}+R^{-1}$ so that

[TABLE]

thus

[TABLE]

which concludes the proof of the lemma. ∎

Putting the last two lemmas side by side in the case $W=0$, we obtain a much simpler update.

Corollary 15.

Consider the extended Kalman filter on the pair $s=(\theta,\hat{y})^{\!\top}\!$ with transition function $f(s)=(\theta,\Phi(\hat{y},\theta,u))^{\!\top}\!$ and no added noise ($Q_{t}=0$). Decompose the covariance $P$ of the state $s$ as in (3.17) using $P^{\theta}$, $G$, $W$. If $W=0$ and $P^{\theta}$ is invertible then performing the Kalman transition update followed by the observation update is equivalent to

[TABLE]

in that order.

From this, the end of the proof of Theorem 12 essentially proceeds as in the non-recurrent case. Since we initialize $W$ to [math] in Theorem 12, we have $W=0$ at all times. As before, for exponential families $R^{-1}$ is equal to the Fisher matrix with respect to $\hat{y}_{t}$, namely, $R^{-1}=\mathbb{E}_{y\sim p(y|\hat{y})}\displaystyle\left[\frac{\partial\ell_{t}(y)}{\partial\hat{y}_{t}}^{\otimes 2}\right]$ (Appendix). Now, the term $\mathbb{E}_{y\sim p(y|\hat{y})}\left[g_{t}(y)^{\otimes 2}\right]$ in the Fisher matrix update (3.9) uses $g_{t}(y)=\displaystyle\frac{\partial\ell_{t}(y)}{\partial\hat{y}_{t}}\,G_{t}$ (3.6) to estimate the derivative of the loss $\ell_{t}(y)$ with respect to $\theta$. So the term $G^{\!\top}\!R^{-1}G$ in (3.38) coincides with the Fisher matrix update term $\mathbb{E}_{y\sim p(y|\hat{y})}\left[g_{t}(y)^{\otimes 2}\right]$ in (3.9). (Compare Lemma 10.) So if we just define $J_{t}\mathrel{\mathop{:}}=\eta_{t}(P^{\theta}_{t})^{-1}$ with $\eta_{t}=1/(t+1)$, the additive update (3.38) on $(P^{\theta})^{-1}$ translates as the online Fisher matrix update (3.9) on $J_{t}$.

Moreover, since the Kalman gradient is an ordinary gradient preconditioned with the covariance matrix $P_{t}$ (Proposition 6), the update of the pair $(\theta,\hat{y})$ is

[TABLE]

(indeed $\ell_{t}$ does not depend explicitly on $\theta$ in recurrent models, only via the current state $\hat{y}_{t}$). Given the decomposition $P_{t}=\begin{pmatrix}P^{\theta}&(GP^{\theta})^{\!\top}\!\\ GP^{\theta}&GP^{\theta}G^{\!\top}\!\end{pmatrix}$, this translates as

[TABLE]

which is the update in Theorem 12.

Appendix A Appendix: reminder on exponential families

An exponential family of probability distributions on a variable $x$ (discrete or continuous), with sufficient statistics $T_{1}(x),\ldots,T_{K}(x)$, is the following family of distributions, parameterized by $\beta\in{\mathbb{R}}^{K}$:

[TABLE]

where $Z(\beta)$ is a normalizing constant, and $\lambda(\operatorname{d}\!{}x)$ is any reference measure on $x$, such as the Lebesgue measure or any discrete measure. The family is obtained by varying the parameter $\beta\in{\mathbb{R}}^{K}$, called the natural or canonical parameter. We will assume that the $T_{k}$ are linearly independent as functions of $x$ (and linearly independent from the constant function); this ensures that different values of $\beta$ yield distinct distributions.

For instance, Bernoulli distributions are obtained with $\lambda$ the uniform measure on $x\in\{0,1\}$ and with a single sufficient statistic $T(0)=0$, $T(1)=1$. Gaussian distributions with a fixed variance are obtained with $\lambda(\operatorname{d}\!{}x)$ the Gaussian distribution centered on [math], and $T(x)=x$.

Another, often convenient parameterization of the same family is the following: each value of $\beta$ gives rise to an average value $\bar{T}$ of the sufficient statistics,

[TABLE]

For instance, for Gaussian distributions with fixed variance, this is the mean, and for a Bernoulli variable this is the probability to sample $1$.

Exponential families satisfy the identities

[TABLE]

by a simple computation [AN00, (2.33)].

These identities are useful to compute the Fisher matrix $J_{\beta}$ with respect to the variable $\beta$, as follows [AN00, (3.59)]:

[TABLE]

or more synthetically

[TABLE]

where the covariance is under the law $p_{\beta}$. That is, for exponential families the Fisher matrix is the covariance matrix of the sufficient statistics. In particular it can be estimated empirically, and is sometimes known algebraically.

In this work we need the Fisher matrix with respect to the mean parameter $\bar{T}$,

[TABLE]

By substituting $\frac{\partial\ln p(x)}{\partial\alpha}=\frac{\partial\ln p(x)}{\partial\beta}\frac{\partial\beta}{\partial\alpha}$, the Fisher matrices $J_{\alpha}$ and $J_{\beta}$ with respect to parameterizations $\alpha$ and $\beta$ are related to each other via

[TABLE]

(consistently with the interpretation of the Fisher matrix as a Riemannian metric and the behavior of metrics under change of coordinates [GHL87, §2.3]). So we need to compute $\partial{\bar{T}}/\partial\beta$. Using the log-trick

[TABLE]

together with (A.3), we find

[TABLE]

so that

[TABLE]

(see [AN00, (3.32)], where $\eta$ denotes the mean parameter) and consequently

[TABLE]

so that we find the Fisher matrix with respect to $\bar{T}$ to be

[TABLE]

that is, the Fisher matrix with respect to $\bar{T}$ is the inverse covariance matrix of the sufficient statistics.

This gives rise to a simple formula for the natural gradient of expectations with respect to the mean parameters. Denoting $\tilde{\nabla}$ the natural gradient,

[TABLE]

which in particular, can be estimated empirically.