Personalized Optimization with User's Feedback

TL;DR

This paper introduces an online algorithm that balances known time-varying costs and unknown user satisfaction, using Gaussian processes to learn and adapt in real-time, with applications demonstrated in vehicle platooning.

Contribution

It develops a novel online optimization method that simultaneously learns user satisfaction and tracks optimal solutions in dynamic environments.

Findings

The algorithm effectively balances performance and user satisfaction.

It achieves no-regret learning of the user's utility function.

Numerical examples demonstrate practical applicability in vehicle platooning.

Abstract

This paper develops an online algorithm to solve a time-varying optimization problem with an objective that comprises a known time-varying cost and an unknown function. This problem structure arises in a number of engineering systems and cyber-physical systems where the known function captures time-varying engineering costs, and the unknown function models user's satisfaction; in this context, the objective is to strike a balance between given performance metrics and user's satisfaction. Key challenges related to the problem at hand are related to (1) the time variability of the problem, and (2) the fact that learning of the user's utility function is performed concurrently with the execution of the online algorithm. This paper leverages Gaussian processes (GP) to learn the unknown cost function from noisy functional evaluation and build pertinent upper confidence bounds. Using the GP formalism, the paper then advocates time-varying optimization tools to design an online algorithm that exhibits tracking of the oracle-based optimal trajectory within an error ball, while learning the user's satisfaction function with no-regret. The algorithmic steps are inexact, to account for possible limited computational budgets or real-time implementation considerations. Numerical examples are illustrated based on a problem related to vehicle platooning.

Equations

$$\text{P}(t):\underset{\text{mathbold}x\in D}{\max }f(\text{mathbold}x;t)=V(\text{mathbold}x;t)+U(\text{mathbold}x),\text{ for all }t\ge 0,$$

$${\text{P}}_k:\underset{\text{mathbold}x\in D}{\max }f(\text{mathbold}x;t_k)=V(\text{mathbold}x;t_k)+U(\text{mathbold}x),$$

$${\mu }_n(\text{mathbold}x)$$

$$k_n(\text{mathbold}x,\text{mathbold}{x}^{\prime })$$

$${\sigma }_n^2(\text{mathbold}x)$$

$${\hat{U}}_n(\text{mathbold}x)={\mu }_{n-1}(\text{mathbold}x)+\sqrt{{\beta }_n}{\sigma }_{n-1}(\text{mathbold}x),$$

$$\text{mathbold}{x}_k\approx \text{arg}\underset{\text{mathbold}x\in D}{\max }\{V(\text{mathbold}x;t_k)+{\hat{U}}_n(\text{mathbold}x)\},$$

$$y_n=U(\text{mathbold}{x}_k)+{\epsilon }_n,$$

$$\text{mathbold}{x}_k=\underset{N_s\text{ compositions}}{\underbrace{\mathcal{M}\circ \cdots \circ \mathcal{M}}}\circ {\phi }_{nk}(\text{mathbold}{x}_{k-1}).$$

$${\phi }_{nk}(\text{mathbold}x):=V(\text{mathbold}x;t_k)+{\mu }_{n-1}(\text{mathbold}x)+\sqrt{{\beta }_n}{\sigma }_{n-1}(\text{mathbold}x).$$

$${\gamma }_T:=\underset{A\subset D:|A|=T}{\max }\frac{1}{2}\log |\mathbf{I}+{\sigma }^{-2}\text{mathbold}{K}_A|,$$

$$\mathbf{Pr}\left\{\underset{\text{mathbold}x\in D}{\sup }|\partial U/\partial x_j|>M\right\}\le a{\mathrm{e}}^{-(M/b{)}^2},j=1,\dots ,d.$$

$${\Delta }_k:=\underset{\text{mathbold}x\in D}{\max }|V(\text{mathbold}x;t_k)-V(\text{mathbold}x;t_{k-1})|\le \Delta .$$

$${\phi }_{nk}^{\ast }-{\phi }_{nk}(\text{mathbold}{x}_k)\le \eta ({\phi }_{nk}^{\ast }-{\phi }_{nk}(\text{mathbold}{x}_{k-1})),\eta <1.$$

$$-V(\text{mathbold}x;t_k)+V(\text{mathbold}y;t_k)\le -{\nabla }_{\text{mathbold}y}V(\text{mathbold}y;t_k)(\text{mathbold}x-\text{mathbold}y)+L/2\Vert \text{mathbold}x-\text{mathbold}y{\Vert }_2^2.$$

$$-V(\text{mathbold}x;t_k)+V(\text{mathbold}y;t_k)\le d{D}_g\Vert \text{mathbold}x-\text{mathbold}y{\Vert }_{\infty }+L/2\Vert \text{mathbold}x-\text{mathbold}y{\Vert }_2^2$$

$$D_g:=\underset{\text{mathbold}y\in D}{\max }\Vert -{\nabla }_{\text{mathbold}y}V(\text{mathbold}y;t){\Vert }_{\infty }.$$

$$\frac{1}{2}\mathcal{D}(\text{mathbold}x,c)\ge \kappa ({\phi }_{nk}^{\ast }-{\phi }_{nk}(\text{mathbold}x)),\text{with}$$

$$\mathcal{D}(\text{mathbold}x,c):=-2c\underset{\text{mathbold}y\in D}{\min }\left\{-{\nabla }_{\text{mathbold}x}{\phi }_{nk}(\text{mathbold}x{)}^{\mathsf{T}}(\text{mathbold}y-\text{mathbold}x)+\frac{c}{2}\Vert \text{mathbold}y-\text{mathbold}x{\Vert }^2\right\},$$

$$R_T:=\sum _{n=1}^T{f}_n^{\ast }-f(\text{mathbold}{x}_n;t_n),$$

$${\beta }_n=2\log (2n^2{\pi }^2/(3\delta ))+2d\log (d{n}^{2}br\sqrt{\log (4da/\delta )}),$$

$$\mathbf{Pr}\Big{\{}R_{T}\leq\sqrt{C_{1}T\beta_{T}\gamma_{T}}+C_{2}+O^{*}(1)+G_{T}\Big{\}}\geq 1-\delta,$$

$$C_2=\frac{2D_g}{b\sqrt{\log (2da/\delta )}}+\frac{2L}{2d{b}^2\log (2da/\delta )}+2,$$

$$G_T:=2\sum _{n=1}^T\sum _{z=1}^{n-1}{\eta }^z{\Delta }_{n-z+1}\le 2\Delta \eta T/(1-\eta ).$$

$${\ell }_n=\underset{\text{mathbold}x\in D}{\max }|{\hat{U}}_n(\text{mathbold}x)-{\hat{U}}_{n-1}(\text{mathbold}x)|.$$

$$\sum _{n=1}^T\sum _{z=1}^{n-1}{\eta }^z{\ell }_{n-z+1}\le \frac{1}{1-\eta }{O}^{\ast }(1)=O^{\ast }(1),$$

$$\mathbold{x}_{k}\!=\!\Pi_{D}\!\left[\mathbold{x}_{k-1}\!+\!\alpha\Big{(}\nabla_{\mathbold{x}}V(\mathbold{x}_{k-1};t_{k})+\nabla_{\mathbold{x}}\hat{U}_{n}(\mathbold{x}_{k-1})\Big{)}\!\right]$$

$${\phi }_{nk}^{\ast }-{\phi }_{nk}(\text{mathbold}{x}_k)\le (1-\alpha \kappa )({\phi }_{nk}^{\ast }-{\phi }_{nk}(\text{mathbold}{x}_{k-1})).$$

$$\text{P}(t):\underset{\text{mathbold}x\in D}{\max }-\frac{1}{2}\Vert \text{mathbold}x-\overline{\text{mathbold}x}(t){\Vert }_{\mathbf{Q}}^2+\gamma \sum _{i=1}^m{U}_i(\text{mathbold}x),$$

$${\text{UC}}_i(t):=U_i(d_i(t))/\underset{d\in D}{\max }{U}_i(d)\in [0,1].$$

$$\text{mathbold}{x}_n\approx \arg \underset{\text{mathbold}x\in D}{\max }{\phi }_{nn}(\text{mathbold}x)=V(\text{mathbold}x;t_n)+{\hat{U}}_n(\text{mathbold}x),$$

Full text

††thanks: This paper was not presented at any conferences. Corresponding author A. Simonetto. Tel. +353-87-3583843. Fax +353-1-8269903.

, , ,

1 Introduction

Optimization is ubiquitous in engineering systems and cyber-physical systems including smart homes, energy grids, and intelligent transportation systems. As an example for the latter, optimization toolboxes are advocated for intelligent traffic light management and for congestion-aware systems for highways, where a networked control system could control vehicles and form a platoon optimized for lowering fuel emission and increasing vehicle spatial density. In many applications involving (and affecting) end-users, optimization problems are formulated with the objective of striking a balance between given engineered performance metrics and user’s satisfaction. Engineered metrics may include, e.g., operational cost and efficiency; on the other hand, metrics to be addressed for the users may be related to comfort (e.g., temperature in a home or building), perceived safety (e.g., distance from preceding vehicles while driving on a freeway), or simply preferences (e.g., taking a route with the least number of traffic lights).

While engineered performance goals may be synthetized based on well-defined metrics emerging from physical models or control structures, the utility function to be optimized for the users is primarily based on synthetic models postulated for comfort and satisfaction; these synthetic functions are constructed based on generic welfare models, averaged or statistical human-perception models estimated over a sufficiently large population of users’ responses, or just simplified mathematical models that make the problem tractable. However, oftentimes, these strategies do not lead to meaningful optimization outcomes. Synthetic models of users’ utilities are difficult to obtain for the associated cost and time of human studies, the data is therefore scarce, biased, and there is a lack of comprehensive theoretical models, see e.g. [1, 2]. In addition, (1) generic welfare or averaged models may not fully capture preferences, comfort, or satisfaction of individual users, and (2) synthetic utility functions may bear no relevance to a number of individuals and users, that is personalization is key (as well-known, e.g., in medicine [3]).

This paper investigates time-varying optimization problems with an objective that comprises a known time-varying cost and an unknown utility function. The known function captures time-varying engineering costs, where time variability emerges from underlying dynamics of the systems. The second term pertains to the users, and models personalized utility functions; these functions are in general unknown, and they must be learned concurrently with the solution of the optimization problem. Learning unknown utility functions for each user is key to overcome the limitations of synthetic models, and to fully capture user satisfaction from possibly limited data.

Typically, three separate timescales are considered in this setting: (i) temporal variations of the engineering cost; (ii) learning rates for the user’s utility function; and, (iii) convergence rate of the optimization algorithm. In the paper, we compress these timescales by advocating time-varying optimization tools to design an online algorithm that tracks time-varying optimal decision variables emerging from time-varying engineering costs within an error bound, while concurrently learning the user’s satisfaction with no-regret. We model the user’s satisfaction as a Gaussian process (GP) and learn its parameters from user’s feedback [4, 5]; in particular, feedback is in the form of noisy functional evaluations of the (unknown) utility function. The user’s satisfaction is learned by implementing the approximate decisions and measuring the user’s feedback to update the user’s GP model.

Overall, the main contributions of the paper are as follows:

$\bullet$ The paper extends the works [4, 5] on GP s to handle a sum of a convex engineering cost and an unknown user utility function by considering approximate online algorithms; the term “approximate” refers to the fact that intermediate optimization sub-problems in the algorithmic steps are not assumed to be solved to convergence due to underlying complexity limits. This is key in many time-varying optimization applications and increasingly important for large-scale systems [6, 7, 8, 9, 10, 11].

$\bullet$ We devise an approximate upper confidence bound algorithm in compact continuous spaces and provide its regret analysis to solve the formulated problem. As in [4], for the time-invariant and squared exponential kernel case we recover a cumulative regret of the form $O^{*}(\sqrt{dT(\log T)^{d+1}})$, where $d$ is the dimension of the problem and $T$ is horizon111The notation $O^{*}$ means a big-$O$ result up to polylog factors.. As in [5], the time-variation of the cost yields an extra $O(T)$ term in the cumulative regret. We explicitly show how the regret depends on the convergence rate of the optimization algorithm that is used, and how it is dependent on the choice of the kernel.

$\bullet$ We further detail the regret analysis to the case where we use a projected gradient method to compute the approximate optimizers of the upper confidence bound algorithm, and we extend the analysis to vanishing cost changes.

The performance of the proposed algorithm is then demonstrated in a numerical example derived from vehicle platooning, where the inter-vehicle distances are computed in real-time based on both the engineering and user’s comfort perspectives.

Incorporating user’s satisfaction in the decision-making process is not a new concept  [12]; yet, it is starting to play a major role in emerging data-driven optimization and control paradigms, especially within the context of cyber-physical-social systems or cyber-physical-and-human systems. For example, [13] employs control-theoretic techniques to model human behavior and drive systems to suitable working conditions. In this context, human behavior has been captured as a stochastic process [14] and discrete decision models [15], among others.

User’s satisfaction and preferences have been modelled as a GP in the machine learning community; see, e.g., [16, 17]. For an account of Gaussian processes we refer to  [18], while for their use in control we refer the reader to [19, 20, 21, 22]. User’s satisfaction and comfort has been taken into account from a control perspective in, e.g., control systems for houses, electric vehicles, and routing [23, 24, 25].

In all the works mentioned above, the users’ functions were modelled based on synthetic models or they were learned a priori (i.e., before the execution of the algorithm). Preferences are also learned among a finite number of possibilities. In this paper, we incorporate learning modules in the online optimization algorithm (with the learning and the optimization tasks implemented in a concurrent timescale), we bypass assumptions on discrete preference sets, and we do not utilize synthetic models. The proposed strategy enables fine-grained personalization.

With respect to works that use GPs in control and reinforcement learning (RL)  [26, 22, 27, 28, 29], they typically use GPs to model unknown system dynamics and state transitions; in this paper, we do not consider system dynamics and state transitions, since we focus on a purely optimization strategy. Even when the reward function is modeled as a GP, the RL machinery is more complex than what we present here; RL-based methods also implicitly assume a timescale separation between optimization and learning, which is not the case in this paper.

The techniques and ideas that are expressed in this paper are related to inverse reinforcement learning [30] -[31] (even though in those settings one still needs to know what is a “good” action via demonstration), restless bandit problems [32] (even though we do not use their machinery and they are typically in discrete spaces, while we are in compact continuous spaces), and (partially) to socially-guided machine learning [33]. Ideas in this paper are also close in spirit to preference elicitation studies; see, e.g., [34, 35, 36].

Online convex optimization and learning are active research areas in the control and signal processing communities, see e.g. [37, 38, 39, 40, 41, 42, 43, 44, 45, 46]; broadly, this paper aims at expanding this line of work by cross-fertilizing time-varying optimization, online learning, and Gaussian Processes by incorporating user’s feedback, and by providing an illustrative application in vehicle control.

Organization. The remainder of the paper is organized as follows. Section 2 presents the mathematical formulation and the proposed approach using an approximate upper confidence bound algorithm. The convergence properties of said algorithm are analyzed in Section 3, along with two specifications of the algorithm in cases where we use a projected gradient method and the cost changes are vanishing. In Section 4, we report numerical examples. Proofs of all the results are in the Appendix.

2 Formulation and Approach

Consider a decision variable $\mathbold{x}\in D\subset\mathbf{R}^{N}$, and a possibly time-varying objective function $f(\mathbold{x};t):\mathbf{R}^{N}\times\mathbf{R}_{+}\to\mathbf{R}$ that has to be maximized. The function $f(\mathbold{x};t)$, where $t>0$ represents time, is given by the sum of two terms: a concave engineering welfare function $V(\mathbold{x};t):\mathbf{R}^{N}\times\mathbf{R}_{+}\to\mathbf{R}$, and a user’s satisfaction function $U(\mathbold{x}):\mathbf{R}^{N}\to\mathbf{R}$. The function $V(\mathbold{x};t)$ is known (or can be easily evaluated at time $t$) and represents a deterministic function that evaluates system rewards or engineering costs; on the other hand, $U(\mathbold{x})$ is unknown and has to be learned. Accordingly, the problem to be solved amounts to:

[TABLE]

that is, the objective is to find an optimal decision $\mathbold{x}^{*}(t)$ for each time $t$. Here, we implicitly assume that the user’s satisfaction function $U(\mathbold{x})$ does not change in time; however, slow changes of the function do not affect reasoning and technical approach. Furthermore, $U(\mathbold{x})$ is not necessarily a concave function, thus rendering $\textsf{P}(t)$ a challenging problem even in the static setting. The structure of (1) naturally suggests that optimal decisions strike a balance between design choices (which may be time-varying, e.g., in tracking problems) and user’s preferences, which are often slowly varying and not known beforehand.

As a first step, consider sampling the problem (1) at discrete time instances $t_{k}$, $k=1,2,\ldots$, with $h$ a given sampling period. This leads to a sequence of time-invariant problems as:

[TABLE]

which we want to solve approximately within the sampling period $h$ to generate a sequence of approximate optimizers $\{\mathbold{x}_{k}\}_{k\in\mathbf{N}}$ (one for each problem $\textsf{P}_{k}$) that eventually converges to an optimal decision trajectory $\{\mathbold{x}^{*}_{k}\}_{k\in\mathbf{N}}$ up to a bounded error. When only the known engineering welfare function $V(\mathbold{x};t)$ is considered, this tracking problem has been considered in various prior works; see, e.g., [47, 48, 8, 49, 50, 10, 11] and pertinent references therein. A key difference in the setting proposed in this paper is that we construct approximate optimizers concurrently with the learning of the (unknown) user’s function $U(\mathbold{x})$. The main operating principles of the algorithm to be explained shortly are, qualitatively, as follows: (i) at time $t_{k}$, an approximate optimizer $\mathbold{x}_{k}$ is computed based on a partial knowledge of $U(\mathbold{x})$; (ii) the optimizer is implemented, and it generates some “feedback” from the user in the form of, for example, $y_{k}=U(\mathbold{x}_{k})+\varepsilon$, where $\varepsilon$ is noise; and, (iii) $y_{k}$ is collected and utilized to “refine” the knowledge of $U(\mathbold{x})$. At time $t_{k+1}$, the process is then repeated.

To this end, the paper leverages a Gaussian process (GP) model for the unknown user function $U(\mathbold{x})$. Such non-parametric model is advantageous in the present setting because of (1) the simplicity of the online updates of both mean and covariance; (2) the inherent ability to handle asynchronous and intermittent updates (which is an important feature in user’s feedback systems); and, (3) the implicit and smooth handling of measurement (i.e., feedback) noise. Accordingly, let $U(\mathbold{x})$ be specified by a GP with mean function $\mu(\mathbold{x}):=\mathbf{E}[U(\mathbold{x})]$ and covariance (or kernel) function $k(\mathbold{x},\mathbold{x}^{\prime}):=\mathbf{E}[(U(\mathbold{x})-\mu(\mathbold{x}))(U(\mathbold{x}^{\prime})-\mu(\mathbold{x}^{\prime}))]$. We assume bounded variance; i.e., $k(\mathbold{x},\mathbold{x})\leq 1,\mathbold{x}\in D$. For GPs not conditioned on data, we assume without loss of generality that $\mu\equiv 0$, i.e., GP$(0,k(\mathbold{x},\mathbold{x}^{\prime}))$. In this context, the user function $U(\mathbold{x})$ is assumed to be a sample path of GP$(0,k(\mathbold{x},\mathbold{x}^{\prime}))$. Let $A_{n}=\{\mathbold{x}_{1}\in D,\ldots,\mathbold{x}_{n}\in D\}$ be a set of $n$ sample points and let $y_{i}=U(\mathbold{x}_{i})+\varepsilon_{i}$, $\varepsilon_{i}\sim N(0,\sigma^{2})$ i.i.d. Gaussian noise, be the noisy measurements at the sample points $\mathbold{x}_{i}$, $i=1,\ldots,n$. Let $\mathbold{y}_{n}=[y_{1},\ldots,y_{n}]^{\mathsf{T}}$. Then the posterior distribution of $(U(\mathbold{x})|A_{n},\mathbold{y}_{n})$ is a GP distribution with mean $\mu_{n}(\mathbold{x})$, covariance $k_{n}(\mathbold{x},\mathbold{x}^{\prime})$, and variance $\sigma_{n}^{2}(\mathbold{x})$ given by:

[TABLE]

where $\mathbold{k}_{n}(\mathbold{x}):=[k(\mathbold{x}_{1},\mathbold{x}),\ldots,k(\mathbold{x}_{n},\mathbold{x})]^{\mathsf{T}}$, and $\mathbold{K}_{n}$ is the positive definite kernel matrix $[k(\mathbold{x},\mathbold{x}^{\prime})]_{\mathbold{x},\mathbold{x}^{\prime}\in A_{n}}$.

With this in place, in order to approximately solve the sequence of optimization problems $\{\textsf{P}_{k}\}$ (where we remind that the objective function is partially unknown), we utilize an online approximate Gaussian process upper confidence bound (AGP-UCB) algorithm as described in the AGP-UCB table in the next page.

2.1 Description of AGP-UCB.

We now describe the main steps of AGP-UCB. First of all, notice that there are two indexes $k$ and $n$; $k$ represents the time/optimization counter (i.e., the time index for the execution of the step (7)), while $n$ is the data counter (see step [S3] in AGP-UCB). Time and optimization steps are fused into $k$, since one step of the algorithm is performed at each time. The data counter $n$ does not generally coincide with $k$, since feedback may not be given at each algorithmic step. After initialization, in (6), a proxy for the unknown function $U(\mathbold{x})$ is built using an upper confidence bound. Upper confidence bound methods are popular in stochastic bandit settings; see, e.g., [4, 51]. The aim of an upper confidence bound method is to trade off exploitation (areas with high mean but low variance) and exploration (areas with low mean but high variance). The upper confidence bound depends on the parameter $\beta_{n}$ which is set as described in Th. 1.

Step (7) is utilized to find $\mathbold{x}_{k}$ based on the (current) upper confidence bound and the new function $V(\cdot;t_{k})$. In [4], it is assumed that (7) can be solved to optimality; on the other hand, we consider a setting where only an approximate optimizer of (7) can be obtained at each time $t_{k}$. In particular, we consider a case where it might not be possible to execute the algorithm until convergence to an optimal solution within a period of time $h$ [47, 48, 8, 49, 50]; this is the case where, for example, the sampling period $h$ is too short (relatively to the time required to perform an algorithmic step and the number of iterations required to converge) or the problem is computationally demanding. Thus, denoting as $\mathcal{M}$ the map of a given algorithmic step (e.g., gradient descent, proximal method, etc.) and assuming that one can perform $N_{s}$ steps within an interval $h$, $\mathbold{x}_{k}$ is obtained as:

[TABLE]

where $\varphi_{nk}(\mathbold{x})$ is defined as:

[TABLE]

For example, if $\mathcal{M}$ is the map of a gradient method, then (8) implies that one runs $N_{s}$ gradient steps222 The choice of the algorithm will be dictated by Assumption 4 and it will influence the overall performance.; this example will be explained shortly in Section 3.2. See also the works [47, 48, 8, 49, 50] (and pertinent references therein) for the case of $N_{s}=1$. From a more utilitarian perspective, this setting allows one to allocate computational resources to solve (7) parsimoniously; spending resources to solve (7) to optimality might not provide performance gains since the user’s function is not known accurately. We remark here that not solving (7) to optimality and still guaranteeing convergence results is one of the key contribution s of this paper.

Finally, step [S3] involves the gathering of the user’s feedback in the form of a measurement of $U(\mathbold{x}_{k})$. This feedback is in general noisy (e.g., it may be collected by sensors or be quantized), and it could be intermittent in the sense that it might not be available at every time step (this explains why we utilize the subscript $k$ for the temporal index and $n$ for the updated of $\hat{U}_{n}(\mathbold{x})$; of course, $k=n$ if $y_{n}$ is collected at each time $k$). If the feedback is available, the GP model is updated via (3)-(4), otherwise not. This is also further explained in Figure 1, where the temporal index $k$ is dictated by the optimizer counter, while the data counter $n$ is dictated by how many feedback data we get from the users. Since $U$ is time-invariant, it does not matter when the feedback is received, as long as it is correctly labeled.

Remark 1

*AGP-UCB assumes the possibility to acquire the feedback $y_{n}$ (associated with a given decision variable $\mathbold{x}_{k}$). Again, feedback may be infrequent, and does not have to happen at every time $k$. Realistic examples of applications include, e.g., a platooning task (see Section 4) with sensors in the seat of the vehicle determining stress levels; a personalized medicine example, where the effects of a treatment can be monitored by user’s feedback via dedicated mobile app; a building temperature control, where the user’s satisfaction can be evaluated by wearables via perspiration levels. *

3 Regret Bounds

In this section, we establish cumulative regret bounds for the AGP-UCB algorithm, which depend on the amount of feedback we receive. A critical quantity affecting these bounds is the maximum information gain $\gamma_{T}$ after $T$ feedback rounds, which is defined as

[TABLE]

with $\mathbold{K}_{A}=[k(\mathbold{x},\mathbold{x}^{\prime})]_{\mathbold{x},\mathbold{x}^{\prime}\in A}$, and where $A$ represents the set of $T$ points $\mathbold{x}\in D$, which maximizes the expression [cf. [4]].

The regret bounds are of the form $O^{*}(\sqrt{T\beta_{T}\gamma_{T}})+O(T)$, where the first term is sub-linear and depends on how fast one can learn the unknown user’s profile, while the second term is due to the time drift (i.e., temporal variability) of the design function. These regrets bound resemble [4] for the time invariant case and [5] for the time-varying one333We leave for future research the use of more sophisticated methods that learn, predict, or bound the variations of the cost function and might deliver no-regret results also in the time-varying setting, see [52, 53]. . The main proof techniques that we employ are a combination of GP regret results, convex analysis, and Kernel ridge regression theory.

3.1 Main result

The following assumptions are imposed.

Assumption 1

The function $-V(\mathbold{x};t)$ is convex and $L$-strongly smooth over $D$, uniformly in $t$.

Assumption 2

Let $D\subseteq[0,r]^{d}$ be compact and convex, $d\in\mathbf{N}$, $r>0$. Let the kernel $k(\mathbold{x},\mathbold{x}^{\prime})$ satisfy the following high probability bound on the derivatives of the GP sample paths $U$: for some constants $a,b>0$:

[TABLE]

Assumption 3

The changes in time of function $V(\mathbold{x};t)$ are bounded, in the sense that at two subsequent sampling times $t_{k}$ and $t_{k-1}$, the following bound holds

[TABLE]

for a given $\Delta<\infty$.

Assumption 4

Recall the definition of $\varphi_{nk}(\mathbold{x})$ in (9) and let its maximum be $\varphi_{nk}^{*}$. There exists an algorithm $\mathcal{M}$ with linear convergence that, when applied to $\varphi_{nk}(\mathbold{x}_{k-1})$, yields a point $\mathbold{x}_{k}$ for which

[TABLE]

Assumption 5

The kernel $k(\cdot)$ is not degenerate; the mean and variance are well-behaved, so that can be represented by generalized Fourier series of the kernel basis.

Most of the assumptions are standard in either time-varying optimization or Gaussian process analysis; see, e.g., [4, 47]. From Assumption 1, one has that the following holds for any $\mathbold{x},\mathbold{y}\in D$ and $t_{k}$:

[TABLE]

Further, since $D$ is compact, one has that:

[TABLE]

where $D_{g}$ is defined (uniformly in $t$) as:

[TABLE]

Assumption 2 is standard in the analysis of GPs (see, e.g., [4]); it holds true for four-time s differentiable stationary kernels, such as the widely-used squared exponential and (some) Matérn kernels, see [54, Eq. (2.33)], or equivalently [55, Theorem 5]. This assumption is needed when $D$ is compact in order to ensure smoothness of the GP samples (Cf. Appendix D for how to compute $a,b$).

Assumption 3 is required in time-varying settings to bound the temporal variability of the cost function; specifically, it presupposes that the decision $\mathbold{x}$ yields similar function values for $V$ at two subsequent time (with their difference being upper bounded by $\Delta$). This assumption is also common in online optimization and machine learning; see, e.g., [52, 53].

Assumption 4 constraints the possible algorithms that are employed to approximately solve (7) in [S2]. Two cases are in order based on whether the assumption holds locally around the optimal trajectory (in the paper, the term optimum refers to the global one) or globally. In the first case, results will be valid around the optimal trajectory (provided that the algorithm is started close enough). In the second case, results will be valid in a general sense. Since the conditions are the same in both cases, the paper will hereafter focus on the second case. In a convex setting, an algorithm $\mathcal{M}$ with $Q$-linear convergence can be found when the function is strongly convex and strongly smooth; when the function $-\varphi_{nk}(\mathbold{x})$ is nonconvex, results are available for special cases. A noteworthy example (which is also explored in the simulation results) is when $-\varphi_{nk}(\mathbold{x})$ has a Lipschitz-continuous gradient and it satisfies the Polyak-Lojasiewicz (PL) inequality; in this case, the (projected) gradient method has a global linear converge rate and Assumption 4 is satisfied [56]. This case implies that the function $-\varphi_{nk}(\mathbold{x})$ is invex (i.e., any stationary point must be a global minimizer).444Other cases in which Assumption 4 is verified for the general nonconvex case could be found by using recent global convergence analysis of ADMM [57, 58], and left for future research. Note that if one has a $C/k$ or $C/\sqrt{k}$ method, then one can transform it to a linear converging method by running at least $k>C/\eta$ or $k>(C/\eta)^{2}$ iterations, respectively.

With this in mind, one could substitute Assumption 4 with the following (more restrictive) assumption.

Assumption 6

The function $-\varphi_{nk}(\mathbold{x})$ in (9) has a Lipschitz continuous gradient with coefficient $\Theta$ and satisfies the Polyak-Lojasiewicz (PL) inequality as,

[TABLE]

[TABLE]

over $\mathbold{x}\in D$, for some $\kappa>0$, $c>0$, uniformly in time (i.e., $\forall\,k$).

Assumption 6 can be also required in a local sense only (around the optimal trajectory), if necessary555Condition (13) imposes that the gradient of the cost function is properly lower bounded. If $D$ is the whole space, $\mathcal{D}$ reduces to $\|\nabla\varphi_{nk}(\mathbold{x})\|^{2}$ and the condition becomes of easier interpretation. Condition (13) also implies invexity. See [56, 59] for detailed discussions. .

Assumption 5 is a mild technical assumption needed to determine the learning rates of the proposed method; see also [18, 60, 61, 62]. The non-degeneracy of the kernel holds true for squared exponential kernels [63]. If the assumption is not satisfied, then one would learn a smoother version of $U$.

Our theoretical regret analysis will be based on the amount of feedback points $n$ that we receive. From data point $n-1$ to data point $n$, a number of time/optimization steps $k$ may have passed. However, to simplify the exposition and the notation, hereafter we set the time index $k$ and data index $n$ in AGP-UCB to be the same; that is, $k=n$ ( and we receive feedback at every time $k$/for every approximate decision variable $\mathbold{x}_{k}$). This is done without loss of generality, since $k=n$ represents the worst case for the regret analysis ( from the data standpoint); in fact, for $k=n$ the time scales of the variation of $\hat{U}_{n}$ and the convergence of the optimization algorithm are the same, which causes coupling of the two. For $k\gg n$, the convergence of the optimization algorithm is basically achieved for any $n$ and we are back to timescale separations as in [4]. We finally notice that if we have, say, at least $p>0$ optimization/time steps within each data point and at most $q\geq p$, the results here presented will be valid substituting $\eta$ with $\eta^{p}$ and $\Delta$ with $\Delta(1-\eta^{q})/(1-\eta)$.

Let then $k=n$, without loss of generality.

Define the cumulative regret $R_{T}$ as:

[TABLE]

where $f^{*}_{n}$ is the maximum of function $f(\cdot;t_{n})$ at time $t_{n}$ provided by an oracle666 The definition of regret here is compatible with the one used in [4, 5] in the static and time-varying case. However, $R_{T}$ could be slightly better interpreted as dynamic cumulative regret against a time-varying comparator [53], or cumulative objective tracking error. . The following result holds for the regret.

Theorem 1** **(Regret bound)

Let Assumptions 1–5 hold. Pick $\delta\in(0,1)$ and set the parameter $\beta_{n}$ as

[TABLE]

for $n\geq 1$, where $a,b,d$, and $r$ are the parameters defined in Assumption 2. Running AGP-UCB with $\beta_{n}$ for a sample $U$ of a GP with mean function zero and covariance function $k(\mathbold{x},\mathbold{x}^{\prime})$, and considering $k=n$, we obtain a regret bound of $O^{*}(\sqrt{dT\gamma_{T}})+O(T)$ with high probability. In particular,

[TABLE]

where $C_{1}=8/\log(1+\sigma^{-1})$,

[TABLE]

and,

[TABLE]

$\diamond$

Proof. See Appendix A.

Theorem 1 asserts that the average regret $R_{T}/T$ converges to an error bound with high probability. It can be noticed that the error bound is a function of the variability of the function $V(\mathbold{x};t)$; if $V(\mathbold{x};t)$ is time-invariant, then a no-regret result can be recovered. The term $\sqrt{C_{1}T\beta_{T}\gamma_{T}}+2$ in the bound can be found in the result of [4, Theorem 2] too; in particular, the term $C_{1}$ depends only on the variance of the measurement noise. On the other hand, the term $C_{2}-2$ is due to the time-varying concave function $V(\mathbold{x};t)$, and it explicitly shows the linear dependence on the Lipschitz constant $L$ of $V(\mathbold{x};t)$ as well as on the bound $D_{g}$ on the gradient of $V(\mathbold{x};t)$. The term $G_{T}$ is also due to the time variation of $V(\mathbold{x};t)$.

The term $O^{*}(1)$ emerges from the approximation error in the step [S2] of the algorithm and it is key in our analysis; that is, because only a limited number of algorithmic steps $\mathcal{M}$ are performed, instead of running the optimization algorithm to convergence. It is also due to the fact that the function $\hat{U}_{n}$ changes every time a new measurement is collected. Define the learning rate error $\ell_{n}$ as

[TABLE]

Then, the error term $O^{*}(1)$ comes from the sum

[TABLE]

whose proof is deferred in the Appendix. From (18), we can see that if the step [S2] is carried out exactly, $\eta\to 0$, then this term vanishes. On the other hand, if $\eta>0$, then the error is weighted by the changes in the surrogate function $\hat{U}_{n}(\mathbold{x})$. Because of the Bayesian update, $\ell_{n}$ converges fast enough so that its cumulative error is constant.

The following result pertains to squared exponential kernels.

Theorem 2** **(Squared exponential kernels)

Under the same assumptions of Theorem 1, consider a squared exponential kernel. Then $\gamma_{T}=O((\log T)^{d+1})$ and the cumulative regret is of the order of $O^{*}(\sqrt{dT(\log T)^{d+1}})+O(T)$. $\diamond$

Proof. See Appendix A.

Theorem 2 is a customized version of Theorem 1 for a squared exponential kernel. Other special cases can be derived for other kernels as in [4], but are omitted here due to space constraints. In the following, we exemplify the results of Theorem 2 for the case where the algorithmic map $\mathcal{M}$ represents a projected gradient method.

3.2 Example: AGP-UCB with projected gradient method

Consider a projected gradient method, applied to a time-varying problem where Assumption 6 holds. Consider further the case where only one step of the projected gradient method can be performed in [S2] (i.e., $N_{s}=1$); then, denoting as $\Pi_{D}[\mathbf{y}]:=\arg\min_{\mathbold{x}\in D}\|\mathbf{y}-\mathbold{x}\|_{2}$ the projection operator, in this case [S2] is replaced with:

[S2’] Update $\mathbold{x}_{k}$ as:

[TABLE]

and implement $\mathbold{x}_{k}$.

Step [S2’] encodes a projected gradient method on the function $\varphi_{nk}(\mathbold{x})$ with stepsize $\alpha$. Note that for the Bayesian framework, the derivative $\nabla_{\mathbold{x}}\hat{U}_{n}(\mathbold{x}_{k-1})$ is straightforward to compute and no approximation has to be made; this is in contrast with online bandit methods where the gradient is estimated from functional evaluations (see e.g., [64, 65, 66] and references therein); this is one of the strengths of Bayesian modelling, which however comes at the cost of an increase in computational complexity due to the GP updates.

Under Assumption 6 and with the choice of stepsize $\alpha\leq 1/\Theta$, a suitable extension of the results in [56] (considering the proximal-gradient method with $g$ being the indicator function, see Appendix B) yields the following convergence result for the iteration (19):

[TABLE]

In this context, the following corollary is therefore in place.

Corollary 1

[AGP-UCB with projected gradient method]* Consider the modified Step [S2’] as in (19), with $\alpha\leq 1/\Theta$. Under the same assumptions of Theorem 1, but with Assumption 4 replaced with Assumption 6 with $c=1/\alpha$, Theorem 1 holds with the specific value $\eta=1-\alpha\kappa$. $\diamond$*

Proof. See Appendix B.

3.3 Vanishing Changes

Consider the case where $\Delta_{k}$ in Assumption 3 vanishes in time; that is, $\Delta_{k}\to 0$ as $k\to\infty$. This case is important when the variations in the engineering cost function eventually vanish (for example, if the cost function is derived by a stationary process that is learned while the algorithm is run, as typically done in online convex optimization [67]).

Theorem 3** **(Vanishing changes)

Under the same assumptions of Theorem 1, if $\Delta_{k}\to 0$ as $k\to\infty$, then $G_{T}$ in Theorem 1 can be upper bounded by a sublinear function in $T$ and we obtain a no-regret result.

Furthermore, if $\Delta_{k}$ decays at least as $O(1/\sqrt{k})$, then the result of Theorem 1 on the regret $R_{T}$ is indistinguishable from the static result of [4, Theorem 2] in a $O^{*}$ sense. $\diamond$

Proof. See Appendix C.

4 Numerical evaluation

This section considers an example of application of the proposed framework in a vehicle platooning problem, whereby vehicles are controlled to maintain a given reference distance between each other and follow a leader vehicle, inspired by [68, 69, 70], and it provides illustrative numerical results. This problem includes all the modeling elements discussed in the paper, and its real-time implementation requirements are aligned with the design principles of the proposed framework.

Consider then $m+1$ automated vehicles that are grouped in a platoon. The leading vehicle is labeled as [math], while the vehicles following the leading one are indexed with increasing numbers from $1$ to $m$. The platoon leading and desired velocity is $v(t)$, while the inter-vehicle distances are denoted as $d_{i}(t)$ for $i=1,\ldots,m$ (cf. Fig. 2). Consider the problem of deciding which are the best inter-vehicle distances such that they are as close as possible to some desired values that are dictated by road, aerodynamics considerations, and platoon’s speed, while being comfortable to the car riders; e.g. the automated vehicles distances are not too different than the distances that users would naturally adopt in human-driven vehicles.

Let $\mathbold{x}:=[d_{1},\ldots,d_{m}]^{\mathsf{T}}$, and denote as $\bar{\mathbold{x}}$ the time-varying vector of distances that one would obtain by considering only the engineering cost. Then, the problem considered in this section is of the form:

[TABLE]

where $D$ is a compact set representing allowed distances between vehicles; $U_{i}$ is the function capturing the “comfort” of user $i$; and, $\|\mathbold{y}\|^{2}_{\mathbf{Q}}:=\mathbold{y}^{\mathsf{T}}\mathbf{Q}\mathbold{y}$ is the weighted norm based on the positive definite matrix $\mathbf{Q}$. We further consider the case $\gamma=1$, and $\mathbf{Q}$ being not diagonal (this way, the decision variables are coupled). We also set $U_{i}(\mathbold{x})=U_{i}(d_{i})$; that is, the comfort of the $i$-th user depends only on the distance with the preceding vehicle $i-1$ (however, more general models can be adopted). The number of vehicles following the leading one is set to $m=2$.

Functions $U_{i}(d_{i})$ are sample paths from a GP$(0,k(\mathbold{x},\mathbold{x}^{\prime}))$ with a squared exponential kernel with length scale parameter $\ell=1$, i.e., $k(\mathbold{x},\mathbold{x}^{\prime})=\exp\left(-\frac{1}{2\ell^{2}}\|\mathbold{x}-\mathbold{x}^{\prime}\|^{2}\right)$, selected to mimic log-normal functions, and specifically, functions of the form $\exp(-\log(d_{i})^{2}/\xi_{i}^{2})/{(\xi_{i}d_{i})}$, with $\xi_{1}=0.6$ and $\xi_{2}=0.7$; the maximizer of $\sum_{i=1}^{m}U_{i}(\mathbold{x})$ does not coincide with the one of the engineering function to avoid a trivial solution. Further, the functions $U_{i}$’s are different for each vehicle. Using log-normal functions for users’ comfort was motivated, see e.g. [71], by observing that inter-vehicle times and distances follow log-normal distributions. Intuitively, smaller distances are more critical than larger ones, and there is a given distance after which the comfort decreases since the users feel that they are too slow relative to the preceding vehicle. In the simulations, we learn functions $U_{i}$ via a GP with the same kernel.

It is important to notice that, with the selected parameters, the approximate function $-\varphi_{nk}(d_{i})$ is invex, it has Lipschitz continuous gradient, and it also satisfies the Polyak-Lojasiewicz inequality over $D$. With this in mind, we can readily apply a projected gradient method with linear convergence. Here, we apply one step of projected gradient method per time step $k$.

For the numerical tests, we set the probability $\delta$ to $\delta=0.1$, and the step size for the gradient method to approximately solve (7) is set to $\alpha=0.1$. The set $D$ is $D=[0,1]^{m}$ (distances are properly scaled so that the set $[0,1]$ maps to a real distance of $[0,3]$). We compute the constants $a,b$ of Assumption 2 as described in Appendix D, and set $a=1.1,b=2$, so $\beta_{n}$ can be determined.

The desired distances $\bar{\mathbold{x}}(t)$ are set to $\bar{d}_{i}=.33+.25\sin(\pi\omega t)$ for all $i$’s, where $\omega$ is a tunable parameter. The rationale for modeling the $\bar{\mathbold{x}}(t)$ in this way is to capture (1) a stationary case ($\omega=0$) and a dynamic case where the distances change because of dynamic traffic conditions (and therefore varying $v(t)$), road changes, etc. Users’ feedback comes as a noisy sample of their comfort function, and the noise is modeled as a zero-mean Gaussian variable with variance $\sigma=0.1$. Feedback in this example can come in different ways: it can come at low frequency, if the users hit the break or the accelerator every time they feel too close or too far from the vehicle in front, or it can come at higher frequency, if the users are equipped with heart rate/breathing rate sensors (which can be in smartwatches or incorporated in the seat of the vehicles [72]) which may be used as proxies of stress and discomfort. Here, we assume that the user’s feedback is either received within the sampling period (in which case the optimization counter $k$ is the same as the data counter $n$, i.e., $n=k$) or that is received only every $4$ optimization steps (in which case $k=4n$).

In Figure 3, we show the performance of the AGP-UCB algorithm varying $\omega$. On the vertical axis, we plot the regret $\frac{1}{T}R_{T}$, averaged on 25 different runs of the algorithm. As expected, we observe that when $\omega=0$, we obtain a no-regret scenario, with $\frac{1}{T}R_{T}$ eventually going to zero; in this case, it goes to zero as $O^{*}(\sqrt{T(\log T)^{2}}/T)$ (note that the learning dimension is $1$ in our example, and in this case $n=k$). In the time-varying scenarios instead, there is an asymptotic error bound, thus corroborating the theoretical results. We also notice that, in the case of $k=4n$, the algorithm is slower to converge, since the number of feedback points is smaller.

In Figure 4, we report how the distance to the platoon leader changes over time for the case $n=k$, $\omega=0.4$, up to $k=400$. We notice how the engineering best (Eng. best) for the different following vehicles has larger excursions, while the computed approximate decision variable – taking into account user’s comfort – is smoother. Roughly speaking, for large distances users feel they are losing the vehicle in front, for smaller distances they feel too close.

In Figure 5, we illustrate the learning of user $2$ satisfaction function as more and more feedback comes in, for the case $n=k$, $\omega=0.4$. In red, we represents the mean, while in green, we represents the 1-$\sigma$ bounds. The grey blobs are the feedback data points.

We notice that in this example, the vehicles are supposed to follow the new set-points (i.e., the new inter-vehicle distances) rapidly i.e., faster than the sampling period $h$. If this is not the case, the analysis of the algorithm could be extended to include an actuation error; we plan to investigate this in future research.

4.1 Comparison with other approaches

We now compare the proposed algorithm to other possible approaches. In particular, we discuss (i) a projected gradient method based on a synthetic one-fits-all model of the user’s satisfaction function; (ii) a projected zeroth-order method based on the user’s feedback at the current and previous point(s) to estimate the gradient of the user’s satisfaction function.

The synthetic model for the user’s satisfaction functions is taken to be $U_{i}(d_{i})=\exp(-\log(d_{i})^{2}/\xi_{i}^{2})/{(\xi_{i}d_{i})}$, with $\xi_{i}=0.9$ for both vehicles (to exactly match the parametric model which the GP samples mimic, with a small parameter inaccuracy). We remark that it is in general complicated to obtain accurate one-fits-all models for human preferences, and even small inaccuracies can lead to a bias and a loss of perceived fairness.

The zeroth-order method utilizes the user’s feedback to obtain an estimate of the gradient of the user’s satisfaction function at each iteration. We implement both a two-point estimate, where we use the current feedback and the previous one, as well as a four-point estimate, where we use the current feedback and three previous ones. We remark that the implementation of zeroth-order methods, which is in line with [73, 74, 75], may be in general impractical; in fact, the user may provide feedback infrequently (whereas the zeroth-order method requires feedback at each iteration). While our GP-based method does work seamlessly with intermittent noisy feedback, a zeroth-order method would not.

In Figure 6, we see how our AGP-UCB algorithm outperforms the others in terms of average regret. This is reasonable, since the synthetic model has no feedback to learn the true user’s satisfaction functions, and the zeroth-order model estimates a noisy gradient (and, likewise, it does not learn the user’s functions).

We report in Figure 7 a metric for the satisfaction of each user. We define this metric as

[TABLE]

The metric $\textsf{UC}_{i}(t)$, for each user for each time $t$, “measures” how satisfied is the user with the current decision variable $d_{i}(t)$. The results are shown for one run, and they are rolling averages of 5 (for the AGP-UCB algorithm and the synthetic functions) and of 400 (for the zeroth-order algorithms) time points to smooth the signal and capture the overall trend. It can be seen that, not only our AGP-UCB algorithm achieves the best $\textsf{UC}_{i}(t)$, but it also achieves a similar one for both vehicles. The other competing algorithms create a bias towards one of the vehicles, degrading the perceived fairness and thereby reducing the possible acceptance of the technology [76]. Finally, the zeroth-order methods appear noisy for this example and may generate frustration in a human user.

5 Conclusions

We developed an online algorithm to solve an optimization problem with a cost function comprising a time-varying engineering function and a user’s utility function. The algorithm is an approximate upper confidence bound algorithm, and it provably generates a sequence of optimizers that is within a ball of the optimal trajectory, while learning the user’s satisfaction function with no-regret. We have illustrated the result with a numerical example derived from vehicle platooning, which offered some additional inspiration for future research. We believe that our algorithm may offer a wealth of opportunities to devise more ways to integrate users’ considerations within the decision making process.

Appendix A Proof of Theorems 1-2

A.1 Preliminaries

We start by defining some quantities that will be subsequently used in the proofs. We rewrite the main algorithmic step (7) as

[TABLE]

where we recall that $k=n$; we simplify the notation from $\varphi_{nn}(\mathbold{x})$ to $\varphi_{n}(\mathbold{x})$, since there is no confusion when $k=n$.

Define the error $e_{n}$ as the difference between the optimum of the function $\varphi_{n}(\mathbold{x})$ and the approximate solution computed via (22), as

[TABLE]

where we use the superscript ${\dagger}$ for quantities that are related to the optimizer/optimum of $\varphi_{n}(\mathbold{x})$, so that the optimum of $\varphi_{n}(\mathbold{x})$ is $\varphi_{n}^{\dagger}=\varphi_{n}(\mathbold{x}_{n}^{\dagger})$.

Define the instantaneous regret $r_{n}$ as the difference between the optimum of function $f(\mathbold{x};t_{n})=V(\mathbold{x};t_{n})+U(\mathbold{x})$, which denoted as $f^{*}_{n}$, and the approximate optimum $f(\mathbold{x}_{n};t_{n})$ where $\mathbold{x}_{n}$ is computed via (22), as

[TABLE]

Further, define the instantaneous optimal regret $r_{n}^{\dagger}$ as the difference between $f^{*}_{n}$ and the optimum $f(\mathbold{x}_{n}^{\dagger};t_{n})$ where $\mathbold{x}_{n}^{\dagger}$ is the optimizer of $\varphi_{n}(\mathbold{x})$, as

[TABLE]

Lemma 5.5 in [4] is applicable here: pick a $\delta\in(0,1)$ and set $\beta_{n}=2\log(\omega_{n}/\delta)$, where $\omega_{n}>0$, $\sum_{n\geq 1}\omega_{n}^{-1}=1$. Then, for all $n\geq 1$,

[TABLE]

hold with probability $\geq 1-\delta$.

Similarly, Lemma 5.6 and Lemma 5.7 of [4] hold. For Lemma 5.7, we choose a regular discretization $D_{n}\subset D$, with $D_{n}\subset[0,r]^{d}$ and the discretization size $(\tau_{n})^{d}$. Similar to [4] we have that

[TABLE]

where $[\mathbold{x}]_{n}$ is the closest point in $D_{n}$ to $\mathbold{x}\in D$. We are now able to bound the instantaneous optimal regret $r_{n}^{\dagger}$.

Lemma 1** **(Extension of Lemma 5.8 of [4])

Let Assumptions 1 and 2 hold true. Pick a $\delta\in(0,1)$ and define a sequence $\{\omega_{n}\}_{n\in\mathbf{N}}$ such that $\sum_{n\geq 1}\omega^{-1}_{n}>0$, $\omega_{n}>0$. Set the parameter $\beta_{n}$ as

[TABLE]

where $a,b,d,r$ are the parameters defined in Assumption 2. Define the parameters $A_{0}$ and $A_{1}$ as

[TABLE]

Then instantaneous optimal regret $r_{n}^{\dagger}$ is bounded for all $n\geq 1$ as follows:

[TABLE]

with probability $\geq 1-\delta$. $\diamond$

Proof.

The proof is a modification of the one of [4, Lemma 5.8]; also in this case we use $\delta/2$ in both Lemmas 5.5 and 5.7 of [4] (which are valid here for discussion above). We report only the parts that are different. Let $\mathbold{x}^{*}_{n}$ be the optimizer of $f(\mathbold{x};t_{n})$, i.e., $f^{*}_{n}=f(\mathbold{x}^{*}_{n};t_{n})$.

By definition of $\mathbold{x}^{\dagger}_{n}$ and optimality, we have that $\varphi_{n}(\mathbold{x}^{\dagger}_{n})\geq\varphi_{n}([\mathbold{x}^{*}_{n}]_{n})$. By [4, Lemma 5.7], we have that $\hat{U}_{n}([\mathbold{x}^{*}_{n}]_{n})+1/n^{2}\geq U(\mathbold{x}^{*}_{n})$. These two combined yield

[TABLE]

Now, by convexity of $-V(\mathbold{x};t_{n})$ and the Lipschitz condition on its gradient (Assumption 1 and Eq. (11)),

[TABLE]

Therefore by substituting $\mathbold{x}=[\mathbold{x}^{*}_{n}]_{n}$ and $\mathbold{y}=\mathbold{x}^{*}_{n}$,

[TABLE]

where, by the definition of $\tau_{n}$ in [4, Lemma 5.7] and Eq.s (29)-(30), we derive the expressions of $A_{0}$ and $A_{1}$ reported in Lemma 1.

Therefore by (31) and (33),

[TABLE]

Putting the above results together, and by using (27),

[TABLE]

from which the claim. $\blacksquare$

Lemma 2

Under the same assumptions and definitions of Lemma 1, the instantaneous regret $r_{n}$ is bounded for all $n\geq 1$ as follows:

[TABLE]

with probability $\geq 1-\delta$. $\diamond$

Proof.

By definition of $e_{n}$, by calling $C=A_{0}/n^{2}+A_{1}/n^{4}$, and by (35), we have

[TABLE]

where we have use (26) in the last inequality. $\blacksquare$

We now move to the analysis of the error sequence $\{e_{n}\}_{n\geq 1}$.

Lemma 3

Let Assumptions 1 until 4 hold true. Define the learning rate error $\ell_{n}$ as

[TABLE]

Then, the error sequence $\{e_{n}\}_{n\geq 1}$ is upper bounded as

[TABLE]

$\diamond$

Proof.

From the definition of $e_{n}$ we obtain,

[TABLE]

where the inequality is due to Assumption 4 on the convergence rate of method $\mathcal{M}$. In addition,

[TABLE]

By algebraic manipulations, the error $(A)$ can be written as,

[TABLE]

The first term in the right hand side can be written as

[TABLE]

The second term is $-(\varphi_{n-1}^{\dagger}-\varphi_{n-1}(\mathbold{x}_{n}^{\dagger}))\leq 0$, by optimality. Therefore, $(A)\leq\Delta_{n}+\ell_{n}$. Furthermore, by Assumption 3,

[TABLE]

Therefore, a bound on $e_{n}$ can be derived as

[TABLE]

From which the claim follows. $\blacksquare$

In the next lemma, we characterize how the learning rate error $\ell_{n}$ evolves.

Lemma 4** **(Cumulative learning rates)

*With the same assumptions of Lemma 3 and additionally Assumption 5, we have that the cumulative learning rate error $\sum_{n=2}^{T}\ell_{n}$ is bounded as $\sum_{n=2}^{T}\ell_{n}\leq O^{*}(1).$

$\diamond$*

Proof.

First we look at the mean $\mu_{n}(\mathbold{x})$ and derive a bound on the difference between $\mu_{n}(\mathbold{x})$ and $\mu_{n-1}(\mathbold{x})$. We proceed in similar fashion for the variance $\sigma_{n}^{2}(\mathbold{x})$. Finally, we put the results together and we prove the lemma. We use the theory of Kernel ridge regression (KRR) estimation to map mean and variance as optimizers of carefully constructed optimization problems over Hilbert spaces, see [18, 62].

Let the covariance $k(\mathbold{x},\mathbold{x}^{\prime})$ be the positive definite kernel associated with the Gaussian process, let $\{\phi_{i}(\mathbold{x})\}$ be the eigenfunctions of said kernel $k$ and $\{\lambda_{i}\}$ be the associated eigenvectors. The posterior mean for GP regression (3) can be obtained as the function which minimizes an appropriately defined functional defined over Hilbert spaces:

[TABLE]

where the norm $\|\cdot\|_{\mathcal{H}}$ is the reproducing kernel Hilbert space (RKHS) norm associated with kernel $k$. Under Assumption 5 the eigenfunctions form a complete orthonormal basis and the true mean can be written as $\mu_{\infty}(\mathbold{x})=\sum_{i=1}^{\infty}\eta_{i}\phi_{i}(\mathbold{x})$, and $\mu_{n}(\mathbold{x})=\sum_{i=1}^{\infty}\zeta_{i}\phi_{i}(\mathbold{x})$. By rewriting (37) in terms of the coefficients $\zeta_{i}$ (see Eq. (7.5) in [18], with our notation), we arrive at the optimal solution:

[TABLE]

The same holds for $\mu_{n-1}(\mathbold{x})$ for some coefficients $\zeta_{i}^{\prime}$. This implies

[TABLE]

Therefore, since $\lambda_{i}>0$, $|\mu_{n}(\mathbold{x})-\mu_{n-1}(\mathbold{x})|=O(1/n^{2})$ and by looking at the cumulative difference,

[TABLE]

We move on to the variance $\sigma^{2}_{n}(\mathbold{x})$. Here we make use of the fact that the variance can be related to the bias of a noise-free KRR estimator (see [62]). In particular, define $\tilde{k}_{n}^{\mathbold{x}^{\prime}}(\mathbold{x})$ as $\mathbold{k}_{n}(\mathbold{x})^{\mathsf{T}}(\mathbold{K}_{n}+\sigma^{2}\mathbold{I}_{n})\mathbold{k}_{n}(\mathbold{x}^{\prime})$. From (4), the convergence rate of the variance $\sigma^{2}_{n}(\mathbold{x})$ is the convergence rate of $\tilde{k}_{n}^{\mathbold{x}^{\prime}}(\mathbold{x})$. Fix $\mathbold{x}^{\prime}$, from the form of $\tilde{k}_{n}^{\mathbold{x}}(\mathbold{x}^{\prime})$, one can interpret $k_{n}(\mathbold{x}_{i},\mathbold{x}^{\prime})$ as noise-free measurements and proceed as in the case of the mean, deriving

[TABLE]

Now, by rewriting the problem in terms of the coefficients of the eigenfunction expansion (so that $\tilde{k}^{\mathbold{x}^{\prime}}_{\infty}(\mathbold{x})=\sum_{i=1}^{\infty}\xi_{i}\phi_{i}(\mathbold{x})$, and $\tilde{k}_{n}^{\mathbold{x}^{\prime}}(\mathbold{x})=\sum_{i=1}^{\infty}\theta_{i}\phi_{i}(\mathbold{x})$), then

[TABLE]

With this in place,

[TABLE]

The numerator of (42) can be treated as done for the mean and is $O(1/n^{2})$, while the denominator is $O(1)$ as $n\to\infty$, therefore

[TABLE]

By combining (39) and (43), one obtains

[TABLE]

from which the claim follows. $\blacksquare$

A.2 Proof of Theorem 1

Proof.

By Lemma 2, we have that with probability greater than $1-\delta$

[TABLE]

whereas $e_{n}$ can be bounded by Lemma 3 and Lemma 4. The bound on the variance, as derived in Lemma 5.4 of [4] is valid here as well, and in particular with probability greater than $1-\delta$,

[TABLE]

with $C_{1}=8/\log(1+\sigma^{-1})$, so that by Cauchy-Schwarz:

[TABLE]

Therefore, for all $T\geq 1$

[TABLE]

where we have used: $\sum 1/n^{2}=\pi^{2}/6$ and $\sum 1/n^{4}={\pi^{4}}/{90}$. As in [4], we use the crude bound ${\pi^{4}}/{90}<\pi^{2}/6<2$, and define

[TABLE]

so that,

[TABLE]

As for the error $e_{n}$ term, by Lemma 3 and Lemma 4

[TABLE]

We bound the right-hand terms as follows

[TABLE]

and

[TABLE]

since for Lemma 4, $\sum_{n=2}^{T}\ell_{n}=O^{*}(1)$; and the claim follows.

Note that the choice of $\beta_{n}$ in the statement is the same choice of Lemma 1 (and subsequent ones), with the special selection of $\omega_{n}=n^{2}\pi^{2}/6$. Note that the error term $e_{1}/(1-\eta)$ is $O(1)$ and it can be incorporated into the $O^{*}(1)$ term. $\blacksquare$

A.3 Proof of Theorem 2

Proof.

We use  [60], or equivalently [4, Theorem 5] to bound the information gain $\gamma_{T}$ for squared exponential kernels as $\gamma_{T}=O((\log T)^{d+1})$, from which the claim is derived. $\blacksquare$

Appendix B Proof of Eq. (20)

Proof.

For completeness and because it is not straightforward, we report here the steps to extend the results of [56] to Eq. (20). First, we notice that in [56], unconstrained gradient methods and proximal-gradient methods are considered. Here we are in the proximal-case, with $g(\mathbold{x})$ representing the indicator function of the compact set $D$. For short-hand notation, we let $\psi_{nk}=-\varphi_{nk}$ and $\nabla=\nabla_{\mathbold{x}}$. In this context, iteration (19) can be equivalently written as

[TABLE]

where $g(\cdot)$ is the indicator function for $D$. Our PL inequality (13), valid for $\mathbold{x}\in D$, is precisely the same as in [56, Eq. (12)-(13)] when $g(\cdot)$ is the indicator function; in fact, the condition in [56] pertains to a $\mathcal{D}(\mathbold{x},c)$ defined as

[TABLE]

Now, by using the Lipschitz continuity of the gradient of $\psi_{nk}$, we have that

[TABLE]

In particular, we have used the following line of reasoning. First line: $g(\mathbold{x}_{k})=0$, since $\mathbold{x}_{k}\in D$; second line: Lipschitz property; third line: upper bound true for any $\alpha\leq 1/\Theta$; fourth line: definition of $\mathcal{D}$; fifth line: PL property (13).

By adding and subtracting $\psi_{nk}^{*}$ to the last inequality and rearranging, the Eq. (20) is proven. $\blacksquare$

Appendix C Proof of Theorem 3

Proof.

The error term $G_{T}$ can be written as

[TABLE]

If $\Delta_{k}\to 0$ as $k\to\infty$, then $\sum_{z=1}^{T-n}\Delta_{z+1}$ is sublinear in $T$, i.e., $\sum_{z=1}^{T-n}\Delta_{z+1}=o(T)$. Therefore

[TABLE]

which is sublinear in $T$ and the no-regret claim is proven.

To obtain a result in terms of regret that is indistinguishable from the static case of[4, Theorem 2], $G_{T}$ needs to scale at worst as $O(\sqrt{T})$, which is now the leading term of the regret. This is the case if $\Delta_{k}\leq O(1/\sqrt{k})$, for which

[TABLE]

which leads to

[TABLE]

Since the term $C_{2}+O^{*}(1)+O(\sqrt{T})$ grows slower than $\sqrt{C_{1}T\beta_{T}\gamma_{T}}$ in $T$, it can be omitted in a $O^{*}$ analysis and the claim is proven. $\blacksquare$

Appendix D Estimating the constants $a$ and $b$.

The constants $a$ and $b$ in Assumption 2 need to be estimated for different kernels; however, the procedure is rather straightforward and empirical bounds for $a$ and $b$ are not difficult to obtain. In the following, we provide an example of procedure for the squared exponential kernel.

Consider then the squared exponential kernel $k(\mathbold{x},\mathbold{x}^{\prime})=\exp\left(-\frac{1}{2\ell^{2}}\|\mathbold{x}-\mathbold{x}^{\prime}\|^{2}\right)$ and a Gaussian Process with sample paths $U\sim$GP$(0,k(\mathbold{x},\mathbold{x}^{\prime}))$. The derivative is a linear operator, therefore derivatives of GPs are still GPs, and specifically

[TABLE]

where the subscript $j$ stands for the $j$-th component. For a squared exponential kernel, then

[TABLE]

see for example [77], and in particular, $k^{\prime}_{j}(\mathbold{x},\mathbold{x})={1}/{\ell^{2}}.$

Now, take any GP, say again $U\sim$GP$(0,k(\mathbold{x},\mathbold{x}^{\prime}))$. For any given $\epsilon>0$ there exists a constant $C_{\epsilon}>0$, such that any sample path $U(\mathbold{x})$ will be bounded by the following:

[TABLE]

where $\sigma^{2}=\sup_{\mathbold{x}\in D}k(\mathbold{x},\mathbold{x})$, see [54, Eq. (2.33)], or equivalently [55, Theorem 5]. Therefore, for the derivative (if exists and bounded as in the squared exponential case),

[TABLE]

Comparing this with the Assumption 2, implies that $b=\sqrt{(2/l^{2}+2\epsilon)}$ and $a=C_{\epsilon}$.

At this point, take $b=\sqrt{(2/\ell^{2}+2\epsilon)}$. In our experiment simulations $\ell=1$ and we take $\epsilon=1$ as well, so that $b=2$.

Compute $C_{\epsilon}$ empirically, at the beginning, by drawing sample paths of the derivative and checking their suprema for different level of $M$ and their empirical frequency. In our case $C_{\epsilon}\approx 1.005$ and we set $a=1.1$ to be on the safe side. Note that theoretical bounds on $C_{\epsilon}$ exist, but they may be loose, so computing $C_{\epsilon}$ empirically is a good strategy.

As discussed, while the value for $b$ is imposed, the value of $a$ is estimated, and in general the estimated $\hat{a}\leq a=C_{\epsilon}$. We see next why this is not an issue.

Appendix E Error in the estimate of $a$

Suppose now that the estimated $\hat{a}\leq a$ (this without loss of generality). Given the bound $\hat{a}$, one can pick $\beta_{n}$ in Theorem 1, as

[TABLE]

and this is the only modification for our algorithm. For the convergence analysis, this new $\hat{\beta}_{n}$ changes our Lemma 1. In particular, $A_{0}$ and $A_{1}$ are modified. By following closely the proof of Lemma 1, and of [4, Lemma 5.7], we can define a discretization size $\hat{\tau}_{n}=dn^{2}br\sqrt{\log(2d\hat{a}/\delta)}$, and by algebraic computations, the new $A_{0}$ and $A_{1}$ with the modified $\hat{\beta}_{n}$ read,

[TABLE]

Since the new $A_{0}$ and $A_{1}$ are larger than the original ones for $\hat{a}\leq a$, then the constant term $C_{2}$ in Theorem 1 is larger than the original one. This however does not change the asymptotics of the average regret.

In addition, in Theorem 1, the term containing $\beta_{T}$ would be substituted with $\hat{\beta}_{T}$, which is smaller for the case $\hat{a}\leq a$. Once again, this does not change the asymptotics but only a multiplicative constant.

Summarizing, while having a good bound for $a$ helps to obtain a good $\beta_{n}$, and therefore a good trade-off between mean and covariance (i.e., exploitation and exploration), even rough bounds for $a$ deliver the same asymptotical behavior for the average regret.