Maximum Likelihood With a Time Varying Parameter

TL;DR

This paper introduces a recursive stochastic gradient descent method for tracking a time-varying parameter in probabilistic models, proving convergence and illustrating its effectiveness with exponential family distributions.

Contribution

It proposes a novel recursive scheme using the log-likelihood as a gain function for tracking time-varying parameters, with proven convergence properties.

Findings

Convergence in mean-square error near the true parameter

Effective tracking demonstrated with exponential family distributions

The method adapts to changing parameters in real-time

Abstract

We consider the problem of tracking an unknown time varying parameter that characterizes the probabilistic evolution of a sequence of independent observations. To this aim, we propose a stochastic gradient descent-based recursive scheme in which the log-likelihood of the observations acts as time varying gain function. We prove convergence in mean-square error in a suitable neighbourhood of the unknown time varying parameter and illustrate the details of our findings in the case where data are generated from distributions belonging to the exponential family.

Equations

$${\lambda }_{t+1}:={\lambda }_t+\alpha {\nabla }_{\lambda }\ln p(X_t|{\lambda }_t),t\in \mathbb{N}.$$

$${\mathbb{R}}^d\ni \lambda \mapsto \ln p(x|\lambda )$$

$${\partial }_{{\lambda }_i}{\partial }_{{\lambda }_j}\mathbb{E}\left[\ln p(X_t|\lambda )\right]=\mathbb{E}\left[{\partial }_{{\lambda }_j}{\partial }_{{\lambda }_j}\ln p(X_t|\lambda )\right],$$

$${\mathbb{R}}^d\ni \lambda \mapsto {\nabla }_{\lambda }\ln p(x|\lambda )$$

$$\Vert {\nabla }_{\lambda }\ln p(x|{\xi }_1)-{\nabla }_{\lambda }\ln p(x|{\xi }_2)\Vert \le L\Vert {\xi }_1-{\xi }_2\Vert ,{\xi }_1,{\xi }_2\in {\mathbb{R}}^d.$$

$$\mathtt{I}({\lambda }_t^{\ast })\le dL,$$

$$\mathtt{I}({\lambda }_t^{\ast })$$

$$=\sum _{j=1}^d\mathbb{E}[{\partial }_{{\lambda }_j}^2(-\ln p(X_t|{\lambda }_t^{\ast }))]=\sum _{j=1}^d\mathbb{E}[⟨{\mathcal{H}}_{\lambda }(-\ln p(X_t|{\lambda }_t^{\ast }))e_j,e_j⟩]$$

$$\le \sum _{j=1}^d\mathbb{E}[⟨L{I}_d{e}_j,e_j⟩]=dL.$$

$$\displaystyle\|\lambda_{t+1}^{*}-\lambda_{t}^{*}\|\leq K\quad\mbox{ for all $t\in\mathbb{N}$}.$$

$$\underset{t\to +\infty }{\mathrm{limsup}}\Vert {\lambda }_{t+1}-{\lambda }_t^{\ast }{\Vert }_{{\mathbb{L}}^2(\Omega )}\le \frac{\phi (\alpha ,L)K+\alpha \sqrt{2dL}}{1-\phi (\alpha ,L)},$$

$$\underset{t\to +\infty }{\mathrm{limsup}}\Vert {\lambda }_{t+1}-{\lambda }_t^{\ast }{\Vert }_{{\mathbb{L}}^2(\Omega )}\le \frac{K\sqrt{{\ell }^2+L^2}+\sqrt{2dL}}{\ell +L-\sqrt{{\ell }^2+L^2}}.$$

$${\tilde{\lambda }}_t:=\arg \underset{\lambda \in {\mathbb{R}}^d}{\max }\mathbb{E}[\ln p(X_t|\lambda )].$$

$$p(x|\lambda )=h(x)\exp \{⟨\lambda ,T(x)⟩-A(\lambda )\},x\in {\mathbb{R}}^m$$

$${\nabla }_{\lambda }\ln p(x|\lambda )=T(x)-{\nabla }_{\lambda }A(\lambda ),$$

$${\mathcal{H}}_{\lambda }[-\ln p(x|\lambda )]=-{\mathcal{H}}_{\lambda }A(\lambda ),$$

$${\bar{\lambda }}_{t+1}:={\Pi }_{\Lambda }\left({\bar{\lambda }}_t+\alpha {\nabla }_{\lambda }\ln p(X_t|{\bar{\lambda }}_t)\right),t\in \mathbb{N},$$

$$\Vert {\bar{\lambda }}_{t+1}-{\lambda }_t^{\ast }{\Vert }^2=\Vert {\Pi }_{\Lambda }({\bar{\lambda }}_t+\alpha {\nabla }_{\lambda }\ln p(X_t|{\bar{\lambda }}_t))-{\lambda }_t^{\ast }{\Vert }^2\le \Vert {\bar{\lambda }}_t+\alpha {\nabla }_{\lambda }\ln p(X_t|{\bar{\lambda }}_t)-{\lambda }_t^{\ast }{\Vert }^2,$$

$$X_t\sim p(x|{\lambda }_t^{\ast }):=h(x)\exp \{{\lambda }_t^{\ast }T(x)-A({\lambda }_t^{\ast })\},x\in \mathbb{R}.$$

$$\ell =\underset{\lambda \in \Lambda }{\min }{A}^{\prime \prime }(\lambda ),L=\underset{\lambda \in \Lambda }{\max }{A}^{\prime \prime }(\lambda ).$$

$$\sqrt{\mathbb{E}[|T(X_t)-{\lambda }_t^{\ast }{|}^2]}=\sqrt{\mathbb{V}[T(X_t)]}=\sqrt{A^{\prime \prime }({\lambda }_t^{\ast })},$$

$$\displaystyle\frac{K\sqrt{\ell^{2}+L^{2}}+\sqrt{2L}}{\ell+L-\sqrt{\ell^{2}+L^{2}}}\leq\sqrt{A^{\prime\prime}(\lambda_{t}^{*})}\quad\mbox{ for all $t\in\mathbb{N}$}.$$

$$\frac{K\sqrt{{\ell }^2+L^2}+\sqrt{2L}}{\ell +L-\sqrt{{\ell }^2+L^2}}\le \sqrt{\ell }.$$

$$\Vert {\lambda }_{t+1}-{\lambda }_t^{\ast }{\Vert }^2=$$

$$=$$

$$+2\alpha ⟨{\lambda }_t-{\lambda }_t^{\ast },{\nabla }_{\lambda }\ln p(X_t|{\lambda }_t^{\ast })⟩+{\alpha }^2\Vert {\nabla }_{\lambda }\ln p(X_t|{\lambda }_t){\Vert }^2$$

$$=$$

$${\mathcal{A}}_1:=2\alpha ⟨{\lambda }_t-{\lambda }_t^{\ast },{\nabla }_{\lambda }\ln p(X_t|{\lambda }_t)-{\nabla }_{\lambda }\ln p(X_t|{\lambda }_t^{\ast })⟩$$

$${\mathcal{A}}_2:={\alpha }^2\Vert {\nabla }_{\lambda }\ln p(X_t|{\lambda }_t){\Vert }^2.$$

$${\mathcal{A}}_1\le$$

$${\mathcal{A}}_2\le 2{\alpha }^2\Vert {\nabla }_{\lambda }\ln p(X_t|{\lambda }_t)-{\nabla }_{\lambda }\ln p(X_t|{\lambda }_t^{\ast }){\Vert }^2+2{\alpha }^2\Vert {\nabla }_{\lambda }\ln p(X_t|{\lambda }_t^{\ast }){\Vert }^2.$$

$$\Vert {\lambda }_{t+1}-{\lambda }_t^{\ast }{\Vert }^2\le$$

Full text

Key words and phrases: stochastic gradient descent, maxim likelihood, exponential family.

AMS 2020 classification: 65K05, 62F12.

1 Introduction

When estimating unknown parameters in a dynamic model the optimum solution to the parameter estimation problem may not remain constant. Specifically, the optimal values of the model parameters may change through time because of the evolution of the underlying process: finding them is, in general, not straightforward. A survey of basic techniques for tracking the time-varying dynamics of a system is provided in [Ljung and Gunnarsson, 1990] where recursive algorithms in non-stationary stochastic optimization are analysed under different assumptions about the true system’s variations, see also [Simonetto et al., 2020] for a review in a purely deterministic setting. In [Delyon and Juditsky, 1995] the problem of tracking the random drifting parameters of a linear regression system is tackled, and [Zhu and Spall, 2016] builds a computable tracking error bound for how a stochastic approximation with constant gain keeps up with a non-stationary target. Successively, [Wilson et al., 2019] introduces a framework for sequentially solving convex stochastic minimization problems, where the distance between successive minimizers is bounded. The minimization problems are then solved by sequentially applying an optimization algorithm, such as stochastic gradient descent (SGD). In a similar setting, [Cao et al., 2019] establishes an upper bound on the regret of a projected SGD algorithm with respect to the drift of the dynamic optima, while [Cutler et al., 2021] provides novel non-asymptotic convergence guarantees for stochastic algorithms with iterate averaging.

We study time-varying stochastic optimization in a general statistical setting where we assume we are given a sequence of independent observations $\{X_{t}\}_{t\in\mathbb{N}}$ with associated densities possessing a parameter that changes through time. In such a framework a problem of interest concerns finding a useful estimator of the time varying parameter at a certain time $t$ - generalizing the classical problem of parameter estimation from the static setting to the time varying parameter setting. Ideally, one would like to find a sequence of estimators that track the time varying parameter through time as closely as possible. We show that, under some assumptions, utilizing the celebrated SGD algorithm [Robbins and Monro, 1951] produces a sequence of estimators that will eventually track the time varying parameter - up to a neighborhood - as the number of observations increase.

Established in a general setting that intersects with the frameworks utilized in [Cao et al., 2019], [Cutler et al., 2021] and [Wilson et al., 2019], our results differ from previous work mainly in one aspect: that our objective functions have the specific form of expected log likelihoods, a dissimilarity that will be exploited by utilizing their informational theoretical properties.

The work we present is also linked to the class of score driven models [Creal et al., 2013]. Score driven models are a class of observation driven models (here we are using the terminology introduced by [Cox et al., 1981]) that update the dynamics of the time varying parameter through the score of the conditional distribution of the observations. Specifically, the same proof technique we utilize to obtain our result can be used to show that a -so called- Newton-score update [Blasques et al., 2015], with the parameter that multiplies the score appropriately chosen, will track the time varying parameter of interest trough time even under possible model misspecificaiton.

A final way to interpret the results we present in this work is as robustness results for a one batch stochastic gradient procedure in the case we are incorrectly assuming that our observations are identically distributed. Indeed, the results show that even if we incorrectly assumed that the true parameter is static (we have IID observations) utilizing a stochastic gradient algorithm with a time dependent single sized batch to optimize the log-likelihood allows us to track the pseudo true time varying parameter up to a neighborhood if it is not moving wildly.

The paper is organised as follows: in Section 2 we list and discuss the assumptions of our framework and state the main result. We then present a class of examples given by the exponential family and discuss the performance of SGD with respect to the one observation maximum likelihood estimator at each time. In the third section we provide a detailed proof of our main result.

2 Statement of the main result

Let $\{X_{t}\}_{t\in\mathbb{N}}$ be a sequence of independent $m$-dimensional random vectors defined on a common probability space $(\Omega,\mathcal{F},\mathbb{P})$. In the sequel we will write $\mathbb{E}[\cdot]$ for the expected value with respect to the probability measure $\mathbb{P}$, $\|\cdot\|$ for the Euclidean norm in $\mathbb{R}^{d}$ and $\|\cdot\|_{\mathbb{L}^{2}(\Omega)}$ for $\mathbb{E}[\|\cdot\|^{2}]^{\frac{1}{2}}$.

We assume that for any $t\in\mathbb{N}$ the random vector $X_{t}$ possesses a joint probability density function which depends on the $d$-dimensional parameter $\lambda_{t}^{*}$, in symbols $X_{t}\thicksim p(\cdot|\lambda_{t}^{*})$. Our aim is to estimate the sequence $\{\lambda_{t}^{*}\}_{t\in\mathbb{N}}$ through the observed values $\{X_{t}\}_{t\in\mathbb{N}}$: To this aim we choose $\lambda_{1}\in\mathbb{R}^{d}$ and utilize the SGD algorithm

[TABLE]

Utilizing SGD to attempt to track $\lambda_{t}^{*}$ is motivated by the principle underlying classical maximum likelihood estimation: in fact, under some canonical assumptions we will present below, $\lambda_{t}^{*}$ will be the maximum of the expected log-likelihood $\lambda\rightarrow\mathbb{E}\left[\ln p(X_{t}|\lambda)\right]$. Thus, finding a sequence of estimators that track the time varying parameter as closely as possible is connected to finding the maxima of a sequence of expected log-likelihoods, a generalization of the classical static framework. Since we have no direct access to the expected log-likelihoods, but only a singe observation for each time $t$, we categorize the problem as a time varying stochastic optimization problem.

The assumptions we will require to obtain our result are the following.

Assumption 2.1** (Smoothness of the log-likelihood).**

The function

[TABLE]

is twice continuously differentiable for all $x\in\mathbb{R}^{m}$; moreover,

[TABLE]

for all $i,j\in\{1,...,d\}$ and $t\in\mathbb{N}$.

Assumption 2.2** (Strong convexity).**

The function in (2.2) is strongly convex uniformly with respect to $x\in\mathbb{R}^{m}$: i.e., there exists a positive constant $\ell$ such that for all $x\in\mathbb{R}^{m}$ the matrix $\mathcal{H}_{\lambda}[-\ln p(x|\lambda)]-\ell I_{d}$ is positive semi-definite. Here, $\mathcal{H}_{\lambda}[-\ln p(x|\lambda)]$ stands for the Hessian matrix of the function in (2.2) while $I_{d}$ denotes the $d\times d$ identity matrix.

Assumption 2.3** (Lipschitz continuity of the gradient).**

The function

[TABLE]

is globally Lipschitz continuous uniformly with respect to $x\in\mathbb{R}^{m}$: i.e., there exists a positive constant $L$ such that for all $x\in\mathbb{R}^{m}$ we have

[TABLE]

Assumptions 2.2 and 2.3 are classical in the optimization literature, see for instance [Boyd and Vandenberghe, 2004] and[Bottou et al., 2018]; we have utilized the versions of [Nesterov, 2014]. We remark that Assumption 2.2 may seem excessively restrictive at first glance, but we will present in Example 2.9 below a large family of examples where it holds.

Remark 2.4**.**

Assumptions 2.1 and 2.3 imply that

[TABLE]

where we have denoted $\mathtt{I}(\lambda_{t}^{*}):=\mathbb{E}[\|\nabla_{\lambda}\ln p(X_{t}|\lambda_{t}^{*})\|^{2}]$, i.e. the trace of Fisher information matrix of $X_{t}$. In fact,

[TABLE]

We will use Remark 2.4 to bound the quantity $\mathbb{E}[\|\nabla_{\lambda}\ln p(X_{t}|\lambda_{t})\|^{2}]$. In the general setting utilized in the optimization literature a bound on $\mathbb{E}[\|\nabla_{\lambda}\ln p(X_{t}|\lambda_{t})\|^{2}]$ requires an extra assumption, see [Bottou et al., 2018] and the discussion in [Nguyen et al., 2018]. In our setting we manage to avoid this type of additional assumption thanks to the properties of the Fisher information matrix.

Our last assumption concerns the evolution of the time varying parameter $\{\lambda_{t}^{*}\}_{t\in\mathbb{N}}$.

Assumption 2.5** (Lipschitz continuity of the true parameter).**

There exists a positive constant $K$ such that

[TABLE]

Assumption 2.5 has been used throughout the literature, see for example [Simonetto et al., 2020], [Cao et al., 2019] and [Wilson et al., 2019], since a limitation on the behavior of the sequence of true parameters values must be imposed to be able to track it.

We can now state our main theorem.

Theorem 2.6**.**

Let Assumptions 2.1, 2.2, 2.3 and 2.5 hold. Then, for $\alpha\in[\frac{1}{\ell+L},\frac{1}{L}[$ running the SGD (2.1) we obtain

[TABLE]

where $\varphi(\alpha,L):=\sqrt{1-2L\alpha+2L^{2}\alpha^{2}}$. Moreover, the minimum of the right hand side in (2.3) is attained at $\alpha=\frac{1}{\ell+L}$ and in this case the last inequality reads

[TABLE]

Remark 2.7**.**

Notice that $\lambda_{t+1}$ depends on $X_{1},X_{2},\dots,X_{t}$, so as an estimator it is natural to compare it with $\lambda_{t}^{*}$.

Remark 2.8**.**

In the case of model misspecification, i.e. when the true distribution of the observations is not included in the parametric model $\{p(\cdot|\lambda)\}_{\lambda\in\mathbb{R}^{d}}$, the same proof technique can be utilized to show that the recursion (2.1) will track the so called pseudo-true time varying parameter $\tilde{\lambda}_{t}$ which is defined as

[TABLE]

*We recall that the pseudo-true time varying parameter $\tilde{\lambda}_{t}$ minimizes the Kullback Leiber divergence between the law of the data generating process and the model densities at each time $t$, see [White, 1982] and [Akaike, 1973] for additional details.

The only technical difference in the proof is that Remark 2.4 can’t be used since $\mathbb{E}[\|\nabla_{\lambda}\ln p(X_{t}|\tilde{\lambda}_{t})\|^{2}]$ is no longer related to the Fisher information matrix of $X_{t}$. Thus, an additional assumption is needed to control $E[\|\nabla_{\lambda}\ln p(X_{t}|\tilde{\lambda}_{t})\|^{2}]$ but this is standard practice in the optimization literature, see [Nguyen et al., 2018] for a discussion on this kind of assumption.*

Example 2.9**.**

The exponential family in canonical form provides a class of natural examples where Theorem 2.6 holds. Take as the parameter of interest the natural parameter of a distribution belonging to the exponential family put in canonical form, i.e.

[TABLE]

*where $h:\mathbb{R}^{m}\rightarrow\mathbb{R}$ is a non-negative function, $T:\mathbb{R}^{m}\rightarrow\mathbb{R}^{d}$ is a sufficient statistic and $A:\mathbb{R}^{d}\rightarrow\mathbb{R}$ must be chosen so that $p(x|\lambda)$ integrates to one.

A standard result for exponential families, see for instance Theorem 1.6.3 in [Bickel and Doksum, 2001], is that $A$ is a convex function of $\lambda$; this fact together with identities*

[TABLE]

and

[TABLE]

*implies that one can find, restricting if necessary the range of $\lambda$ (and hence of $\{\lambda_{t}^{*}\}_{t\in\mathbb{N}}$) to a suitable convex compact set $\Lambda$, the positive constants $l$ and $L$ required for the validity of Assumptions 2.2-2.3.

Note that the restriction of the range of $\lambda$ to the convex compact set $\Lambda$ is carried out by simply modifying (2.1) as*

[TABLE]

where $\Pi_{\Lambda}$ denotes the orthogonal projection onto the set $\Lambda$. This alternative scheme doesn’t affect the validity of Theorem 2.6; in fact, from the contraction property of $\Pi_{\Lambda}$ we get

[TABLE]

and this corresponds to the first step in the proof of Theorem 2.6 (see Section 3 below for more details).

An important question concerning applied settings is whether the estimator $\lambda_{t}$ defined in (2.1) performs asymptotically better than the maximum likelihood estimator $\hat{\lambda}_{t}$ calculated by optimizing the one observation log-likelihood $\ln p(X_{t}|\lambda_{t})$. The following example will showcase that there are indeed cases when utilizing (2.1) is beneficial.

Example 2.10**.**

Referring to Example 2.9 and setting $m=d=1$ for easiness of notation, we consider a sequence of independent observations $\{X_{t}\}_{t\in\mathbb{N}}$ with

[TABLE]

We assume in addition that $\lambda\mapsto A^{\prime\prime}(\lambda)$ is continuous and we restrict the parameter space to $\Lambda=[\lambda_{m},\lambda_{M}]$ for suitable real numbers $\lambda_{m}<\lambda_{M}$. Observe that Assumptions 2.2 and 2.3 hold in this case with

[TABLE]

In Theorem 2.6 we obtained an upper bound for the asymptotic mean-square error of $\lambda_{t}$ as defined in (2.1). We now want to compare it with the mean-square error of the sufficient statistic $T(X_{t})$, which we assume to be unbiased; this means considering the quantity

[TABLE]

where the last equality follows from Theorem 1.6.2 in [Bickel and Doksum, 2001]. Therefore, our estimator $\lambda_{t}$, performs asymptotically better than $T(X_{t})$ if

[TABLE]

Here, the left hand side corresponds to right hand side in (2.4) with $d=1$ while the right hand side follows from (2.5). We want this inequality to hold for all possible values of the sequence $\{\lambda_{t}^{*}\}_{t\in\mathbb{N}}$ and this is achieved by taking the infimum of the right hand side of (3), i.e., we want

[TABLE]

A simple investigation of the previous inequality shows that the left hand side increases for small values of $\ell$ or large values of $L$; hence, there exist $\bar{\ell}$ and $\bar{L}$ such that for all $\bar{\ell}\leq\ell\leq L\leq\bar{L}$ the asymptotic mean-square error of $\lambda_{t}$ is lower than the mean-square error of the sufficient statistic $T(X_{t})$. Figures (1) and(2) provide an illustration of this fact. Finally, notice that there are cases when the sufficient statistic of the exponential family is unbiased and coincides with the one observation maximum likelihood estimator, as is the case if we choose as the parameter of interest the variance of a Gaussian.

3 Proof of the main result

Using (2.1) and expanding the squared Euclidian norm we can write

[TABLE]

where we set

[TABLE]

and

[TABLE]

To treat $\mathcal{A}_{1}$ we employ Theorem 2.1.12 from [Nesterov, 2014]; with $C_{1}:=\frac{\ell L}{\ell+L}$ and $C_{2}=\frac{1}{\ell+L}$ this gives

[TABLE]

moreover, using inequality $\|a+b\|^{2}\leq 2\|a\|^{2}+2\|b\|^{2}$ we get

[TABLE]

Combining (3) with (3.2) and (3.3) we obtain

[TABLE]

Imposing that $2\alpha(\alpha-C_{2})\geq 0$, or equivalently $\alpha\geq C_{2}$, we can utilize the Lipschitz continuity of the gradient in the second line above to get

[TABLE]

Notice that according to the definitions of $C_{1}$ and $C_{2}$ we can write

[TABLE]

therefore, setting $\varphi(\alpha,L):=\sqrt{1-2L\alpha+2L^{2}\alpha^{2}}$ inequality (3) now reads

[TABLE]

Taking the conditional expectation with respect to the sigma-algebra $\mathcal{F}_{t-1}:=\sigma(X_{1},...,X_{t-1})$ of both sides above we obtain

[TABLE]

Here, we have utilized that

• •

$\lambda_{t}$ is by construction $\mathcal{F}_{t-1}$-measurable for all $t\in\mathbb{N}$;

• •

the $X_{t}$’s are independent;

• •

the expectation of the score is zero;

• •

Remark 2.4.

We now compute the expectation of the first and last members of (3) to get

[TABLE]

which together with inequality $\sqrt{a+b}\leq\sqrt{a}+\sqrt{b}$ gives

[TABLE]

The last step involves using Assumption 2.5 in the previous estimate to obtain

[TABLE]

which upon iteration yields

[TABLE]

If $\alpha<\frac{1}{L}$, then $\varphi(\alpha,L)<1$; we can therefore take the limit as $t$ tends to infinity of both sides to get

[TABLE]

moreover, the minimum of the right hand side above is attained at $\alpha=\frac{1}{l+L}$ (in view of the constraints needed on $\alpha$ to recover inequality (3)).