Data assimilation and online optimization with performance guarantees

TL;DR

This paper introduces the ONDA Algorithm for real-time stochastic optimization that adaptively incorporates streaming data, ensuring high-probability performance guarantees and convergence under controlled data rates.

Contribution

The paper presents a novel online data assimilation algorithm within a distributionally robust framework that guarantees performance and convergence in streaming data scenarios.

Findings

Guarantees high-probability out-of-sample performance

Converges with sufficiently slow data streaming rate

Provides criteria for algorithm termination

Abstract

This paper considers a class of real-time stochastic optimization problems dependent on an unknown probability distribution. In the considered scenario, data is streaming frequently while trying to reach a decision. Thus, we aim to devise a procedure that incorporates samples (data) of the distribution sequentially and adjusts decisions accordingly. We approach this problem in a distributionally robust optimization framework and propose a novel Online Data Assimilation Algorithm (ONDA Algorithm) for this purpose. This algorithm guarantees out-of-sample performance of decisions with high probability, and gradually improves the quality of the decisions by incorporating the streaming data. We show that the ONDA Algorithm converges under a sufficiently slow data streaming rate, and provide a criteria for its termination after certain number of data have been collected. Simulations illustrate the results.

Equations

$$d_W({\mathbb{Q}}_1,{\mathbb{Q}}_2):=\underset{f\in \mathcal{L}}{\sup }{\int }_{\mathcal{Z}}f(\xi ){\mathbb{Q}}_1(d\xi )-{\int }_{\mathcal{Z}}f(\xi ){\mathbb{Q}}_2(d\xi ),$$

$$\underset{\mathbf{x}{\in }^d}{\inf }{\mathbb{E}}_{\mathbb{P}}[f(\mathbf{x},\xi )],$$

$${\mathbf{P}}^n({\mathbb{E}}_{\mathbb{P}}[f(\mathbf{x},\xi )]\le J_n(\mathbf{x}))\ge 1-{\beta }_n,$$

$${\mathbf{P}}^n({\mathbb{E}}_{\mathbb{P}}[f({\mathbf{x}}_n^{{ϵ}_2},\xi )]\le J_n^{{ϵ}_1}({\mathbf{x}}_n^{{ϵ}_2})+{ϵ}_1)\ge 1-{\beta }_n.$$

$$\underset{0\le x\le 1}{\min }{\mathbb{E}}_{\mathbb{P}}[-{\xi }_{1}x-{\xi }_2(1-x)].$$

$$\underset{x}{\min }{\mathbb{E}}_{\mathbb{P}}[-{\xi }_{1}x-{\xi }_2(1-x)-{\xi }^{\top }\xi ]-\rho \left(\log (x)+\log (1-x)\right),$$

$$f(\mathbf{x},\xi ):=-{\xi }_{1}x-{\xi }_2(1-x)-\rho \left(\log (x)+\log (1-x)\right)-{\xi }^{\top }\xi .$$

$$J_n^{\text{sae}}(\mathbf{x}):={\mathbb{E}}_{{\hat{\mathbb{P}}}^n}[f(\mathbf{x},\xi )]=\frac{1}{n}\sum _{k=1}^{n}f(\mathbf{x},{\xi }_k).$$

$$\mathbf{P}^{n}\{d_{W}(\mathbb{P},\hat{\mathbb{P}}^{n})\geq\epsilon\}\leq\left\{{\begin{array}[]{*{20}{l}}c_{1}e^{-c_{2}n\epsilon^{\max\{2,m\}}},&\textrm{if}\;\epsilon\leq 1,\\ c_{1}e^{-c_{2}n\epsilon^{a}},&\textrm{if}\;\epsilon>1,\end{array}}\right.$$

$$\epsilon(\beta_{n}):=\left\{{\begin{array}[]{*{20}{l}}\left(\frac{\log(c_{1}\beta_{n}^{-1})}{c_{2}n}\right)^{1/\max\{2,m\}},&\textrm{if}\;n\geq\frac{\log(c_{1}\beta_{n}^{-1})}{c_{2}},\\ \left(\frac{\log(c_{1}\beta_{n}^{-1})}{c_{2}n}\right)^{1/a},&\textrm{if}\;n<\frac{\log(c_{1}\beta_{n}^{-1})}{c_{2}},\end{array}}\right.$$

$$J_n(\mathbf{x}):=\underset{\mathbb{Q}\in {\mathcal{P}}_n}{\sup }{\mathbb{E}}_{\mathbb{Q}}[f(\mathbf{x},\xi )].$$

$$J_n(\mathbf{x}):=\underset{\mathbf{y}}{\sup }$$

$$s.t.$$

$${\mathbb{Q}}_n^{(l)}(\mathbf{x}):=\frac{1}{n}\sum _{k=1}^n{\delta }_{\{{\xi }_k-{\mathbf{y}}_k^{(l)}\}}.$$

$$J_n(\mathbf{x}):=\underset{\mathbf{u},\mathbf{v}}{\max }$$

$$s.t.$$

$$\underset{\mathbf{u},\mathbf{v}}{\max }$$

$${\mathbf{u}}_k-{\mathbf{v}}_k-{\mathbf{y}}_k^{(l-1)}⟩,$$

$$s.t.$$

$$\underset{\gamma {\in }^T}{\max }$$

$$s.t.$$

$$\varphi (n)\le (2mn){\log }_{\kappa }(\frac{{ϵ}_1}{J_n(\mathbf{x})-J_n^{(0)}(\mathbf{x})}).$$

$$\bar{\varphi }(n):=(2mn){\log }_{\bar{\kappa }}(\frac{{ϵ}_1}{J_N(\mathbf{x})-J_N^{\text{sae}}(\mathbf{x})})$$

$$\nabla h(\mathbf{y}):=\left(\begin{array}{c}x+2({\xi }_{k,1}+{\xi }_{k,2}-y_1-y_2)\\ 1-x+2({\xi }_{k,1}+{\xi }_{k,2}-y_1-y_2)\end{array}\right).$$

$$J_n^{\star }:=\underset{\mathbf{x}{\in }^d}{\inf }{J}_n(\mathbf{x}),$$

$$g_n^r(\mathbf{x}):=\frac{d}{d\mathbf{x}}{\mathbb{E}}_{{\mathbb{Q}}_n^{{ϵ}_1}({\mathbf{x}}^{(r)})}[f(\mathbf{x},\xi )].$$

$$g_n^r({\mathbf{x}}^{(r)})\in {\partial }_{ϵ}{J}_n({\mathbf{x}}^{(r)}),$$

$$J_n(\mathbf{z})\ge J_n({\mathbf{x}}^{(r)})+{g_n^r({\mathbf{x}}^{(r)})}^{\top }(\mathbf{z}-{\mathbf{x}}^{(r)})-ϵ.$$

$$g_n^r(\tilde{\mathbf{x}})\in {\partial }_{ϵ}{J}_n(\tilde{\mathbf{x}}).$$

$${\mathbf{x}}^{(r+1)}={\mathbf{x}}^{(r)}-{\alpha }^{(r)}\frac{g_n^r({\mathbf{x}}^{(r)})}{\max \{\Vert g_n^r({\mathbf{x}}^{(r)})\Vert ,1\}},$$

$$\underset{k\in \{r_n,\dots ,r+r_n\}}{\min }\{J_n({\mathbf{x}}^{(k)})\}-J_n({\mathbf{x}}_n^{\star })\le {ϵ}_2,\forall r\ge \bar{r},$$

$${\alpha }^{(i)}:=\frac{M}{\sqrt{\bar{r}+1}},\forall i\in \{r_n,\dots ,\bar{r}+r_n\},$$

$${\alpha }^{(i)}:=\frac{M}{i-r_n+1},\forall i\in \{r_n,\dots ,\bar{r}+r_n\},$$

Full text

I Introduction

Online data assimilation can benefit many applications that require real-time decision making under uncertainty, such as optimal target tracking, sequential planning problems, and robust quality control. In these problems, uncertainty is often represented by a multivariate random variable that has an unknown distribution. To quantify uncertainty and make reliable decisions, one often needs to gather a large number of samples in advance. Such requirement, however, is hard to achieve under scenarios where acquiring samples is expensive, or when real-time decisions must be made. Alternatively, recent methods such as distributionally robust optimization (DRO) have attracted recent attention due to their capability to provide out-of-sample performance guarantees with a finite number of samples. However, when the data is collected over time, it remains unclear what the best procedure is to assimilate the data in the ongoing optimization process. Motivated by this, this work studies how to incorporate finitely streaming data into a DRO problem, while guaranteeing out-of-sample performance via the generation of time-varying certificates.

Literature Review. Optimization under uncertainty is a vast research area, and as such, available methods include stochastic optimization [1] and robust optimization [2]. Recently, data-driven distributionally robust optimization has regained popularity thanks to its out-of-sample performance guarantees, see e.g. [3, 4] and [5, 6], for a distributed algorithm counterpart, and references therein. In this setup, one defines a set of distributions or ambiguity set, which contains the true distribution of the data-generating system with high probability. Then, the out-of-sample performance of the data-driven decision is obtained as the worst-case optimization over the ambiguity set. An attractive way of designing these sets is to consider a ball in the space of probability distributions centered at a reference or most-likely distribution constructed from the available data. In the space of distributions, the popular distance metric is the Prokhorov metric [7], $\phi$-divergence [8] and the Wasserstein distance [3]. In particular, the work [3] presents a tractable reformulation of DRO via Wasserstein balls, and is extended in [5] to a distributed setting. However, the available problem reformulation in [3] and the distributed algorithm in [5] do not consider the update of the decision over time and streaming data, which is the focus of this work. To design a tractable algorithm incorporating streaming data, our work connects to various convex optimization methods [9, 10] such as the Frank-Wolfe (FW) Algorithm (e.g., conditional gradient algorithm), the Subgradient Algorithm, and their variants, see e.g. [11, 12, 13] and references therein. Our emphasis on the convergence of the data-driven decision obtained through a sequence of optimization problems contrasts with typical algorithms developed for single (non-updated) problems.

Statement of Contributions. In this paper, we propose a new Online Data Assimilation Algorithm (OnDA Algorithm) to solve decision-making problems subject to uncertainty. The distribution of the uncertainty is unknown and the algorithm adjusts decisions based on realizations of the stochastic variable sequentially revealed over time. The new algorithm addresses four challenges: 1) the evaluation of the out-of-sample performance of every possible online decision; 2) the adaptation to online, increasingly-larger data sets to reach a decision with performance guarantees with increasingly higher probabilities; 3) the availability of an online decision vector with performance guarantees at any time; 4) the capability of handling sufficiently large streaming data sets.

To address 1), we start from a DRO problem setting. This leads to a worst-case optimization over an ambiguity set or neighborhood of the empirical distribution constructed from a data set. To solve this intractable problem, we reformulate it into an equivalent convex optimization over a simplex. This enables us to explore the simplex vertex set to find a certificate (a value bounding the cost) of a given decision with certain confidence. When the data is streaming, we consider a sequence of DRO problems and their equivalent convex reformulations employing increasingly larger data sets. Thus, as the data streams, the associated problems are defined over simplices of increasingly larger dimension. The similarities of the feasible sets allow us to assimilate the data via specialized Frank-Wolfe Algorithm variants, thus solving 2) via a Certificate Generation Algorithm (C-Gen Algorithm) described in Section V. Further, to seek for decisions that approach to the minimizers of the optimization problem, the OnDA Algorithm adapts its iterations online via a Subgradient Algorithm as described in Section VI. We show in Section VII that the resulting OnDA Algorithm is finitely convergent in the sense that the confidence of the out-of-sample performance guarantee for the generated data-driven decision converges to $1$ as the number of data samples increases to a sufficiently large but finite value. Under this scheme, a data-driven decision with certain performance guarantee is also available any time as soon as the algorithm finishes generating the first certificate for the initial decision, which resolves the challenge 3). To expedite the algorithm and deal with challenge 4), we develop in Section VIII an Incremental Covering Algorithm (I-Cover Algorithm) to obtain low-dimensional ambiguity sets. These new sets are based on a weighted version of the empirical distribution and thus close to the full empirical distribution of the data. We finally illustrate the performance of the proposed OnDA Algorithm in Section IX, with and without the I-Cover Algorithm. A preliminary study of this work has appeared in [14].

II Preliminaries

This section introduces basic notations and convexity definitions111Let m, $\mathbb{R}_{\geq 0}^{m}$ and m×d denote respectively the $m$-dimensional Euclidean space, the $m$-dimensional nonnegative orthant, and the space of $m\times d$ matrices. We let $\boldsymbol{x}\in^{m}$ denote a column vector of dimension $m$, while ${\boldsymbol{x}}^{\top}$ represents its transpose. We say $\boldsymbol{x}\geq 0$, if all its entries are nonnegative. We use the shorthand notations $\boldsymbol{0}_{m}$ for the column vector ${(0,\cdots,0)}^{\top}\in^{m}$, $\boldsymbol{1}_{m}$ for the column vector ${(1,\cdots,1)}^{\top}\in^{m}$, and $\boldsymbol{I}_{m}\in^{m\times m}$ for the identity matrix. We use either subscripts or parentheses superscripts to index vectors, i.e., $\boldsymbol{x}_{k}\in^{m}$ or $\boldsymbol{x}^{(k)}\in^{m}$, for $k\in\{1,2,\ldots,n\}$. We use $(\boldsymbol{x},\boldsymbol{y})\in^{m+d}$ to indicate the concatenated column vector from $\boldsymbol{x}\in^{m}$ and $\boldsymbol{y}\in^{d}$. The $1$-norm of the vector $\boldsymbol{x}\in^{m}$ is denoted by $\|\boldsymbol{x}\|$. We define the $m$-dimensional Euclidean ball centered at $\boldsymbol{x}\in^{m}$ with radius $\omega$ as the set $B_{\omega}(\boldsymbol{x}):=\{\boldsymbol{y}\in^{m}\;|\;\|\boldsymbol{y}-\boldsymbol{x}\|\leq\omega\}$. Given a set of points $I$ in m, we let $\operatorname{conv}(I)$ indicate its convex hull. The gradient of a real-valued function $f:^{m}\rightarrow$ is written as $\nabla_{\boldsymbol{x}}f(\boldsymbol{x})$. The $i^{\textup{th}}$ component of the gradient vector is denoted by $\nabla_{i}f(\boldsymbol{x})$. We use $\operatorname{dom}f$ to denote the domain of the function $f$, i.e., $\operatorname{dom}f:=\{\boldsymbol{x}\in^{m}\;|\;-\infty<f(\boldsymbol{x})<+\infty\}$. We call the function $f$ proper if $\operatorname{dom}f\neq\O$. We say a function $F:\mathcal{X}\times\mathcal{Y}\rightarrow$ is convex-concave on $\mathcal{X}\times\mathcal{Y}$ if, for any point $(\tilde{\boldsymbol{x}},\tilde{\boldsymbol{y}})\in\mathcal{X}\times\mathcal{Y}$, $\boldsymbol{x}\mapsto F(\boldsymbol{x},\tilde{\boldsymbol{y}})$ is convex and $\boldsymbol{y}\mapsto F(\tilde{\boldsymbol{x}},\boldsymbol{y})$ is concave. We use the notation $\operatorname{sgn}:\rightarrow$, $x\mapsto\{-1,0,1\}$ denote the sign function. Finally, the projection operator $\operatorname{proj}_{\mathcal{Y}}(\mathcal{X}):\mathcal{X}\rightarrow\mathcal{Y}$ projects the set $\mathcal{X}$ onto $\mathcal{Y}$ under the Euclidean norm., including some from Probability Theory to describe the distributionally robust optimization framework following [3].

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space, with $\Omega$ the sample space, $\mathcal{F}$ a $\sigma$-algebra on $\Omega$, and $\mathbb{P}$ the associated probability distribution. Let $\xi:\Omega\rightarrow^{m}$ be an induced multivariate random variable. We denote by $\mathcal{Z}\subseteq^{m}$ the support of the random variable $\xi$ and denote by $\mathcal{M}(\mathcal{Z})$ the space of all probability distributions supported on $\mathcal{Z}$ with finite first moment. In particular, $\mathbb{P}\in\mathcal{M}(\mathcal{Z})$. To measure the distance between distributions in $\mathcal{M}(\mathcal{Z})$, in this paper we use the dual characterization of the Wasserstein metric [15] $d_{W}:\mathcal{M}(\mathcal{Z})\times\mathcal{M}(\mathcal{Z})\rightarrow\mathbb{R}_{\geq 0}$, defined by

[TABLE]

where $\mathcal{L}$ is the space of all Lipschitz functions defined on $\mathcal{Z}$ with Lipschitz constant 1. A closed Wasserstein ball of radius $\omega$ centered at a distribution $\mathbb{P}\in\mathcal{M}(\mathcal{Z})$ is denoted by $\mathbb{B}_{\omega}(\mathbb{P}):=\{\mathbb{Q}\in\mathcal{M}(\mathcal{Z})\;|\;d_{W}(\mathbb{P},\mathbb{Q})\leq\omega\}$. We define the Dirac measure at $x_{0}\in\Omega$ as $\delta_{\{x_{0}\}}:\Omega\rightarrow\{0,1\}$. For any set $A\in\mathcal{F}$, we let $\delta_{\{x_{0}\}}(A)=1$, if $x_{0}\in A$, otherwise [math].

III Problem Description

Consider a decision-making problem under of the form

[TABLE]

where $\boldsymbol{x}\in^{d}$ is the decision variable, the uncertainty random variable $\xi:\Omega\rightarrow^{m}$ is induced by the probability space $(\Omega,\mathcal{F},\mathbb{P})$, and the expectation of $f$ is taken w.r.t. the unknown distribution $\mathbb{P}\in\mathcal{M}(\mathcal{Z})$. We aim to develop an Online Data Assimilation Algorithm (OnDA Algorithm) that efficiently adapts iterations on decisions $\boldsymbol{x}$ of (P) with streaming data. The streaming data are sequentially available iid realizations of the random variable $\xi$ under $\mathbb{P}$, denoted by ${\xi}_{n}$, $n=1,2,\ldots$. This defines a sequence of streaming data sets, ${\Xi}_{n}\subseteq{\Xi}_{n+1}$, for each $n$. W.l.o.g. assume that each ${\Xi}_{n+1}$ consists of just one more new data point, i.e., ${\Xi}_{n+1}={\Xi}_{n}\cup\{{\xi}_{n+1}\}$ and ${\Xi}_{1}=\{{\xi}_{1}\}$. In the following, we refer to the time slot between the updates ${\Xi}_{n}$ and ${\Xi}_{n+1}$ as the $n^{\textup{th}}$-time period and to its rate of change as the data-streaming rate.

In practice, we cannot evaluate the objective function of (P) because $\mathbb{P}$ is unknown. We call a decision $\boldsymbol{x}\in^{d}$ a proper data-driven decision of (P), if its out-of-sample performance, defined by ${\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x},\xi)]}$, satisfies the performance guarantee

[TABLE]

where the expected cost upper bound or certificate ${J}_{n}(\boldsymbol{x})$ is a function that indicates the goodness of $\boldsymbol{x}$ under the data set ${\Xi}_{n}$. If $\boldsymbol{x}$ is adopted during the $n^{\textup{th}}$ time period, then ${\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x},\xi)]}\leq{J}_{n}(\boldsymbol{x})$ is an event that depends on the $n$ samples in $\Xi_{n}$, and $\mathbf{P}^{n}$ denotes the probability with respect to these. The confidence $1-\beta_{n}\in(0,1)\subset$ governs the choice of $\boldsymbol{x}$ and the resulting certificate ${J}_{n}(\boldsymbol{x})$. In words, the inequality (1) establishes that, given finite data $n$, the performance of the decision under the unknown distribution will not surpass the upper-bound certificate $J_{n}(\boldsymbol{x})$ with high probability. In the following section, we will determine the values ${J}_{n}$ via the solution of a parameterized maximization problem over $\boldsymbol{x}$. Therefore, finding an approximate certificate will be much easier than finding the exact one. Based on this, we call $\boldsymbol{x}$ $\epsilon_{1}$-proper, if it satisfies (1) with ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x})$ such that ${J}_{n}(\boldsymbol{x})\leq{J}_{n}^{\epsilon_{1}}(\boldsymbol{x})+\epsilon_{1}$. Thus, the approximates ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x})$ also provide upper bounds to the optimal value of (P) with high confidence $1-\beta_{n}$.

To sum up, for any $n$, given a confidence level $1-\beta_{n}$, our goal is to approach to an $\epsilon_{1}$-proper data-driven decision with a low certificate. Later in Section VI, we will show, under assumptions on $f$, that both Problem (P) and these certificates ${J}_{n}$ are convex. To find a decision with a low certificate, we will call any proper data-driven decision $\epsilon_{2}$-optimal, labeled as $\boldsymbol{x}_{n}^{\epsilon_{2}}$, if ${J}_{n}(\boldsymbol{x}_{n}^{\epsilon_{2}})-{J}_{n}(\boldsymbol{x})\leq\epsilon_{2}$ for all $\boldsymbol{x}\in^{d}$. Then, for any $\epsilon_{2}$-optimal and $\epsilon_{1}$-proper data-driven decision $\boldsymbol{x}_{n}^{\epsilon_{2}}$ with certificate ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})$ and $\epsilon_{1}\ll\epsilon_{2}$, we have the guarantee

[TABLE]

Then, any decision $\boldsymbol{x}_{n}^{\epsilon_{2}}$ ensures a high-confidence, potentially-low objective value of (P), upper bounded by ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})+\epsilon_{1}$.

III-1 Motivating Example in Portfolio Optimization

Consider an agent who does short-term trading, i.e., she is to select a minute-based portfolio weight $\boldsymbol{x}:=(x,1-x)$, $0\leq x\leq 1$, for two risky assets that give random returns, with return rates $\xi:=(\xi_{1},\xi_{2})$ following some unknown distribution $\mathbb{P}$. The agent aims to select $\boldsymbol{x}$ such that the expected profit is maximized, or equivalently, she seeks to solve

[TABLE]

Assume that $\mathbb{P}$ is unknown and independent from the selection of the portfolio and that, at every minute, the agent has access to return rates, which are iid samples of $\mathbb{P}$. Due to the independence of $\mathbb{P}$ and $\boldsymbol{x}$, Problem (P0) is equivalent to $\min\limits_{0\leq x\leq 1}\mathbb{E}_{\mathbb{P}}[{-\xi_{1}x-\xi_{2}(1-x)-{\xi}^{\top}\xi}].$ To adapt this problem to our unconstrained setting, we consider an approximation

[TABLE]

where $\rho>0$ is some penalty for the constraint terms. This problem is in form (P), fitting in our problem setting with

[TABLE]

Our proposed algorithm allows the agent to make online decisions $\boldsymbol{x}_{n}^{\epsilon_{2}}$ that minimize the objective with high confidence.

III-A Main Algorithm Goal and its High-level Procedure

We describe now the goal of the OnDA Algorithm  that handles a streaming sequence of data sets with $n\in\{1,\dots,N\}$. Let tolerance parameters $\epsilon_{1}$ and $\epsilon_{2}$ be given and let us choose strictly increasing confidence levels $\{1-\beta_{n}\}_{n=1}^{N}$ such that $\sum_{n=1}^{\infty}\beta_{n}<\infty$ whenever $N\rightarrow\infty$. The algorithm aims to find a sequence of $\epsilon_{2}$-optimal and $\epsilon_{1}$-proper decisions $\{\boldsymbol{x}_{n}^{\epsilon_{2}}\}_{n=1}^{N}$ associated with the sequence of the certificates $\{{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})\}_{n=1}^{N}$ so that (2) holds for all $n\in\{1,\dots,N\}$. Additionally, as the data set streams to infinite cardinality, i.e., $N\rightarrow\infty$, there exists a large enough but finite $n_{0}$ such that the algorithm returns a final $\boldsymbol{x}_{n_{0}}^{\epsilon_{2}}$ after processing the data set ${\Xi}_{n_{0}}$. The final decision $\boldsymbol{x}_{n_{0}}^{\epsilon_{2}}$ guarantees performance almost surely, that is, $\mathbf{P}^{n_{0}}({\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x}_{n_{0}}^{\epsilon_{2}},\xi)]}\leq{J}_{n_{0}}^{\epsilon_{1}}(\boldsymbol{x}_{n_{0}}^{\epsilon_{2}})+\epsilon_{1})=1$, with a certificate ${J}_{n_{0}}^{\epsilon_{1}}(\boldsymbol{x}_{n_{0}}^{\epsilon_{2}})$ close to the optimal objective value of Problem (P).

To achieve this goal, the algorithm will output a sequence of decisions, $\{\boldsymbol{x}^{(r)}\}_{r=1}^{\infty}$, on a time scale that is faster than the data-streaming rate. We refer to this as the decision-update rate and we use a parenthetical superscript $(r)$ to denote its iterations. New data arrival will reset the OnDA Algorithm’s subroutines to update the sub-sequences of decisions in each $n^{\textup{th}}$ time period efficiently. We denote these sub-sequences as $\{\boldsymbol{x}^{(r)}\}_{r=r_{n}}^{r_{n+1}}$, where $\boldsymbol{x}^{(r_{n})}$ is the initial decision adapted from the the time period $n-1$. These updates will require the computation of certificates (cost upper bounds) and the progressive reduction of these bounds.

The computation of certificates is carried out by the Certificate Generation Algorithm (or C-Gen Algorithm, for short). Given a current decision value, $\boldsymbol{x}^{(r)}$, and data set $\Xi_{n}$, the C-Gen Algorithm finds an $\epsilon_{1}$ certificate $J_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})$ and a worst-case distribution associated with the data set. It operates on a faster time scale than the decision-update rate, the so-called certificate-generation rate. Upon the receipt of new data, this algorithm will reset as described in Section V.

The second process of the OnDA Algorithm relies on iterating decisions to reduce the values of the functions $J_{n}^{\epsilon_{1}}(\boldsymbol{x})$. This employs the Subgradient Algorithm and is described in Section VI. A more thorough description of how new data triggers a reset in the algorithm is described in the following sections. A summary of the OnDA Algorithm can be found in Section VII as well as a descriptive table.

IV Certificate Design

In this section, we present a tractable formulation of certificates ${J}_{n}(\boldsymbol{x})$ and its approximation ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x})$ for a fixed $\boldsymbol{x}=\boldsymbol{x}^{(r)}$, as described in (1) and (2), respectively. To achieve this, we first follow [5, 6, 3] on DRO to find certificates ${J}_{n}$. This defines a parameterized maximization problem for ${J}_{n}$, called Problem (P1n). Then we reformulate (P1n) as Problem (P2n), a convex optimization problem over a simplex, for efficient solutions of approximated certificates ${J}_{n}^{\epsilon_{1}}$ in the next section.

To design ${J}_{n}$, a reasonable attempt is to use the data ${\Xi}_{n}$ to estimate an empirical distribution, $\hat{\mathbb{P}}^{n}$, and let ${\mathbb{E}_{\mathbb{\hat{P}^{\it n}}}[f(\boldsymbol{x},\xi)]}$ be the candidate certificate for the performance guarantee (1). More precisely, assume that the data set ${\Xi}_{n}$ are uniformly sampled from $\mathbb{P}$. The discrete empirical probability measure associated with ${\Xi}_{n}$ is the following $\hat{\mathbb{P}}^{n}:=\frac{1}{n}\sum_{k=1}^{n}\delta_{\{{\xi}_{k}\}},$ where $\delta_{\{{\xi}_{k}\}}$ is the Dirac measure at ${\xi}_{k}$. The candidate certificate is

[TABLE]

The above approximation $\hat{\mathbb{P}}^{n}$ of $\mathbb{P}$, also known as the sample-average estimate, makes $J^{\textup{sae}}_{n}$ easy to compute. However, such value only results in an approximation of the unknown out-of-sample performance $\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x},\xi)]$. Following [5, 3], we are to determine an ambiguity set $\mathcal{P}_{n}$ containing all the possible probability distributions supported on $\mathcal{Z}\subseteq^{m}$ that can generate ${\Xi}_{n}$ with high confidence. Then, one can consider the worst-case expectation of $f(\boldsymbol{x},\xi)$ with respect to all distributions contained in $\mathcal{P}_{n}$. The solution to such problem offers an upper bound for the out-of-sample performance with high probability in the form of (1), and we refer to this upper bound as the certificate of the decision $\boldsymbol{x}$.

In order to quantify the ambiguity set and certificate for an $\epsilon_{1}$-proper data-driven decision, we denote by $\mathcal{M}_{\textup{lt}}(\mathcal{Z})\subset\mathcal{M}(\mathcal{Z})$ the set of light-tailed probability measures in $\mathcal{M}(\mathcal{Z})$, and introduce the following assumption for $\mathbb{P}$ {assumption}[Light tailed unknown distributions] It holds that $\mathbb{P}\in\mathcal{M}_{\textup{lt}}(\mathcal{Z})$, i.e., there exists an exponent $a>1$ such that $b:=\mathbb{E}_{\mathbb{P}}[\exp(\|\xi\|^{a})]<\infty$.

Remark IV.1** (Class of distributions satisfying Assumption IV).**

Any distribution with an exponentially decaying tail satisfies this assumption, such as Gaussian, subGaussian, exponential, and geometric distributions. Any distribution with a compact support $\mathcal{Z}$ will trivially satisfy the assumption. In engineering problems, the values of random variables are usually truncated to a compact set and hence the Assumption IV is automatically satisfied.* *

Assumption IV validates the following modern measure concentration result, which provides an intuition for considering the Wasserstein ball $\mathbb{B}_{\epsilon}(\hat{\mathbb{P}}^{n})$ of center $\hat{\mathbb{P}}^{n}$ and radius $\epsilon$ as the ambiguity set $\mathcal{P}_{n}$.

Theorem IV.1** (Measure concentration [16, Theorem 2]).**

If $\mathbb{P}\in\mathcal{M}_{\textup{lt}}(\mathcal{Z})$, then

[TABLE]

*for all $n\geq 1$, $m\neq 2$, and $\epsilon>0$, where $c_{1}$, $c_{2}$ are positive constants that only depend on $m$, $a$ and $b$. $\square$

Equipped with this result, we are able to provide the certificate that ensures the performance guarantee in (1), for any decision $\boldsymbol{x}\in^{d}$.

Lemma IV.1** (Certificate in Performance

Guarantee (1)).**

Given ${\Xi}_{n}:=\{{\xi}_{k}\}_{k=1}^{n}$, $\beta_{n}\in(0,1)$ and $\boldsymbol{x}\in^{d}$, let

[TABLE]

and $\mathcal{P}_{n}:=\mathbb{B}_{\epsilon(\beta_{n})}(\hat{\mathbb{P}}^{n})$. Then the following certificate satisfies the performance guarantee in (1) for all $\boldsymbol{x}\in^{d}$

[TABLE]

To get ${J}_{n}$ in (5), one needs to solve an infinite-dimensional optimization problem. Luckily, Problem (5) can be reformulated into a finite-dimensional convex problem as follows.

Theorem IV.2** (Convex reduction of (5) [3, Application of Theorem

4.4]).**

Under Assumption IV, on $\mathbb{P}$ being light-tailed, for all $\beta_{n}\in(0,1)$ the value of the certificate in (5) for a given decision $\boldsymbol{x}$ under a data set ${\Xi}_{n}$ is equal to the optimal value of the following optimization problem

[TABLE]

where each component of the concatenated variable $\boldsymbol{y}:=(\boldsymbol{y}_{1},\ldots,\boldsymbol{y}_{n})$ is in m, and the parameter $\epsilon(\beta_{n})$ is the radius of $\mathbb{B}_{\epsilon(\beta_{n})}$ calculated from (4). Moreover, given any feasible point $\boldsymbol{y}^{(l)}:=(\boldsymbol{y}^{(l)}_{1},\ldots,\boldsymbol{y}^{(l)}_{n})$ of (P1n), indexed by $l$, define a finite atomic probability measure at $\boldsymbol{x}$ in the Wasserstein ball $\mathbb{B}_{\epsilon(\beta_{n})}$ of the form

[TABLE]

*Now, denote by $\mathbb{Q}_{n}^{\star}(\boldsymbol{x})$ the distribution in (6) constructed by an optimizer $\boldsymbol{y}^{\star}:=(\boldsymbol{y}^{\star}_{1},\ldots,\boldsymbol{y}^{\star}_{n})$ of (P1n) and evaluated over $\boldsymbol{x}$. Then, $\mathbb{Q}_{n}^{\star}$ is a worst-case distribution that can generate the data set ${\Xi}_{n}$ with (high) probability no less than $1-\beta_{n}$. $\square$

Remark IV.2**.**

Theorem IV.2 provides a way of computing certificates of (1) as the solution to a parameterized optimization problem for a decision $\boldsymbol{x}$. In addition, it constructs a worst-case distribution that achieves the worst-case bound.* *

To enable efficient online solutions of approximated certificates ${J}_{n}^{\epsilon_{1}}$, let us define parameterized functions $h_{k}:^{m}\rightarrow$

$h_{k}(\boldsymbol{y}):=f(\boldsymbol{x},{\xi}_{k}-\boldsymbol{y}),\quad{{k\in\{1,\dots,n\}}}$,

and consider the following convex optimization problem over a simplex

[TABLE]

where the concatenated variable $(\boldsymbol{u},\boldsymbol{v})$ is composed of $\boldsymbol{u}:=(\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{n})$ and $\boldsymbol{v}:=(\boldsymbol{v}_{1},\ldots,\boldsymbol{v}_{n})$ with $\boldsymbol{u}_{k}$, $\boldsymbol{v}_{k}\in^{m}$ for all $k\in\{1,\dots,n\}$; and the scalar $n\epsilon(\beta_{n})$ regulates the size of the feasible set via scaling of the unit simplex $\Delta_{2mn}:=\{(\boldsymbol{u},\boldsymbol{v})\in^{2mn}\;|\;{\boldsymbol{1}_{2mn}}^{\top}(\boldsymbol{u},\boldsymbol{v})=1,\;\boldsymbol{u}\geq 0,\;\boldsymbol{v}\geq 0\}$. We denote by ${\Lambda}_{2mn}$ the set of all the extreme points for the simplex $n\epsilon(\beta_{n})\Delta_{2mn}$.

The following lemma shows that Problem (P1n) and Problem (P2n) are equivalent. Thus, we can approximately solve (P2n) to find ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x})$ and $\mathbb{Q}_{n}^{\epsilon_{1}}(\boldsymbol{x})$.

Lemma IV.2** (Equivalence of the problem formulation).**

Solving (P1n) is equivalent to solving (P2n) in the sense that

• 1

For any feasible solution $(\tilde{\boldsymbol{u}},\tilde{\boldsymbol{v}})$ of (**P2n***), let $\tilde{\boldsymbol{y}}:=\tilde{\boldsymbol{u}}-\tilde{\boldsymbol{v}}$. Then $\tilde{\boldsymbol{y}}$ is feasible for (P1n).*

• 2

For any feasible solution $\tilde{\boldsymbol{y}}$ of (**P1n***), there exists a feasible point $(\tilde{\boldsymbol{u}},\tilde{\boldsymbol{v}})$ of (P2n).*

• 3

Assume that the point $(\tilde{\boldsymbol{u}}^{\star},\tilde{\boldsymbol{v}}^{\star})$ is an optimizer of (**P2n***). Then by letting $\tilde{\boldsymbol{y}}^{\star}:=\tilde{\boldsymbol{u}}^{\star}-\tilde{\boldsymbol{v}}^{\star}$, the point $\tilde{\boldsymbol{y}}^{\star}$ is also an optimizer of (P1n), with the same optimal value.*

For this section, we provide proofs of Lemma 5 and Lemma IV.2 in the Appendix, and from now on, all the proofs of the lemmas and theorems can be found in the Appendix.

V Certificate Generation Algorithm

Given a tolerance $\epsilon_{1}$, sequentially available data sets $\{{\Xi}_{n}\}_{n=1}^{N}$ and decisions $\{\boldsymbol{x}^{(r)}\}_{r=1}^{\infty}$, we present in this section the Certificate Generation Algorithm  (C-Gen Algorithm) to obtain approximated certificates $\{{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})\}_{n,r}$ and associated $\epsilon_{1}$-worst-case distributions $\{\mathbb{Q}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})\}_{n,r}$. To achieve this, we first design for each fixed $\boldsymbol{x}=\boldsymbol{x}^{(r)}$ the C-Gen Algorithm to solve (P2n) to ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x})$ efficiently. This is developed via Frank-Wolfe Algorithm variants, e.g., the Simplicial Algorithm [12] and the AFWA as described in the Appendix. Then we analyze the convergence of the C-Gen Algorithm under $\{{\Xi}_{n}\}_{n=1}^{N}$.

V-A The C-Gen Algorithm

For each fixed $\boldsymbol{x}=\boldsymbol{x}^{(r)}\in^{d}$ the algorithm is run at a fast time scale (the certificate-update rate), over iterations $l=0,1,2,\dots$. The algorithm is then employed inside the OnDA Algorithm, so its execution rate is the fastest within this algorithm. At each iteration $l$, the C-Gen Algorithm generates $(\boldsymbol{u}^{(l)},\boldsymbol{v}^{(l)})$, the candidate optimizer of (P2n). Let the objective value of (P2n) at $(\boldsymbol{u}^{(l)},\boldsymbol{v}^{(l)})$ be ${J}_{n}^{(l)}(\boldsymbol{x})$, and, equivalently, write the candidate optimizer in form of $\boldsymbol{y}^{(l)}:=\boldsymbol{u}^{(l)}-\boldsymbol{v}^{(l)}$ (exploiting the equivalence in Lemma IV.2). Each candidate $\boldsymbol{y}^{(l)}$ is associated with a set of search points denoted by $I_{n}^{(l)}:=\{\tilde{\boldsymbol{y}}^{[i]}:=\tilde{\boldsymbol{u}}^{[i]}-\tilde{\boldsymbol{v}}^{[i]},\;i\in\{1,\dots,T\}\}$, where we use bracket superscript $[i]$ to index its elements. As we will see later, the set $I_{n}^{(l)}$ plays a key role in generating the certificate when assimilating data, and is called the candidate vertex set.

Given a data set $\Xi_{n}$, and until new data arrives, the C-Gen Algorithm solves the following problems alternatively

[TABLE]

[TABLE]

(Note how the solution to one problem parameterizes the other.) In this way, $\boldsymbol{y}^{(l-1)}$ (the solution to the $\text{CP}_{n}^{(l-1)}$) parameterizes the linear problem (LP${}_{n}^{(l)}$). The solution to (LP${}_{n}^{(l)}$) is then used to refine the set point $I_{n}^{(l)}=\{\tilde{\boldsymbol{y}}^{[i]}\}_{i}$, which spans the constraint set $\Delta_{T}\equiv\operatorname{conv}(I_{n}^{(l)})$ in problem (CP${}_{n}^{(l)}$). A solution to (CP${}_{n}^{(l)}$) then determines the new $\boldsymbol{y}^{(l)}$ of the next LP problem. This process corresponds to lines 3: to 9: in the following C-Gen Algorithm table.

More precisely, at each iteration $l=1,2,\ldots$, the C-Gen Algorithm first solves (LP${}_{n}^{(l)}$) using the Point Search Algorithm, which returns the optimal objective value $\eta^{(l)}$ and the set of maximizers $\Omega^{(l)}$ such that $\eta^{(l)}\geq{J}_{n}(\boldsymbol{x})-{J}_{n}^{(l)}(\boldsymbol{x})$ and $\Omega^{(l)}\subset\Lambda_{2mn}$. The value $\eta^{(l)}$ is then used to determine the $\epsilon_{1}$-suboptimality condition to the optimal objective of Problem (P2n) (see below). Meanwhile, the set $\Omega^{(l)}$ is used to update candidate vertex set to $I_{n}^{(l)}:=I_{n}^{(l-1)}\cup\Omega^{(l)}$, which is used in problem (CP${}_{n}^{(l)}$). In particular, the Point Search Algorithm computes all optimizers by iteratively choosing a sparse vector with only a positive entry. That is, an extreme point of the feasible set of (LP${}_{n}^{(l)}$), such that the nonzero component of $(\tilde{\boldsymbol{u}}^{(l)},\tilde{\boldsymbol{v}}^{(l)})$ has the largest absolute gradient component in the linear cost function of (LP${}_{n}^{(l)}$). Using the obtained $I_{n}^{(l)}$, the algorithm solves the Problem (CP${}_{n}^{(l)}$) over the simplex $\Delta_{T}:=\{\gamma\in^{T}\;|\;{\boldsymbol{1}_{T}}^{\top}\gamma=1,\;\gamma\geq 0\}$, where $T$ is the cardinality of $I_{n}^{(l)}$ and each component $\gamma_{i}$ of $\gamma\in\Delta_{T}$ represents the convex combination coefficient of a candidate vertex \tilde{\boldsymbol{y}}^{[i]}$$\in I_{n}^{(l)}. After solving (CP${}_{n}^{(l)}$) to $\epsilon_{1}$-optimality via the AFWA (see Appendix), an $\epsilon_{1}$-optimal weighting $\gamma^{\epsilon_{1}}\in\Delta_{T}$ with the objective value ${J}_{n}^{(l)}(\boldsymbol{x})$ is obtained. A new candidate optimizer $\boldsymbol{y}^{(l)}$ is then calculated by $\boldsymbol{y}^{(l)}=\sum_{i=0}^{T}\gamma_{i}^{\epsilon_{1}}\tilde{\boldsymbol{y}}^{[i]}$. The algorithm repeats the process and increments $l$ if the optimality gap $\eta^{(l)}$ is greater than $\epsilon_{1}$, otherwise it returns the certificate ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}):={J}_{n}^{(l)}(\boldsymbol{x})$, an $\epsilon_{1}$-optimal solution ${\boldsymbol{y}}^{\epsilon_{1}}:=\boldsymbol{y}^{(l)}$ and an $\epsilon_{1}$-worst-case distribution $\mathbb{Q}_{n}^{\epsilon_{1}}(\boldsymbol{x}):=\frac{1}{n}\sum_{k=1}^{n}\delta_{\{{\xi}_{k}-{\boldsymbol{y}}^{\epsilon_{1}}_{k}\}}$.

When new data arrives, the algorithm will reset by adapting the Problem (P2n+1) from Problem (P2n) (line 2: in the table of the C-Gen Algorithm (update from $\Xi_{n}$ to $\Xi_{n+1}$). Note that adapting the C-Gen Algorithm to online data sets $\{{\Xi}_{n}\}_{n=1}^{N}$ is inherently difficult due to the changes in the Problems (P2n). As the size of ${\Xi}_{n}$ grows by $1$, the dimension of the Problem (P2n) increases by $2m$. To obtain ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x})$ and $\mathbb{Q}_{n}^{\epsilon_{1}}(\boldsymbol{x})$ sufficiently fast, we exploit the relationship among Problems (P2n), for different $n$, by adapting the candidate vertex sets $I_{n}^{(l)}$. Specifically, we initialize the set $I_{n+1}^{(0)}$ for the new Problem (P2n+1) by $I_{n}^{(l)}$, constructed from the previous (P2n). Suppose that the C-Gen Algorithm receives a new data set ${\Xi}_{n+1}\supset{\Xi}_{n}$ at some intermediate iteration $l$ with candidate vertex set $I_{n}^{(l)}$. At this stage, the subset $\operatorname{conv}(I_{n}^{(l)})$ has been explored by the previous optimization problem, and the gradient information of the objective function based on the data set ${\Xi}_{n}$ has been partially integrated. Then, by projecting the set $I_{n}^{(l)}$ onto the set of extreme points of the new Problem (P2n), i.e., $I_{n+1}^{(0)}:=\operatorname{proj}_{\Lambda_{2m(n+1)}}(\{(\tilde{\boldsymbol{y}}^{[i]},\boldsymbol{0}_{m})\;|\;\tilde{\boldsymbol{y}}^{[i]}\in I_{n}^{(l)}\})$, the subset $\operatorname{conv}(I_{n+1}^{(0)})$ of the feasible set of (P2n) is already explored. Such integration contributes to the reduction of the number of iterations in the C-Gen Algorithm for Problems (P2n). This insight gives us a sense of the worst-case efficiency to update a certificate under the streaming data.

V-B Convergence Analysis of the C-Gen Algorithm

We make the following assumptions on the local strong concavity of the function $f$ and the computation of its gradient {assumption}[Local strong concavity] For any ${\boldsymbol{x}}\in^{d}$ and ${\xi}\in^{m}$, the function $h:^{m}\rightarrow$, $\boldsymbol{y}\mapsto f({\boldsymbol{x}},{\xi}-\boldsymbol{y})$ is differentiable, concave with a curvature constant $C_{h}$, and with a positive geometric strong concavity constant $\mu_{h}$ on $\Delta_{2mN}$222 For a concave function $h:^{m}\rightarrow$ on $\Delta$, we define $\displaystyle{C_{h}:=\sup-\frac{2}{\gamma^{2}}\left(h(\boldsymbol{y}^{\star})-h(\boldsymbol{y})-\left\langle\nabla h(\boldsymbol{y}),\boldsymbol{y}^{\star}-\boldsymbol{y}\right\rangle\right),}$ s.t. $\boldsymbol{y}^{\star}=\boldsymbol{y}+\gamma(\boldsymbol{s}-\boldsymbol{r}),$ $\gamma\in[0,1]$, $\boldsymbol{y},\boldsymbol{s},\boldsymbol{r}\in{\Delta}.$ and $\displaystyle{\mu_{h}:=\inf\limits_{\boldsymbol{y}\in\Delta}\inf\limits_{\boldsymbol{y}^{\star}\in\Delta}-\frac{2}{\Gamma(\boldsymbol{y},\boldsymbol{y}^{\star})^{2}}\times}$ $\left(h(\boldsymbol{y}^{\star})-h(\boldsymbol{y})-\left\langle\nabla h(\boldsymbol{y}),\boldsymbol{y}^{\star}-\boldsymbol{y}\right\rangle\right),$ s.t. $\left\langle\nabla h(\boldsymbol{y}),\boldsymbol{y}^{\star}-\boldsymbol{y}\right\rangle>0,$ where $\Gamma(\boldsymbol{y},\boldsymbol{y}^{\star})$ is a step-size measure in AFWA. See, e.g., [13] for details. We say $h$ is locally strongly concave, if $\mu_{h}>0$..

{assumption}

[ Accessible gradients] For any decision ${\boldsymbol{x}}\in^{d}$, we denote by $\nabla h(\boldsymbol{y})$ the gradient of the function $h:^{m}\rightarrow$, $\boldsymbol{y}\mapsto f({\boldsymbol{x}},\boldsymbol{y})$ and assume it is accessible.

Under Assumptions 2 and V-B, we show the convergence properties of the C-Gen Algorithm.

Theorem V.1** (Convergence of the C-Gen Algorithm).**

Let a tolerance $\epsilon_{1}$ and a decision $\boldsymbol{x}$ be given. Let us choose $\boldsymbol{y}^{(0)}=\boldsymbol{0}_{m}$ and $I_{1}^{(0)}=\O$ as the initial candidate optimizer and candidate vertex set for the C-Gen Algorithm, respectively. Consider the online data sets $\{{\Xi}_{n}\}_{n=1}^{N}$ and the set of parameterized functions $\{h_{n}\}_{n=1}^{N}$. Under Assumption 2 and Assumption V-B, we have that for all data set ${\Xi}_{n}$, there exists a parameter $\kappa\in(0,1)\subset$ such that the worst-case computational bound $\phi(n)$ of the C-Gen Algorithm, depending on $n$, is

[TABLE]

Moreover, consider that data sets $\{{\Xi}_{n}\}_{n=1}^{N}$ are streaming and consider function $J^{\textup{sae}}_{N}(\boldsymbol{x})$ defined as in Section IV. Then there exists a parameter $\bar{\kappa}\in(0,1)\subset$ and a computational bound

[TABLE]

*such that, if the data-streaming rate is slower or equal than $(\bar{\phi}(1))^{-1}$, then the C-Gen Algorithm is guaranteed to obtain the certificates $\{{J}_{n}^{\epsilon_{1}}(\boldsymbol{x})\}_{n=1}^{N}$ and $\{\mathbb{Q}_{n}^{\epsilon_{1}}(\boldsymbol{x})\}_{n=1}^{N}$. *

Theorem V.1 relates the worst-case computational bound of the C-Gen Algorithm, executed on the certificate-update rate (the fastest of the time scales considered), to the data-streaming rate. Note that, as $\epsilon_{1}$ decreases, the bound $\bar{\phi}(1)$ increases and therefore the smaller the data-streaming rate has to be so that the certificates can be generated by the algorithm. In practice, the C-Gen Algorithm tends to find the smallest implicit feasible set that contains an optimal solution of (P2n). This means that the computation of the C-Gen Algorithm generally performs better than its worst-case bound as in Theorem V.1 and so it can handle data-streaming rates faster than $(\bar{\phi}(1))^{-1}$. In the sequel, we assume that the C-Gen Algorithm converges with a rate that is faster than the worst-case bound in Theorem V.1.

Remark V.1** (Effects of Assumption 2 and V-B).**

The essential ingredients for convergence of the Certificate Generation Algorithm are 1) the concavity of $h$, which ensures that (P2n) is a convex problem, and 2) accessible gradients of $h$, which allows for computations to a solution of (P2n). To obtain a fast, linear convergence rate as in Theorem V.1, we assume that $h$ is strongly concave on the simplex $\Delta_{2mN}$ located at each data point $\xi\in{\Xi}_{n}$. Intuitively, as ${\Xi}_{n}$ comes from $\mathbb{P}$, Assumption 2 eventually requires $h$ to be strongly concave on a subset of the support $\mathcal{Z}$ of $\mathbb{P}$ where the high-probability outcomes are concentrated onto. Otherwise, if $h$ is concave but not locally strongly concave, or if the gradients of $h$ are inaccessible (e.g., when only non-biased gradient estimate of $h$ are available), the convergence of AFWA, as described in Theorem .1, reduces to a sublinear rate. This, in turn, reduces the computational bound $\bar{\phi}(n)$ in Theorem V.1 to a bound of order $\mathcal{O}(1/\epsilon_{1})$.* *

Remark V.2** (Example in Portfolio Optimization).**

The portfolio problem in Section III results in a strongly concave $h$ which implies the local strong concavity as required by Assumption 2. Let $\boldsymbol{y}:=(y_{1},y_{2})$ and, for any given data point $\xi_{k}:=(\xi_{k,1},\xi_{k,2})$, $\nabla h$ is accessible and computed by

[TABLE]

VI Sub-optimal decisions with guarantees

In this section, we aim to construct a sub-sequence of $\epsilon_{2}$-optimal data-driven decisions $\{\boldsymbol{x}_{n}^{\epsilon_{2}}\}_{n=1}^{N}$, associated with the $\epsilon_{2}$-lowest certificates $\{{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})\}_{n=1}^{N}$ over time. We achieve this by means of the Subgradient Algorithm to derive an $\epsilon_{1}$-proper decision sequence $\{\boldsymbol{x}^{(r)}\}_{r=r_{n}}^{r_{n+1}}$; and the concatenation of $\{\boldsymbol{x}^{(r)}\}_{r=r_{n}}^{r_{n+1}}$ for different $n$ to obtain $\{\boldsymbol{x}_{n}^{\epsilon_{2}}\}_{n=1}^{N}$.

To construct an $\epsilon_{1}$-proper decision sub-sequence $\{\boldsymbol{x}^{(r)}\}_{r=r_{n}}^{r_{n+1}}$ let us consider the following problem

[TABLE]

where the function ${J}_{n}(\boldsymbol{x})$ is defined as in either (5) or (P1n), and we assume the approximation of ${J}_{n}(\boldsymbol{x})$, ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x})$, can be evaluated as in Section V.

To solve this Problem to ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})$, we have the following assumption on the convexity of $f$ {assumption}[Convexity in $\boldsymbol{x}$] The function $f_{\xi}:^{d}\rightarrow$ ${\boldsymbol{x}}\mapsto{f(\boldsymbol{x},\xi)}$ is convex for all $\xi\in^{m}$.

Assumption VI results in convexity of ${J}_{n}(\boldsymbol{x})$ as follows.

Lemma VI.1** (Convexity of ${J}_{n}(\boldsymbol{x})$).**

*If Assumption VI (convexity in $\boldsymbol{x}$) holds, then for each $n\in\{1,\dots,N\}$ the certificate ${J}_{n}(\boldsymbol{x})$ defined by (5) is convex in $\boldsymbol{x}$.

Lemma VI.1 allows us to apply the Subgradient Algorithm [17, 18, 19] to obtain $\boldsymbol{x}_{n}^{\epsilon_{2}}$ via $\{\boldsymbol{x}^{(r)}\}_{r=r_{n}}^{r_{n+1}}$ and the following lemma.

Lemma VI.2** (Easy estimate of the $\epsilon$-subgradients

of ${J}_{n}(\boldsymbol{x})$).**

Let the tolerance $\epsilon_{1}$ and time period $n$ be given. For any decision $\boldsymbol{x}^{(r)}$, we denote an $\epsilon_{1}$-optimal solution and $\epsilon_{1}$-worst-case distribution of (P1n) by $\boldsymbol{y}^{\epsilon_{1}}$ and $\mathbb{Q}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})$, respectively. Let us consider the function $g_{n}^{r}:^{d}\rightarrow^{d}$, defined as

[TABLE]

Denote an $\epsilon$-subdifferential of ${J}_{n}(\boldsymbol{x})$ at $\boldsymbol{x}$, by $\partial_{\epsilon}{J}_{n}(\boldsymbol{x})$. Then, for all $\epsilon\geq\epsilon_{1}$ we have the following

[TABLE]

or equivalently, for every $\boldsymbol{z}\in{\operatorname{dom}}\>{J}_{n}$ and $\epsilon\geq\epsilon_{1}$, we have

[TABLE]

Moreover, for any $\tilde{\boldsymbol{x}}\in^{d}$, there exists $\eta>0$ such that for all $\epsilon\geq\eta$ the following relation holds

[TABLE]

Note how Lemma VI.2 employs the discrete distribution $\mathbb{Q}_{n}^{\epsilon_{1}}$ generated from the C-Gen Algorithm in the computation of an $\epsilon$-subgradient function of ${J}_{n}$. Thus, the Subgradient Algorithm can be employed to reach an $\epsilon_{1}$-proper data-driven decision with a lower certificate.

To do this, we make use of the scaled $\epsilon$-subgradient direction for the update of decisions $\{\boldsymbol{x}^{(r)}\}_{r=r_{n}}^{r_{n+1}}$, as follows

[TABLE]

where the nonnegative step size rule $\{\alpha^{(r)}\}_{r}$ is determined in advance. Later in the next subsection we will see how the choice of a step size rule affects the convergence of the Subgradient Algorithm to an $\boldsymbol{x}_{n}^{\epsilon_{2}}$.

The Subgradient Algorithm requires access of $\{g_{n}^{r}\}_{r=r_{n}}^{r_{n+1}}$, which are obtained from C-Gen Algorithm. To reduce the number of computations, we estimate the candidate subgradient functions $\{g_{n}^{r}\}_{r=r_{n}}^{r_{n+1}}$ as follows. Let $\epsilon_{\textup{SA}}\geq\epsilon_{1}$ be a specified tolerance. At some iteration $r\geq r_{n}$, assume that an $\epsilon_{1}$-optimizer $\boldsymbol{y}^{\epsilon_{1}}$ and $\epsilon_{1}$-worst-case distribution $\mathbb{Q}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})$ are obtained from the C-Gen Algorithm. Using $\mathbb{Q}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})$, we calculate the function $g_{n}^{r}$ at $\boldsymbol{x}^{(r)}$ and perform the subgradient iteration (7). At iteration $r+1$ with $\boldsymbol{x}^{(r+1)}$, we firstly check for the suboptimality of Problem (P1${}_{n}^{(r+1)}$) using the initial candidate optimizer $\boldsymbol{y}^{(0)}:=\boldsymbol{y}^{\epsilon_{1}}$ in the Point Search Algorithm. If the optimality gap $\eta^{(1)}$ is less than $\epsilon_{\textup{SA}}$, we estimate the candidate subgradient function $g_{n}^{r+1}$ using $g_{n}^{r}$ and proceed with the subgradient iteration. Otherwise, we obtain $g_{n}^{r+1}$ from the C-Gen Algorithm, which is again an ${\epsilon}_{1}$-subgradient function at $\boldsymbol{x}^{(r+1)}$. Thus, we construct a sequence of $\epsilon_{\textup{SA}}$-subgradient functions $\{g_{n}^{r}\}_{r=r_{n}}^{r_{n+1}}$ that achieve an $\boldsymbol{x}_{n}^{\epsilon_{2}}$ efficiently.

Remark VI.1** (Effect of tolerance

$\epsilon_{\textup{SA}}$).**

The tolerance $\epsilon_{\textup{SA}}$ quantifies whether the function $g_{n}^{r}$ generated by the current worst-case distribution can also provide a good estimate of the ${\epsilon}$-subgradient at the next iteration point. If the function $g_{n}^{r}$ is an ${\epsilon}$-subgradient for ${\epsilon}$ small enough, there is no need of employing the C-Gen Algorithm to obtain a new subgradient function, which will be again an ${\epsilon_{1}}$-subgradient. In practice, we suggest to choose $\epsilon_{\textup{SA}}\gg\epsilon_{1}$ as it reduces the number of computations from the C-Gen Algorithm.* *

VI-A Convergence Analysis for the $\epsilon_{2}$-optimal decisions

The following lemma follows from the convergence of the Subgradient Algorithm applied to our problem scenario.

Lemma VI.3** (Convergence of $\epsilon_{\textup{SA}}$-Subgradient

Algorithm).**

For each time period $n$ with an initial data-driven decision $\boldsymbol{x}^{(r_{n})}$, assume that the subgradients defined in Lemma VI.2 are uniformly bounded, i.e., there exists a constant $L>0$ such that $\|g_{n}^{r}\|\leq L$ for all $r\geq r_{n}$.

Given a predefined $\epsilon_{2}>0$, let the certificate tolerance $\epsilon_{1}$ and the subgradient tolerance $\epsilon_{\textup{SA}}$ be such that $0<\epsilon_{1}\leq\epsilon_{\textup{SA}}<\epsilon_{2}/\mu$ with $\mu:=\max\{L,\,1\}$. Let $\boldsymbol{x}_{n}^{\star}\in\operatorname*{argmin}_{\boldsymbol{x}\in^{d}}{J}_{n}(\boldsymbol{x})$. Then, there exists a large enough number $\bar{r}$, depending on $\epsilon_{\textup{SA}}$ and the step size rule $\{\alpha^{(r)}\}_{r}$, such that the above designed Subgradient Algorithm in (7) has the following performance bounds

[TABLE]

and terminates at the iteration $r_{n+1}:=\bar{r}{{+r_{n}}}$ with an $\epsilon_{2}$-optimal decision by choosing $\boldsymbol{x}_{n}^{\epsilon_{2}}\in\operatorname*{argmin}\limits_{k\in\{{r_{n}},\dots,{\bar{r}{{+r_{n}}}}\}}\{{J}_{n}(\boldsymbol{x}^{(k)})\}$. In particular, there exists a large enough parameter $M$ such that we can select 1) a constant step-size rule given by

[TABLE]

where $\bar{r}:=M^{2}\left(\frac{\epsilon_{2}}{\mu}-\epsilon_{\textup{SA}}\right)^{-2}$; or 2) a divergent, but square-summable, step-size rule given by

[TABLE]

*where $\bar{r}=\min\{r\in\mathbb{N}\;|\;M(3-\frac{1}{r+1})\leq 2(\frac{\epsilon_{2}}{\mu}-\epsilon_{\textup{SA}})\ln(r+1)\}$.

In other words, Lemma VI.3 specifies that there is a finite, large enough iteration step at which the $\epsilon_{\textup{SA}}$-Subgradient Algorithm terminates using the estimated $\epsilon_{\textup{SA}}$-subgradient functions.

To quantify the effect of the subgradient estimation on the convergence rate under $\Xi_{n}$, we have the following theorem.

Theorem VI.1** (Worst-case computational bound for an

$\boldsymbol{x}_{n}^{\epsilon_{2}}$).**

For each time period $n$ with an initial $\boldsymbol{x}^{(r_{n})}$, let us consider the algorithm setting as in Lemma VI.3. Then, there exist parameters $\kappa\in(0,1)$, and $t>\epsilon_{1}$ such that the computational steps $\varphi(n,\bar{r})$ to reach $\boldsymbol{x}_{n}^{\epsilon_{2}}$ are bounded by

[TABLE]

*where $\bar{r}$ are the subgradient steps of Lemma VI.3. The value $\phi(n)$ is the worst-case computational bound as in Theorem V.1 and one should use $\bar{\phi}(1)$ in the bound in place of $\phi(n)$ if considering a data-streaming scenario. *

Theorem VI.1 integrates together the obtained bounds for the Subgradient Algorithm as well as the C-Gen Algorithm. As the result, a worst-case computational bound for the OnDA Algorithm, executed on the decision-update rate (the second time scale in $(r)$), is related to the data-streaming rate. Whenever the data-streaming rate is greater than the worst-case bound, the OnDA Algorithm provides an $\boldsymbol{x}_{n}^{\epsilon_{2}}$ decision together with its estimated certificate ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})$ for each $n$.

From the Subgradient Algorithm, we provide, for each $\Xi_{n}$, a sequence $\{\boldsymbol{x}^{(r)}\}_{r=r_{n}}^{r_{n+1}}$ that approaches an $\boldsymbol{x}_{n}^{\epsilon_{2}}$. If the new data set ${\Xi}_{n+1}$ is received before reaching $\boldsymbol{x}_{n}^{\epsilon_{2}}$, we initialize the next sub-sequence obtained by applying the Subgradient Algorithm, using the best decision at current iteration $r$, i.e., $\boldsymbol{x}^{(r_{n+1})}:=\boldsymbol{x}^{\textup{best}}_{n}\in\operatorname*{argmin}_{k\in\{{r_{n}},\dots,{r}\}}\{{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(k)})\}$. Then by connecting these sequences over $n$, our goal is achieved.

VII Data Assimilation via OnDA Algorithm

This section summarizes and analyzes our Online Data Assimilation Algorithm (OnDA Algorithm) for online data sets $\{{\Xi}_{n}\}_{n=1}^{N}$. Specifically, we present the algorithm procedure, its transient behavior and the convergence result.

The OnDA Algorithm starts from some random initial decision $\boldsymbol{x}^{(r)}\in^{d}$ and a data set ${\Xi}_{n}$, with $r=1$ and $n=1$. Then, it first generates the certificate ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})$ via C-Gen Algorithm, then it executes the Subgradient Algorithm to obtain the decisions $\{\boldsymbol{x}^{(r+1)},\boldsymbol{x}^{(r+2)},\ldots\}$ with lower and lower certificates $\{{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r+1)}),{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r+2)}),\ldots\}$. This algorithm has the anytime property, meaning that the performance guarantee is provided anytime, as soon as the first $\epsilon_{1}$-proper data-driven decision with certificate ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})$ is found. If no new data set ${\Xi}_{n+1}$ comes in, the algorithm terminates as soon as the Subgradient Algorithm terminates at iteration $r_{n+1}$. Otherwise, the algorithm resets the C-Gen Algorithm and the Subgradient Algorithm to update the decision using more data. This achieves lower certificates with higher confidence until we obtain the lowest possible certificate and guarantee the performance almost surely. The details of the whole algorithm procedure are summarized in the table of OnDA Algorithm.

The transient behavior of the OnDA Algorithm is affected by the data-streaming rate and the rate of convergence of the intermediate algorithms (decision-update rate and certificate-update rate). To further describe these effects in each time period $n$, we say that the data-streaming rate is slow with respect to the decision-update rate, if we can find an $\boldsymbol{x}_{n}^{\epsilon_{2}}$ via the OnDA Algorithm, where the worst-case scenario is described in Theorem VI.1. Further, we call it slow with respect to the certificate-update rate, if we can find at least one certificate during this time period, where the worst-case scenario is described in Theorem V.1. When the data-streaming rate is slow w.r.t. the decision-update rate for all time periods, the OnDA Algorithm guarantees to find $\{\boldsymbol{x}_{n}^{\epsilon_{2}}\}_{n=1}^{N}$. When the data-streaming rate is slow w.r.t. the decision-update rate for at least one time period $n_{0}$, it guarantees to find an $\boldsymbol{x}_{n_{0}}^{\epsilon_{2}}$. When the data-streaming rate is slow w.r.t. the certificate-update rate for at least one time period $n_{0}$, the OnDA Algorithm guarantees to find a ${J}_{n_{0}}^{\epsilon_{1}}$ for an $\boldsymbol{x}^{(r)}$ and $r\geq r_{n_{0}}$. When the data-streaming rate is not slow w.r.t. the certificate-update rate for any time period, the OnDA Algorithm will hold on the newly streamed data set, to make the data-streaming rate slow w.r.t. the decision-update rate and achieve a better data-driven decision efficiently.

Next, we state the convergence result of the OnDA Algorithm when the data streams are slow w.r.t. the decision-update rate for all time periods.

Theorem VII.1** (Finite convergence of the OnDA Algorithm).**

Consider tolerances $\epsilon_{1}$, $\epsilon_{2}>0$ and streaming data sets $\{{\Xi}_{n}\}_{n=1}^{N}$ with $N<\infty$ for a decision making problem (P). Assume that the data streams are slow w.r.t. the to decision-update rate for all $n$, i.e., assume the length of each time period $n$ is no shorter than $\varphi(1,\bar{r})$, where $\varphi$ and $\bar{r}$ are described as in Theorem VI.1 and Lemma VI.3, respectively. Then, the OnDA Algorithm guarantees to find a sequence of $\epsilon_{2}$-optimal $\epsilon_{1}$-proper data-driven decisions $\{\boldsymbol{x}_{n}^{\epsilon_{2}}\}_{n=1}^{N}$ associated with the sequence of the certificates $\{{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})\}_{n=1}^{N}$ so that the performance guarantee (2) holds for all $n$. Furthermore, the values of these certificates are guaranteed to be low in high probability. That is, for each $n$

[TABLE]

holds, where $J^{\star}:=\inf\limits_{\boldsymbol{x}\in^{d}}{\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x},\xi)]}$ is the optimal objective value for the original problem (P), the parameter $\hat{L}$ depends on steepness of the function $f$, and the parameter $\epsilon(\beta_{n})$ is determined as in Lemma 5.

In addition, given any tolerance $\epsilon_{3}$, data stream that is slow w.r.t. the decision-update rate for all $n\in\{1,\dots,N\}$ with $N\rightarrow\infty$, and $\sum_{n=1}^{\infty}\beta_{n}<\infty$, there exists a large enough number $n_{0}(\epsilon_{3})>0$, such that the algorithm terminates in finite time with a guaranteed $\epsilon_{2}$-optimal and $\epsilon_{1}$-proper data-driven decision $\boldsymbol{x}_{n_{0}}^{\epsilon_{2}}$ and a certificate ${J}_{n_{0}}^{\epsilon_{1}}(\boldsymbol{x}_{n_{0}}^{\epsilon_{2}})$ such that the performance guarantee holds almost surely. That is,

[TABLE]

and meanwhile the quality of the designed certificate ${J}_{n_{0}}^{\epsilon_{1}}(\boldsymbol{x}_{n_{0}}^{\epsilon_{2}})$ is guaranteed. In other words, for all the rest of the data sets $\{{\Xi}_{n}\}_{n=n_{0}}^{\infty}$, any element in the desired certificate sequence $\{{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})\}_{n=n_{0}}^{\infty}$ satisfies

[TABLE]

Theorem 10 quantifies the goodness of the certificates that are achievable via the OnDA Algorithm, under the condition that the data-streaming rate is slow w.r.t. the decision-update rate. Intuitively, the smaller the tolerances $\epsilon_{i}$ are, the lower the certificates become. Further, as $N\rightarrow\infty$, the smaller the parameter $\epsilon(\beta_{N})$ is and the higher the confidence $1-\beta_{N}\rightarrow 1$. When infinitely many data are streamed in, the theorem implies that we can get arbitrarily close to the optimal decision with probability one.

Remark VII.1** (Selection of tolerances $\epsilon_{2}$, $\epsilon_{1}$ and

$\epsilon_{\textup{SA}}$).**

In practice, the tolerance $\epsilon_{2}$ determines the performance bound of ${J}_{n}$, which governs the whole algorithm. With a given $\epsilon_{2}$, tolerance $\epsilon_{1}$ and $\epsilon_{\textup{SA}}$ can be chosen following the rule in Lemma VI.3. Intuitively, $\epsilon_{1}$ can be chosen to be two orders of magnitude smaller than $\epsilon_{2}$, while $\epsilon_{\textup{SA}}$ can be an order of magnitude smaller than $\epsilon_{2}$. These tolerances can also be chosen in a data-driven fashion, to achieve asymptotic convergence, or a better transient behavior of the algorithm.* *

Remark VII.2** (Asymptotic behavior of the OnDA Algorithm).**

Theorem 10 claims the convergence of the OnDA Algorithm to a decision with a desired certificate using large but a finite data set. The smaller $\epsilon_{3}$ is, the larger data set is needed to achieve the desired certificate. Because tolerances $\epsilon_{1}$,$\epsilon_{2}$ and $\epsilon_{3}$ can be chosen arbitrarily small, the certificate can indeed approach to $J^{\star}$. However, $\epsilon_{1}$,$\epsilon_{2}$ may affect the transient behavior of the algorithm and the data-streaming rate. In practice, to reach $J^{\star}$, these tolerances can be chosen in a data-driven fashion, for example diminishing sequences.* *

VIII Data incremental covering

In this section, we aim to handle large streaming data sets for efficient Online Data Assimilation Algorithm (OnDA Algorithm). To achieve this, we firstly propose an Incremental Covering Algorithm (I-Cover Algorithm). This algorithm leverages the pattern of the data points to obtain a new ambiguity set, denoted by $\mathcal{\tilde{P}}_{n}$. Then, we adapt $\mathcal{\tilde{P}}_{n}$ for a variant of the OnDA Algorithm. The resulting algorithm enables us to construct subproblems which have a lower dimension than those generated without it, and we verify its capability of handling large data sets in simulation.

VIII-A The I-Cover Algorithm

Let $\zeta$ and $\omega$ denote the center and radius of the Euclidean ball $B_{\omega}(\zeta)$, respectively. For each data set $\Xi_{n}$ and a given $\omega$, let $\mathcal{C}_{n}\subset\Xi_{n}$ denote the set of points such that $\Xi_{n}\subset\cup_{\zeta\in\mathcal{C}_{n}}B_{\omega}(\zeta)$. Let $p:=|\mathcal{C}_{n}|$ denote the number of these Euclidean balls. To account for the number of data points that are covered by a specific ball, we associate each ball $B_{\omega}(\zeta_{k})$ a weighting parameter $\theta_{k}$. We denote by $\mathcal{Q}_{n}:=\{\theta_{k}\}_{k=1}^{p}$ the set of these parameters. Then, as data sets $\{{\Xi}_{n}\}_{n=1}^{N}$ are sequentially accessible, we are to incrementally cover data sets by adapting $\mathcal{C}_{n}$ and $\mathcal{Q}_{n}$.

Formally, the I-Cover Algorithm works as follows. Let $\mathcal{C}_{0}=\O$ and $\mathcal{Q}_{0}=\O$. For the $n^{\textup{th}}$ time period with set $\Xi_{n}$, we initialize sets as $\mathcal{C}_{n}:=\mathcal{C}_{n-1}$ and $\mathcal{Q}_{n}:=\mathcal{Q}_{n-1}$. To generate a random cover for $\Xi_{n}$, we randomly and sequentially evaluate each newly streamed data point. Let $\varsigma\in\Xi_{n}\setminus\Xi_{n-1}$ denote the data point under consideration. If $\varsigma\notin B_{\omega}(\zeta_{k})$ for all $\zeta_{k}\in\mathcal{C}_{n}$, we update $\mathcal{C}_{n}\leftarrow\mathcal{C}_{n}\cup\{\zeta_{p+1}:=\varsigma\}$, $\mathcal{Q}_{n}\leftarrow\mathcal{Q}_{n}\cup\{\theta_{p+1}:=1\}$ and $p\leftarrow|\mathcal{C}_{n}|$. If $\varsigma$ is covered by some (at least one) Euclidean balls, i.e., $\varsigma\in B_{\omega}(\zeta_{k})$ for some $k$ with $\zeta_{k}\in\mathcal{C}_{n}$, we only update $\mathcal{Q}_{n}$. Let $\ell_{\varsigma}$ denote the number of the balls that cover $\varsigma$ and let $I_{\varsigma}\subset\{1,\dots,p\}$ denote the index set of these balls. Then we update elements of $\mathcal{Q}_{n}$ via $\theta_{k}\leftarrow\theta_{k}+\ell_{\varsigma}^{-1}$ for all $k\in I_{\varsigma}$. After all the new data points have been evaluated in this way, we achieve a cover of $\Xi_{n}$. Then, as the data set streams over time, the algorithm incrementally updates the cover and weights. By construction, we see that $|\mathcal{C}_{n}|\leq n$.

Next, we use $\mathcal{C}_{n}$ and $\mathcal{Q}_{n}$ to construct a new ambiguity set that results in potentially low dimensional subproblems in the OnDA Algorithm.

VIII-B Integration of I-Cover Algorithm into OnDA Algorithm

Following the I-Cover Algorithm, we consider a distribution $\tilde{\mathbb{P}}^{n}$ associated with ${\Xi}_{n}$, as follows

[TABLE]

where $\delta_{\{\zeta_{k}\}}$ is a Dirac measure at the center of the covering ball $B_{\omega}(\zeta_{k})$ and $\theta_{k}$ is the associated weight of $B_{\omega}(\zeta_{k})$. We claim the distribution $\tilde{\mathbb{P}}^{n}$ is close to the empirical distribution $\hat{\mathbb{P}}^{n}$ under the Wasserstein metric, using the following lemma.

Lemma VIII.1** (Distribution $\tilde{\mathbb{P}}^{n}$ is a good estimate of

$\hat{\mathbb{P}}^{n}$).**

*Let the radius $\omega$ of the Euclidean ball be chosen. Then the distribution $\tilde{\mathbb{P}}^{n}$ constructed by the I-Cover Algorithm on ${\Xi}_{n}$ is close to $\hat{\mathbb{P}}^{n}$ under the Wasserstein metric, i.e., $d_{W}(\hat{\mathbb{P}}^{n},\tilde{\mathbb{P}}^{n})\leq\omega$.

Then equipped with Lemma VIII.1 and Theorem IV.1 on the measure of concentration result, we can provide the certificate that ensures the performance guarantee in (1).

Lemma VIII.2** (Tractable certificate generation for $\boldsymbol{x}$ with

Performance Guarantee (1) using $\tilde{\mathbb{P}}^{n}$).**

Given ${\Xi}_{n}:=\{{\xi}_{k}\}_{k=1}^{n}$, $\beta_{n}\in(0,1)$, $\boldsymbol{x}\in^{d}$, and the radius $\omega$ of the covering balls. Define the new ambiguity set $\mathcal{\tilde{P}}_{n}:=\mathbb{B}_{\tilde{\epsilon}(\beta_{n})}(\tilde{\mathbb{P}}^{n})$ where the center of the Wasserstein ball $\tilde{\mathbb{P}}^{n}$ is defined in (11) and the radius $\tilde{\epsilon}(\beta_{n}):=\epsilon(\beta_{n})+\omega$. Then the following certificate satisfies (1) for all $\boldsymbol{x}\in^{d}$

[TABLE]

Further, under the same assumptions required in Theorem IV.2 we have the new version of (P1n) as follows

[TABLE]

and the associated worst-case distribution $\tilde{\mathbb{Q}}_{n}^{\star}(\boldsymbol{x})$ is a weighted version of $\mathbb{Q}_{n}^{\star}(\boldsymbol{x})$ in Theorem IV.2, i.e.,

[TABLE]

*where $\boldsymbol{y}^{\star}:=(\boldsymbol{y}^{\star}_{1},\ldots,\boldsymbol{y}^{\star}_{p})$ is an optimizer of (${{\tilde{\rm{P1}}_{n}}}$).

Remark VIII.1** (New version of (P2n)).**

The equivalent formulation of Problem (${{\tilde{\rm{P1}}_{n}}}$) is a new version of (P2n), defined as follows

[TABLE]

where for each $k\in\{1,\dots,p\}$, $\zeta_{k}\in\mathcal{C}_{n}$ and $\boldsymbol{x}\in^{d}$, we define ${{\tilde{h}_{k}}}:^{m}\rightarrow$ as

${{\tilde{h}_{k}}}(\boldsymbol{y}):=\theta_{k}f(\boldsymbol{x},\zeta_{k}-\boldsymbol{y})$.

With the constructed ambiguity set $\mathcal{\tilde{P}}_{n}$ and certificate function ${J}_{n}$, the the developed algorithms in Section V and Section VI are valid to solve Problem (${{\tilde{\rm{P2}}_{n,p}}}$). And the main Theorem 10 on the finite convergence of the OnDA Algorithm is valid for the certificate function ${J}_{n}$ where the only difference is that the quality of the certificate for $\boldsymbol{x}_{n}^{\epsilon_{2}}$ in (10) is replaced by

[TABLE]

where $n_{0}$ is the number of the data set in ${\Xi}_{n_{0}}$ and $p_{n_{0}}$ indicates the number of Euclidean balls that cover ${\Xi}_{n_{0}}$.

IX Simulation results

In this section, we demonstrate the application of the proposed algorithms on two case studies, with a potentially large streaming data set.

IX-A Study 1: The Effect of the I-Cover Algorithm

In order to visualize the effect of the I-Cover Algorithm, here we solve a toy problem in form of (P) using OnDA Algorithm, with and without the I-Cover Algorithm respectively. Let $x\in$ be the variable for Problem (P). Assume there are $N=200$ data points $\{{\xi}_{k}\}_{k=1}^{N}$ streaming into the algorithm. Assume each time period is one second, and for each second $k$ we only stream in one data point ${\xi}_{k}\in^{3}$, where ${\xi}_{k}$ is a realization of the unknown distribution $\mathbb{P}$. The $\mathbb{P}$ we use for simulation is a multivariate weighted Gaussian mixture distribution with three centers, where each center has mean $\mu_{1}=(2,-4,3)$, $\mu_{2}=(-3,5,0)$, $\mu_{3}=(0,0,-6)$, variance $\sum_{1}=\rm{diag}(1,3,2)$, $\sum_{2}=2\cdot\boldsymbol{I}_{3}$, $\sum_{3}=\boldsymbol{I}_{3}$, and weights 0.25, 0.5, 0.25, respectively. Let the cost function to be $f({x},\xi):=x^{2}-{\xi}^{\top}\xi$, the confidence be $1-\beta_{n}:=1-0.95e^{1-\sqrt{n}}$ and we use the parameter $c_{1}=2$, $c_{1}=1$ to design the radius $\epsilon(\beta_{n})$ of the Wasserstein ball in (4). The radius of the Euclidean ball for the I-Cover Algorithm is $w=1.5$. We sample the initial decision $x^{(0)}$ from the uniform distribution $[-10,10]$. The tolerance for the algorithm is $\epsilon_{1}=10^{-5}$, $\epsilon_{2}=10^{-4}$.

Figure 3a and Figure 3b demonstrate the effect of the I-Cover Algorithm in the OnDA Algorithm. Specifically, Figure 3a shows the incremental data covering at the end of the $200^{\textup{th}}$ time period in the $(\xi_{1},\xi_{2},\xi_{3})$ coordinates. The large shaded area are $59$ Euclidean balls $B_{\omega}$ with their centers $\{\zeta_{k}\}_{k=1}^{59}$ denoted by some of the small circles, where all these small circles constitute the streamed data set ${\Xi}_{200}:=\{{\xi}_{k}\}_{k=1}^{200}$. In Figure 3b, the gray dashed line represents the number of the data points used as centers of the empirical distribution $\hat{\mathbb{P}}^{n}$ over time and the black dashed line is that for distribution $\tilde{\mathbb{P}}^{n}$. Clearly as the data streams over time, the number $p:=|\mathcal{C}_{n}|$ is significantly smaller than $n:=|{\Xi}_{n}|$, which results in the size of Problem (${{\tilde{\rm{P2}}_{n,p}}}$) being much smaller than that of (P2n). Further, the gray solid line counts the total number of subproblems (CP${}_{n}^{(l)}$) solved to generate certificates over time and the black solid line represents that for subproblems ($\tilde{\rm{CP}}_{n}^{(l)}$) in solution to (${{\tilde{\rm{P2}}_{n,p}}}$). These subproblems search the explicit solution for the $\epsilon_{1}$-worst-case distribution and consume the major computing resources in the OnDA Algorithm. It can be seen that the number of ($\tilde{\rm{CP}}_{n}^{(l)}$) solved over time is on average only half of the (CP${}_{n}^{(l)}$) in each time period. Together, the dimension and total number of subproblems ($\tilde{\rm{CP}}_{n}^{(l)}$) solved with the I-Cover Algorithm is significantly smaller than that without it.

To evaluate the quality of the obtained $\epsilon_{1}$-proper data-driven decision with the streaming data, we estimate the optimizer of (P), $\boldsymbol{x}^{\star}$, by minimizing the average value of the cost function $f$ for a validation data set of $N_{\textup{val}}=10^{4}$ data points randomly generated from the distribution $\mathbb{P}$ (in the simulation case $\mathbb{P}$ is known). We take the resulting objective value as the estimated optimal objective value for Problem (P), i.e., $J^{\star}:={J}^{\star}({\boldsymbol{x}}^{\star})$. We calculate ${J}^{\star}({\boldsymbol{x}}^{\star})$ using the underline distribution $\mathbb{P}$, serving as the true but unknown scale to evaluate the goodness of the certificate obtained throughout the algorithm.

Figure 3c and Figure 3d show the evolution of the certificate sequence $\{{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})\}_{n=1,r=1}^{N,\infty}$ with the I-Cover Algorithm and that without the I-Cover Algorithm, respectively. Here, the optimal decision of (P) is trivially ${\boldsymbol{x}}^{\star}=0$, and for both algorithms the subgradient counterpart of the OnDA Algorithm returns the optimal decision after the first data point ${\xi}_{1}$ is used. Therefore, after a very short period within the first second, both figures start reflecting the certificate evolution under the decision sequence $\{\boldsymbol{x}^{(r)}\approx 0\}_{r=r_{2}}^{\infty}$. The gray solid line in both Figure 3c and Figure 3d show the relative goodness of the certificates for the currently used $\epsilon_{1}$-proper data-driven decision $\boldsymbol{x}^{(r)}\approx 0$ calibrated by the estimated optimal value $J^{\star}$ over time. The black segments on the gray solid line indicate that the C-Gen Algorithm is executing for certificates update, while at these time intervals the old certificate ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})$, associated with the $\epsilon_{2}$-optimal and $\epsilon_{1}$-proper data-driven decision $\boldsymbol{x}_{n}^{\epsilon_{2}}$, is still valid to guarantee the performance under the old confidence $1-\beta_{n}$. This situation commonly happens when a new data set ${\Xi}_{n+1}$ is streamed in and a new certificate ${J}_{n+1}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})$ is yet to be obtained. It can be seen that after a few samples streamed, both the obtained certificate becomes close (within $10\%$) to the estimated true optimal value $J^{\star}$. In Figure 3d however, as the data streams over $50$ seconds, the computing cost for updating certificates becomes significant for the algorithm without the I-Cover Algorithm. After $100^{\textup{th}}$ data point has been assimilated, the certificate ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})$ stops updating for all $n\geq 100$. And, further, after all the data points streamed (in $200$ seconds), the algorithm took about $70$ seconds to terminate the algorithm with certificate ${J}_{200}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})$. This is a clear disadvantage compared to the algorithm with the I-Cover Algorithm, which terminates as soon as all the data points were taken in.

IX-B Study 2: OnDA Algorithm with Large Streaming Data Sets

Here, we are to find an $\epsilon_{2}$-optimal, $\epsilon_{1}$-proper decision $\boldsymbol{x}\in^{30}$ for Problem (P). We consider $N=500$ iid sample points $\{{\xi}_{k}\}_{k=1}^{N}$ streaming randomly in between every $1$ to $3$ seconds with each data point ${\xi}_{k}\in^{10}$ a realization of $\mathbb{P}$. We assume that the unknown distribution is a multivariate Gaussian mixture distribution with three centers where the components of the mean of each center are uniformly chosen between $[-10,10]$, and the variance matrix is $\boldsymbol{I}_{m}$ for each center. We assume the cost function $f:^{30}\times^{10}\rightarrow$ to be $f(\boldsymbol{x},\xi):={\boldsymbol{x}}^{\top}A\boldsymbol{x}+{\boldsymbol{x}}^{\top}B\xi+{\xi}^{\top}C\xi$ with random values for the positive semi-definite matrix $A\in^{30\times 30}$, $B\in^{30\times 10}$ and negative definite matrix $C\in^{10\times 10}$. The radius of the Euclidean ball for the I-Cover Algorithm is $w=5$.

Similarly to Figure 3b, Figure 3e demonstrates the incremental construction of the distribution $\tilde{\mathbb{P}}^{n}$ and the accumulated number of Problem (${{\tilde{\rm{P2}}_{n,p}}}$) solved over time. Clearly, after certain amount of data have been assimilated, the structure of the data set was inferred by the I-Cover Algorithm and the number of Euclidean balls used to cover the data set is about $20$. Also, after the $100^{\textup{th}}$ time period (from $100$ to $200$ seconds in this case), the algorithm can validate new certificate without solving any Problem (${{\tilde{\rm{P2}}_{n,p}}}$). This feature dramatically improves the performance of the OnDA Algorithm and makes the algorithm flexible for online settings.

Similarly to Figure 3c, Figure 3f shows the evolution of the certificate sequence $\{{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})\}_{n=1,r=1}^{N,\infty}$ for the decision sequence $\{\boldsymbol{x}^{(r)}\}_{r=1}^{\infty}$. In the same way as in the last case study, the obtained certificate becomes close to the estimated true optimal value $J^{\star}$ (within $10\%$) after about $25$ seconds with the assimilation of $10$ data sets. Also, as more data sets are assimilated, the update of the certificate ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{\epsilon_{2}}_{n})$ remains fast and the algorithm terminates within a second after the last data set was streamed in.

X Conclusions

In this paper, we have proposed the Online Data Assimilation Algorithm (the OnDA Algorithm) to solve the problem in the form of (P), where the realizations of the unknown distribution (i.e., the streaming data) are collected over time in order for the real-time data-driven decision of (P) to have guaranteed out-of-sample performance. The data-driven decision with the certificate that guarantees out-of-sample performance are available any time during the execution of the algorithm, and the optimal data-driven decision are approached with a (sub)linear convergence rate. The algorithm terminates after collecting a sufficient amount of data to make good decision. To facilitate the decision making, an enhanced version of the proposed algorithm is further constructed, by using an Incremental Covering Algorithm (the I-Cover Algorithm) to estimate new ambiguity sets over time. We provided sample problems and showed the actual performance of the proposed OnDA Algorithm with the I-Cover Algorithm over time. Future work will generalize the results for weaker assumptions of the problem and potentially extend the algorithm to scenarios that include system dynamics.

There are mainly two types of Numerical methods that serve as the main ingredients of our OnDA Algorithm. One type is given by Frank-Wolfe Algorithm (FWA) variants and another is the Subgradient Algorithm. In this Section, we describe FWA and the Away-step Frank-Wolfe Algorithm (AFWA) for the sake of completeness. We combine AFWA with another variant, the Simplicial Algorithm, in Section V. For the Subgradient Algorithm, please refer to [17, 18, 20].

Frank-Wolfe Algorithm over a unit simplex

To solve convex programs over a unit simplex, we introduce the FWA and AFWA following [13, 12]. Let us denote the $m$-dimensional unit simplex by $\Delta_{m}:=\{\lambda\in^{m}\;|\;{\boldsymbol{1}_{m}}^{\top}\lambda=1,\;\lambda\geq 0\}$. Let ${\Lambda}_{m}$ be the set of all extreme points for the simplex $\Delta_{m}$. Consider the maximization of a concave function $f(\boldsymbol{x})$ subject to $\boldsymbol{x}\in\Delta_{m}$; we refer to this problem by $(\star)$ and denote by $\boldsymbol{x}^{\star}$ an optimizer of $(\star)$. We call ${\boldsymbol{x}}^{\epsilon}$ an $\epsilon$-optimal solution of $(\star)$, if ${\boldsymbol{x}}^{\epsilon}\in\Delta_{m}$ and $f(\boldsymbol{x}^{\star})-f({\boldsymbol{x}}^{\epsilon})\leq\epsilon$. The classical FWA solves problem $(\star)$ to an ${\boldsymbol{x}}^{\epsilon}$ via the iterative process as follows. Let $\boldsymbol{x}^{(0)}\in\Delta_{m}$ denote a random initial point for FWA. For each iteration $k$ with an $\boldsymbol{x}^{(k)}\in\Delta_{m}$, the concavity of $f$ enables $f(\boldsymbol{x}^{\star})\leq f(\boldsymbol{x}^{(k)})+{\nabla f(\boldsymbol{x}^{(k)})}^{\top}(\boldsymbol{x}^{\star}-\boldsymbol{x}^{(k)})$, which implies $f(\boldsymbol{x}^{\star})\leq f(\boldsymbol{x}^{(k)})+\max_{\boldsymbol{x}\in\Delta_{m}}{\nabla f(\boldsymbol{x}^{(k)})}^{\top}(\boldsymbol{x}-\boldsymbol{x}^{(k)})$. Using this property, we define a FW search point $\boldsymbol{s}^{(k)}$ by an extreme point such that $\boldsymbol{s}^{(k)}\in\operatorname*{argmax}_{\boldsymbol{x}\in\Delta_{m}}{\nabla f(\boldsymbol{x}^{(k)})}^{\top}(\boldsymbol{x}-\boldsymbol{x}^{(k)})$. With this search point we define the FW direction at $\boldsymbol{x}^{(k)}$ by $d_{\textup{FW}}^{(k)}:=\boldsymbol{s}^{(k)}-\boldsymbol{x}^{(k)}$. The classical FWA then iteratively finds a FW direction and solves a line search problem over this direction until an ${\epsilon}$-optimal solution ${\boldsymbol{x}}^{\epsilon}:=\boldsymbol{x}^{(k)}$ is found, certified by $\eta^{(k)}:={\nabla f(\boldsymbol{x}^{(k)})}^{\top}d_{\textup{FW}}^{(k)}\leq\epsilon$.

It is known that the classical FWA has linear convergence rate if the cost function $f$ is $\mu$-strongly concave and the optimum is achieved in the relative interior of the feasible set $\Delta_{m}$. If the optimal solution lies on the boundary of $\Delta_{m}$, then this algorithm only has a sublinear convergence rate, due to a zig-zagging phenomenon [13]. AFWA is an extension of the FWA that guarantees the linear convergence rate of the problem $(\star)$ under some conditions related to the local strong concavity. The main difference between AFWA and the classical FWA is that the latter solves the line-search problem after obtaining a ascent direction by considering all extreme points, while the AFWA chooses a ascend direction that prevents zig-zagging. We summarize the convergence properties of the AFWA here. For complete descriptions of the AFWA, we refer the reader to [13, 21]. The detailed FWA and AFWA are shown in the following Algorithm tables.

Theorem .1** (Linear convergence of AFWA [13, Theorem 8]).**

*Suppose the function $f$ has a curvature constant $C_{f}$ and a geometric strong concavity constant $\mu_{f}$ on $\Delta_{m}$, as defined in footnote 2 Let us define the decay rate $\kappa:=1-{\mu_{f}}/({4C_{f}})\in(0,1)\subset$. Then the suboptimality bound at the iteration point $\boldsymbol{x}^{(k)}$ of the AFWA decreases geometrically as $f(\boldsymbol{x}^{\star})-f(\boldsymbol{x}^{(k+1)})\leq\kappa(f(\boldsymbol{x}^{\star})-f(\boldsymbol{x}^{(k)}))$.

Proofs

Lemma 5

Proof:

Following [5, 3] and from Theorem IV.1, we prove that ${J}_{n}(\boldsymbol{x})$ is a valid certificate for (1). Knowing that (4) is obtained by letting the right-hand side of (3) to be equal to a given $\beta_{n}$, for each $n$ we substitute (4) into the right-hand side of (3), yielding $\mathbf{P}^{n}\{d_{W}(\mathbb{P},\hat{\mathbb{P}}^{n})\geq\epsilon(\beta_{n})\}\leq\beta_{n}$ for each $n$. This means that a data set ${\Xi}_{n}$ we can construct an empirical probability measure $\hat{\mathbb{P}}^{n}$ such that $d_{W}(\mathbb{P},\hat{\mathbb{P}}^{n})\leq\epsilon(\beta_{n})$ with probability at least $1-\beta_{n}$. Namely, $\mathbf{P}^{n}\{\mathbb{P}\in\mathbb{B}_{\epsilon(\beta_{n})}(\hat{\mathbb{P}}^{n})\}\geq 1-\beta_{n}$. Thus, for all $\boldsymbol{x}\in^{d}$, we have $\mathbf{P}^{n}\{\mathbb{P}\in\mathbb{B}_{\epsilon(\beta_{n})}(\hat{\mathbb{P}}^{n})\}=\mathbf{P}^{n}\{\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x},\xi)]\leq\sup\limits_{\mathbb{Q}\in\mathcal{P}_{n}}{\mathbb{E}_{\mathbb{Q}}[f(\boldsymbol{x},\xi)]}\}=\mathbf{P}^{n}\{\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x},\xi)]\leq{J}_{n}(\boldsymbol{x})\}\geq 1-\beta_{n}$.

Lemma IV.2

In the following, we for matrices $A_{1}\in^{m\times d}$ and $A_{2}\in^{p\times q}$, we let $A_{1}\oplus A_{2}$ denote their direct sum. The shorthand notation $\oplus_{i=1}^{m}A_{i}$ represents $A_{1}\oplus\cdots\oplus A_{m}$.

Proof:

To prove 1, for any feasible solution $(\tilde{\boldsymbol{u}},\tilde{\boldsymbol{v}})$ of (P2n), we compute $\frac{1}{n}\sum_{k=1}^{n}\|\tilde{\boldsymbol{y}}_{k}\|=\frac{1}{n}\sum_{k=1}^{n}\|\tilde{\boldsymbol{u}}_{k}-\tilde{\boldsymbol{v}}_{k}\|\leq\frac{1}{n}\sum_{k=1}^{n}\|\tilde{\boldsymbol{u}}_{k}\|+\frac{1}{n}\sum_{k=1}^{n}\|\tilde{\boldsymbol{v}}_{k}\|=\frac{1}{n}{\boldsymbol{1}_{mn}}^{\top}\tilde{\boldsymbol{u}}+\frac{1}{n}{\boldsymbol{1}_{mn}}^{\top}\tilde{\boldsymbol{v}}=\frac{1}{n}{\boldsymbol{1}_{2mn}}^{\top}(\tilde{\boldsymbol{u}},\tilde{\boldsymbol{v}})=\epsilon(\beta_{n})$. Therefore $(\tilde{\boldsymbol{y}}_{1},\ldots,\tilde{\boldsymbol{y}}_{n})$ is feasible for (P1n).

For 2, we exploit that any feasible solution $\tilde{\boldsymbol{y}}$ of (P1n) is a linear combination of the extreme points of the constraint set in (P1n). Let us denote the matrix $A_{n}:=[\oplus_{i=1}^{n}I_{m},-\oplus_{i=1}^{n}I_{m}]\in^{mn\times 2mn}$. By construction of Problem (P2n), we see that each column vector of the matrix $n\epsilon(\beta_{n})A_{n}$ is a concatenated vector of an extreme point of Problem (P1n), and that all the extreme points of (P1n) are included. Then, any feasible solution of (P1n) can be written as $\tilde{\boldsymbol{y}}=n\epsilon(\beta_{n})A_{n}(\hat{\boldsymbol{u}},\hat{\boldsymbol{v}})$ where $(\hat{\boldsymbol{u}},\hat{\boldsymbol{v}})$ is a vector of the convex combination coefficients of the extreme points of the constraint set in (P1n). Clearly, we have $(\hat{\boldsymbol{u}},\hat{\boldsymbol{v}})\in\Delta_{2mn}$, i.e., $n\epsilon(\beta_{n})(\hat{\boldsymbol{u}},\hat{\boldsymbol{v}})$ is in the feasible set of the Problem (P2n). Then, by construction $(\tilde{\boldsymbol{u}},\tilde{\boldsymbol{v}}):=n\epsilon(\beta_{n})(\hat{\boldsymbol{u}},\hat{\boldsymbol{v}})$ is feasible for (P2n).

For 3, since (P1n) and (P2n) are the same in the sense of (1) and (2), then if $(\tilde{\boldsymbol{u}}^{\star},\tilde{\boldsymbol{v}}^{\star})$ is an optimizer of (P2n), by letting $\tilde{\boldsymbol{y}}_{k}^{\star}:=\tilde{\boldsymbol{u}}_{k}^{\star}-\tilde{\boldsymbol{v}}_{k}^{\star}$ for each $k\in\{1,\dots,n\}$ we know the objective values of the two problems coincide. We claim that the optimum of (P1n) is achieved via the optimizer $\tilde{\boldsymbol{y}}^{\star}$. If not, then there exists $\hat{\boldsymbol{y}}^{\star}\neq\tilde{\boldsymbol{y}}^{\star}$ such that the optimum is achieved with higher value. Then, from the construction in (2) we can find a feasible solution $(\hat{\boldsymbol{u}},\hat{\boldsymbol{v}})$ of (P2n) that results in a higher objective value. This contradicts the assumption that $(\tilde{\boldsymbol{u}}^{\star},\tilde{\boldsymbol{v}}^{\star})$ is an optimizer of (P2n).

Theorem V.1

Proof:

Given tolerance $\epsilon_{1}$, decision $\boldsymbol{x}$ and any data set ${\Xi}_{n}$ with $n\in\{1,\dots,N\}$, let $H_{n}:^{mn}\rightarrow$, $H_{n}:=\frac{1}{n}\sum_{k=1}^{n}h_{k}$ denote the objective function of (P2n) and let $\mathcal{S}_{n}$ denote the family of subsets of $\Lambda_{2mn}$. In the procedure of C-Gen Algorithm, let us consider a sequence of generated candidate vertex sets: $I_{n}^{(l)}\subset I_{n}^{(l+1)}$, $l=0,1,2,\ldots$ with $I_{n}^{(l)}\in\mathcal{S}_{n}$. We show the convergence of C-Gen Algorithm for any data set ${\Xi}_{n}$, by two steps.

Step 1) The sequence $\{I_{n}^{(l)}\}_{l}$ is finite and the number of iterations is at most $2mn$. For each $l$ and candidate optimizer $\boldsymbol{y}^{(l-1)}$, we generate a nonempty set of search points $\Omega^{(l)}$ with suboptimality gap $\eta^{(l)}$ via (LP${}_{n}^{(l)}$). If $\eta^{(l)}\leq\epsilon_{1}$, then we solved (P2n) to $\epsilon_{1}$-optimality and $l$ is therefore finite, otherwise we update $I_{n}^{(l)}:=I_{n}^{(l-1)}\cup\Omega^{(l)}$. Given that the maximal cardinality of each $I_{n}^{(l)}\in\mathcal{S}_{n}$ is bounded by $2mn$, then it is sufficient to show $\Omega^{(l)}\cap I_{n}^{(l-1)}=\O$. Because $\boldsymbol{y}^{(l-1)}$ is an $\epsilon_{1}$-optimal of (CP${}_{n}^{(l)}$) under ${\operatorname{conv}}(I_{n}^{(l-1)})$, then for any ${\boldsymbol{y}}\in{\operatorname{conv}}(I_{n}^{(l-1)})$, it holds that $\frac{1}{n}\sum_{k=1}^{n}\langle\nabla h_{k}(\boldsymbol{y}^{(l-1)}_{k}),\boldsymbol{y}_{k}-\boldsymbol{y}^{(l-1)}_{k}\rangle\leq\epsilon_{1}$. Since any element in $\Omega^{(l)}$ is such that $\eta^{(l)}>\epsilon_{1}$, then for any ${\boldsymbol{y}}\in{\operatorname{conv}}(I_{n}^{(l-1)})$, we have ${\boldsymbol{y}}\notin\Omega^{(l)}$, which concludes $\Omega^{(l)}\cap I_{n}^{(l-1)}=\O$. Further, the cardinality of $\Omega^{(l)}$ is at least one for every iteration $l$, then after at most $2mn$ steps the cardinality of $I_{n}^{(l)}$ becomes $2mn$, which implies the $\epsilon_{1}$-optimality of (P2n) by the $\epsilon_{1}$-optimality of (CP${}_{n}^{(l)}$).

Step 2) The computational bound of C-Gen Algorithm is quantified. To see this, consider the problems $\{\eqref{eq:LP}\}_{l}$ and $\{\eqref{eq:CP}\}_{l}$. By Assumption V-B on the cheap access of the gradients, the computation of (LP${}_{n}^{(l)}$) is negligible. Thus, the computational bound is given by the sum of the steps to solve the $\{\eqref{eq:CP}\}_{l}$, where the number of iterations $l$ is $2mn$ in the worst case.

For each (CP${}_{n}^{(l)}$) solved by AFWA, index the AFWA iterations by $i=0,1,2,\ldots$, let ${\operatorname{obj}}_{i}^{(l)}$ be the objective value at each iteration, and assume the optimal objective value is ${\operatorname{obj}}_{\star}^{(l)}$. As in Theorem .1, let $\kappa_{{{n}},l}\in(0,1)\subset$ be the decay parameter related to local strong concavity of $H_{n}$ over $\operatorname{conv}(I_{n}^{(l)})$. Then using the linear convergence rate of the AFWA, each (CP${}_{n}^{(l)}$) achieves the following computational bound

[TABLE]

where the initial condition ${\operatorname{obj}}_{0}^{(l)}$ results from an $\epsilon_{1}$-optimal optimizer of CP at iteration $l-1$, i.e., we can equivalently denote ${\operatorname{obj}}_{0}^{(l)}$ by ${J}_{n}^{(l-1)}(\boldsymbol{x})$, for all $l\in\{1,\dots,2mn\}$.

Let us consider sequence $\{\eqref{eq:CP}\}_{l}$ with feasible sets $\{\operatorname{conv}(I_{n}^{(l)})\}_{l}$. Then we have

$\operatorname{conv}(I_{n}^{(0)})\subset\operatorname{conv}(I_{n}^{(1)})\subset\cdots\subset\operatorname{conv}(I_{n}^{(2mn)})$.

This results into monotonically decaying parameters and ($\epsilon_{1}$-)optimal objective values, as given in the following

$0<\kappa_{{{n}},1}\leq\kappa_{{{n}},2}\leq\cdots\leq\kappa_{{{n}},2mn}<1$,

${J}_{n}^{(0)}(\boldsymbol{x})\leq{J}_{n}^{(1)}(\boldsymbol{x})\leq\cdots\leq{J}_{n}^{(2mn)}(\boldsymbol{x})$,

${\operatorname{obj}}^{(0)}_{\star}\leq{\operatorname{obj}}_{\star}^{(1)}\leq\cdots\leq{\operatorname{obj}}_{\star}^{(2mn)}$.

Using the previous notation, we can identify ${J}_{n}^{(2mn)}(\boldsymbol{x})\equiv{J}_{n}^{\epsilon_{1}}(\boldsymbol{x})$, ${\operatorname{obj}}^{(0)}_{\star}\equiv{J}_{n}^{(0)}(\boldsymbol{x})$, and ${\operatorname{obj}}_{\star}^{(2mn)}\equiv{J}_{n}(\boldsymbol{x})$. Let us denote $\kappa:=\max_{{{n}},l}\{\kappa_{{{n}},l}\}$. Then, by solving each $\eqref{eq:CP}$ to $\epsilon_{1}$-optimality, it leads to the accumulated computational steps $\phi(n):=\sum_{l}i_{l}$, where each $i_{l}$ is the computation step for $\epsilon_{1}$-optimal (CP${}_{n}^{(l)}$) that satisfies the following inequality

[TABLE]

Finally, in the worst-case scenario, the computational bound of the C-Gen Algorithm is

[TABLE]

Next, we show the convergence of the C-Gen Algorithm under online data sets $\{{\Xi}_{n}\}_{n=1}^{N}$. Similarly to the proof for the computational bound for a given $n$, we can compute the worst-case bound under $\{{\Xi}_{n}\}_{n=1}^{N}$, by summing over the steps required to solve the $\{\eqref{eq:CP}\}_{n,l}$. This leads to the stated bound $\bar{\phi}(n)$, where the empirical cost $J^{\textup{sae}}_{N}(\boldsymbol{x})$ serves as the cost of initial condition $\boldsymbol{y}^{(0)}:=\boldsymbol{0}_{2mN}$. In this way, when the data-streaming rate is slower or equal than $(\bar{\phi}(1))^{-1}$, we claim that C-Gen Algorithm can always find the certificate for each data set ${\Xi}_{n}$. This is because in each time period $n$, we only have $2mn$ extreme points, and $2m(n-1)$ has been explored due to the adaptation of the candidate vertex set $I_{n}^{(0)}$.

Lemma VI.1

Proof:

For any $\boldsymbol{x}$, $\boldsymbol{y}\in^{d}$ and $t\in[0,1]\subset$, we have $\boldsymbol{z}=t\boldsymbol{x}+(1-t)\boldsymbol{y}\in^{d}$ and an optimizer of (5), $\mathbb{Q}_{n}^{\star}(\boldsymbol{z})$, such that

[TABLE]

Lemma VI.2

Proof:

Let us consider the function ${\mathbb{E}_{\mathbb{Q}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})}[f(\boldsymbol{x},\xi)]}$. Using Assumption VI on convexity of $f$ in $\boldsymbol{x}$, we have for any $\boldsymbol{z}\in{\operatorname{dom}}\>{J}_{n}$ the following relation

[TABLE]

Knowing that ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})\geq{J}_{n}(\boldsymbol{x}^{(r)})-\epsilon_{1}$ and ${J}_{n}(\boldsymbol{z})\geq{\mathbb{E}_{\mathbb{Q}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})}[f(\boldsymbol{z},\xi)]}$, this concludes the first part of the proof.

To show the second part, similarly, we also have for any $\tilde{\boldsymbol{x}}$, $\boldsymbol{z}\in{\operatorname{dom}}\>{J}_{n}$ the following relation

[TABLE]

Using Point Search Algorithm, we achieve an $\eta>0$ such that ${\mathbb{E}_{\mathbb{Q}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r)})}[f(\tilde{\boldsymbol{x}},\xi)]}\geq{J}_{n}(\tilde{\boldsymbol{x}})-\eta$. Finally, by similar statement as in the first part, we claim $g_{n}^{r}(\tilde{\boldsymbol{x}})\in\partial_{\epsilon}{J}_{n}(\tilde{\boldsymbol{x}})$.

Lemma VI.3

Proof:

In the $n^{\textup{th}}$ time period, let us consider subgradient iterates $i$ for all $r_{n}\leq i\leq r$

[TABLE]

From Lemma VI.2, we know that ${J}_{n}(\boldsymbol{x}_{n}^{\star})\geq{J}_{n}(\boldsymbol{x}^{(i)})+{g_{n}^{i}(\boldsymbol{x}^{(i)})}^{\top}(\boldsymbol{x}_{n}^{\star}-\boldsymbol{x}^{(i)})-\epsilon_{\textup{SA}}$ for all $\boldsymbol{x}^{(i)}$. Then, we have

[TABLE]

Combining the inequalities over iterations from $r_{n}$ to $r$ gives

[TABLE]

Then, using the fact that

[TABLE]

and the previous iteration, we have

[TABLE]

Next, it remains to select a step-size rule $\{\alpha^{(i)}\}_{i=r_{n}}^{r}$ such that 1) the above right hand side term is upper bounded by $\epsilon_{2}$, and 2) the number of subgradient iterations $\bar{r}$ as described in the lemma is bounded. Note that the selection procedure is not unique, so we propose two step-size rules to obtain an explicit expression of $\bar{r}$.

For any data set $\Xi_{n}$, let us select a sufficiently large value $M$ to be the diameter of the decision domain of interest, i.e.,

[TABLE]

Then the step size rule and $\bar{r}$ has to satisfy the following

[TABLE]

We first consider a constant step-size rule, and select the step size as follows

[TABLE]

Then, to satisfy (13), we determine

[TABLE]

Alternatively, consider the divergent but square-summable step-size rule, i.e. $\sum_{i=r_{n}}^{\infty}{\alpha^{(i)}}=\infty$, $\sum_{i=r_{n}}^{\infty}{(\alpha^{(i)})^{2}}<\infty$. For this class of step-size rules, as $\bar{r}$ increases to $\infty$, we have the left-hand side term of (13) goes to $\mu\epsilon_{\textup{SA}}<\epsilon_{2}$, then there exists a large enough but finite number $\bar{r}$, such that (13) holds. To see this explicitly, we select the step-size rule to be the harmonic sequence as follows

[TABLE]

Now we upper bound the numerator and lower bound the denominator of (13) using the following fact

[TABLE]

Then we determine $\bar{r}$ to be the following

[TABLE]

This concludes the proof.

Theorem VI.1

Proof:

The computational bound to achieve an $\boldsymbol{x}_{n}^{\epsilon_{2}}$ strongly depends on the subgradient iterations $\bar{r}:=r_{n+1}-r_{n}$ in Lemma VI.3 and the number of subgradient functions $\{g_{n}^{r}\}_{r=r_{n}}^{r_{n+1}}$ constructed via the C-Gen Algorithm. To characterize this bound, we quantify the computational steps for $\{g_{n}^{r}\}_{r=r_{n}}^{r_{n+1}}$ next.

For each time period $n$, let us assume the C-Gen Algorithm has explored the feasible set of (P2n) when obtaining the initial certificate ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r_{n})})$. This procedure consumes a worst-case computational time $\phi(n)$, (or $\bar{\phi}(1)$ if a data-streaming scenario), as stated in Theorem V.1. After this initial step, every time the Subgradient Algorithm needs to execute C-Gen Algorithm at some $r\leq r_{n+1}$, C-Gen Algorithm will solve a unique (CP${}_{n}^{(l)}$) and return $\mathbb{Q}_{n}^{\epsilon_{1}}$ for an ${\epsilon}_{1}$-subgradient function $g_{n}^{r}$ at $\boldsymbol{x}^{(r)}$. Let ${\textrm{CP}}_{r}$ denote the unique (CP${}_{n}^{(l)}$) solved at $\boldsymbol{x}^{(r)}$. Then, to quantify the computational steps for $\{g_{n}^{r}\}_{r=r_{n}}^{r_{n+1}}$, we compute the sum of the steps to solve $\{{\textrm{CP}}_{r}\}_{r}$.

Let us denote the number of steps solving ${\textrm{CP}}_{r}$ by $i_{r}$, for all $r\in\{r_{n},\ldots,r_{n+1}\}$. Then, we aim to quantify $i_{r+1}$ for $g_{n}^{r+1}$. To achieve this, let us assume a subgradient function $g_{n}^{r}$ is computed at an iteration $r$. Then we perform a subgradient iteration (7) and obtain an $\boldsymbol{x}^{(r+1)}$. By using a subgradient estimation strategy, we obtain the optimality gap $\eta^{(1)}$ via Point Search Algorithm, denoted by $\bar{\eta}_{r+1}:=\eta^{(1)}$. This gap $\bar{\eta}_{r+1}$ enables us to quantify the distance between the initial objective value and the optimal objective value of ${\textrm{CP}}_{r+1}$. When $\bar{\eta}_{r+1}\leq\epsilon_{\textup{SA}}$, the algorithm uses the estimated subgradient function and $i_{r+1}=0$. Otherwise, the computational steps can be calculated via convergence of AFWA for ${\textrm{CP}}_{r+1}$, by $\kappa^{i_{r+1}}\bar{\eta}_{r+1}\leq\epsilon_{1},$ where $\kappa$, or using $\bar{\kappa}$ for the data-streaming case, is determined as in Theorem V.1. Let us consider a threshold value $t_{r}$

[TABLE]

Then we can represent each value $i_{r}$ by $i_{r}={\log}_{\kappa}(\frac{\epsilon_{1}}{t_{r}}),\;r\in\{r_{n},\ldots,r_{n+1}\}.$ Let us denote $t:=\max_{r}\{t_{r}\}$. Then, the computational steps for $\{g_{n}^{r}\}_{r=r_{n}}^{r_{n+1}}$, $\sum_{r}i_{r}$, are bounded by $\sum_{r}i_{r}\leq\bar{r}{\log}_{\kappa}(\frac{\epsilon_{1}}{t}).$ Finally, the computational steps to achieve an $\boldsymbol{x}_{n}^{\epsilon_{2}}$, denoted by $\varphi(n,\bar{r}):=\phi(n)+\sum_{r}i_{r}+\bar{r}$, are bounded as $\varphi(n,\bar{r})\leq\phi(n)+\bar{r}\left({\log}_{\kappa}(\frac{\epsilon_{1}}{t})+1\right)$. Again, one should use $\bar{\phi}(1)$ in the bound in place of $\phi(n)$ if considering a data-streaming scenario.

Theorem 10

Proof:

The first part of the proof is an application of Theorem V.1 and Theorem VI.1. For any data set ${\Xi}_{n}$ and the initial data-driven decision $\boldsymbol{x}^{(r_{n})}$, by Theorem V.1 we can show $\boldsymbol{x}^{(r_{n})}$ to be $\epsilon_{1}$-proper, via finding ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r_{n})})$ such that $\mathbf{P}^{n}({\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x}^{(r_{n})},\xi)]}\leq{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{(r_{n})})+\epsilon_{1})\geq 1-\beta_{n}$. Then using Theorem VI.1, an $\epsilon_{2}$-optimal $\epsilon_{1}$-proper data-driven decision $\boldsymbol{x}^{\epsilon_{2}}_{n}$ with certificate ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})$ can be achieved. Therefore the performance guarantee (2) holds for $\boldsymbol{x}^{\epsilon_{2}}_{n}$, i.e., $\mathbf{P}^{n}({\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x}_{n}^{\epsilon_{2}},\xi)]}\leq{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})+\epsilon_{1})\geq 1-\beta_{n}$.

In the following, we show the certificate ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})$ can be upper bounded in high probability, for each $n$.

First, let $\boldsymbol{x}^{\delta}$ denote the $\delta$-optimal solution of (P), i.e., ${\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x}^{\delta},\xi)]}\leq J^{\star}+\delta$. By construction of the certificate in the algorithm we have ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})\leq{J}_{n}(\boldsymbol{x}_{n}^{\epsilon_{2}})\leq{J}_{n}(\boldsymbol{x}_{n}^{\star})+\epsilon_{2}\leq{J}_{n}(\boldsymbol{x}^{\delta})+\epsilon_{2}\leq{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{\delta})+\epsilon_{1}+\epsilon_{2}$ for all $n$, where the first inequality holds because ${J}_{n}$ is the function that achieves the supreme of Problem (5) while ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})$ is the objective value for a feasible distribution $\mathbb{Q}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{\epsilon_{2}})$, the second inequality holds because $\boldsymbol{x}_{n}^{\epsilon_{2}}$ is $\epsilon_{2}$-optimal, the third inequality holds because $\boldsymbol{x}_{n}^{\star}$ is a minimizer of the certificate function ${J}_{n}$, the last inequality holds because the C-Gen Algorithm for certificate generation guarantees the existence of ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{\delta})$ such that ${J}_{n}(\boldsymbol{x}^{\delta})\leq{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{\delta})+\epsilon_{1}$, with an distribution $\mathbb{Q}_{n}^{\epsilon_{1}}({\boldsymbol{x}}^{\delta})$ satisfying $d_{W}(\hat{\mathbb{P}}^{n},\mathbb{Q}_{n}^{\epsilon_{1}}({\boldsymbol{x}}^{\delta}))\leq\epsilon(\beta_{n})$.

Next, we exploit the connection between ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{\delta})$ and $J^{\star}$. By Assumption 2 on the concavity of $f$ in $\xi$, there exists a constant $\hat{L}>0$ such that $f(x,\xi)\leq\hat{L}(1+\|\xi\|)$ holds for all $x\in^{d}$ and $\xi\in\mathcal{Z}$. Then by the dual representation of the Wasserstein metric from Kantorovich and Rubinstein [15, 3] we have ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{\delta}):={\mathbb{E}_{\mathbb{Q}_{n}^{\epsilon_{1}}({\boldsymbol{x}}^{\delta})}[f(\boldsymbol{x}^{\delta},\xi)]}\leq{\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x}^{\delta},\xi)]}+\hat{L}d_{W}(\mathbb{P},\mathbb{Q}_{n}^{\epsilon_{1}}({\boldsymbol{x}}^{\delta}))$. In order to quantify the last term, we apply the triangle inequality, which gives us $d_{W}(\mathbb{P},\mathbb{Q}_{n}^{\epsilon_{1}}({\boldsymbol{x}}^{\delta}))\leq d_{W}(\mathbb{P},\hat{\mathbb{P}}^{n})+d_{W}(\hat{\mathbb{P}}^{n},\mathbb{Q}_{n}^{\epsilon_{1}}({\boldsymbol{x}}^{\delta}))$. Then by the performance guarantee we have $\mathbf{P}^{n}\{d_{W}(\mathbb{P},\hat{\mathbb{P}}^{n})\leq\epsilon(\beta_{n})\}\geq 1-\beta_{n}$, and by the the way of constructing $\mathbb{Q}_{n}^{\epsilon_{1}}({\boldsymbol{x}}^{\delta})$ we have $d_{W}(\hat{\mathbb{P}}^{n},\mathbb{Q}_{n}^{\epsilon_{1}}({\boldsymbol{x}}^{\delta}))\leq\epsilon(\beta_{n})$. These inequalities result in $\mathbf{P}^{n}\{d_{W}(\mathbb{P},\mathbb{Q}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{\delta}))\leq 2\epsilon(\beta_{n})\}\geq 1-\beta_{n}$. We use now this bound to deal with the last term in the upper bound of ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{\delta})$. In particular, we have $\mathbf{P}^{n}\{{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{\delta})\leq{\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x}^{\delta},\xi)]}+2\hat{L}\epsilon(\beta_{n})\}\geq 1-\beta_{n}$ for all $n$. Using the obtained inequality ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})\leq{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}^{\delta})+\epsilon_{1}+\epsilon_{2}$ and knowing $\delta$ can be arbitrary small, we achieved the goal as in (8).

Now, it remains to find an $n_{0}$, associated with an $\epsilon_{2}$-optimal and $\epsilon_{1}$-proper data-driven decision $\boldsymbol{x}_{n_{0}}^{\epsilon_{2}}$, such that the almost sure guarantee (9) and bound (10) of the certificate ${J}_{n_{0}}^{\epsilon_{1}}(\boldsymbol{x}_{n_{0}}^{\epsilon_{2}})$ can be guaranteed for the termination of the OnDA Algorithm as $N\rightarrow\infty$. We achieve this by two steps.

First, we show the almost sure performance guarantee when the data set is sufficiently large. For any time period $n$, the algorithm finds $\boldsymbol{x}^{\epsilon_{2}}_{n}$ with the performance guarantee (2), which can be equivalently written as $\mathbf{P}^{n}({\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x}_{n}^{\epsilon_{2}},\xi)]}\geq{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})+\epsilon_{1})\leq\beta_{n}$. As $\sum_{n=1}^{\infty}\beta_{n}<\infty$, from the $1^{\textup{st}}$ Borel-Cantelli Lemma we have that $\mathbf{P}^{\infty}\{{\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x}_{n}^{\epsilon_{2}},\xi)]}\geq{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})+\epsilon_{1}\;\textrm{occurs infinitely many often}\}=0$. That is, almost surely we have that ${\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x}_{n}^{\epsilon_{2}},\xi)]}\geq{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})+\epsilon_{1}$ occurs at most for finite number of $n$. Thus, there exists a sufficiently large $n_{1}$, such that for all $n\geq n_{1}$, we have ${\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x}_{n}^{\epsilon_{2}},\xi)]}\leq{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})+\epsilon_{1}$ occurs almost surely, i.e., $\mathbf{P}^{n}({\mathbb{E}_{\mathbb{P}}[f(\boldsymbol{x}_{n}^{\epsilon_{2}},\xi)]}\leq{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})+\epsilon_{1})=1$ for all $n\geq n_{1}$. Later if we pick $n_{0}\geq n_{1}$, then the almost sure performance guarantee holds for such $\boldsymbol{x}_{n_{0}}^{\epsilon_{2}}$ and ${J}_{n_{0}}^{\epsilon_{1}}(\boldsymbol{x}_{n_{0}}^{\epsilon_{2}})$.

Second, we show a tight certificate bound can be achieved almost surely. Consider performance bound (8). As $\epsilon(\beta_{n})$ decreases and goes to [math] as $n\rightarrow\infty$, there exists $n_{2}$ such that $2\hat{L}\epsilon(\beta_{n})\leq\epsilon_{3}$ holds for all $n\geq n_{2}$. Therefore, we have $\mathbf{P}^{n}\{{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})\leq J^{\star}+\epsilon_{1}+\epsilon_{2}+\epsilon_{3}\}\geq 1-\beta_{n}$ for all $n\geq n_{2}$, or equivalently, $\mathbf{P}^{n}\{{J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})\geq J^{\star}+\epsilon_{1}+\epsilon_{2}+\epsilon_{3}\}\leq\beta_{n}$. As $\sum_{n=1}^{\infty}\beta_{n}<\infty$, then the $1^{\textup{st}}$ Borel-Cantelli Lemma applies to this situation. Thus we claim that there exists a sufficiently large $n_{3}$ such that for all $n\geq\max\{n_{2},n_{3}\}$ we have almost surely, ${J}_{n}^{\epsilon_{1}}(\boldsymbol{x}_{n}^{\epsilon_{2}})\leq J^{\star}+\epsilon_{1}+\epsilon_{2}+\epsilon_{3}$.

Then, by letting $n_{0}:=\max\{n_{1},n_{2},n_{3}\}$ we have almost sure performance guarantee (9) and almost surely, the bound (10).

Lemma VIII.1

Proof:

The proof is an application of the dual characterization of the Wasserstein distance. Let us consider

[TABLE]

By partitioning the data set ${\Xi}_{n}$ into $\mathcal{C}_{n}$ and ${\Xi}_{n}\setminus\mathcal{C}_{n}$ for each summation term, we have

[TABLE]

Canceling the first summation term gives us the following

[TABLE]

where the first inequality is derived taking component-wise absolute values; the second inequality is due to $f$ being in the space of Lipschitz functions defined on $\mathcal{Z}$ with Lipschitz constant 1; and the third inequality is due to $\varsigma\in B_{\omega}(\zeta_{k})$.

Lemma VIII.2

Proof:

From Lemma 5 we have $\mathbf{P}^{n}\{d_{W}(\mathbb{P},\hat{\mathbb{P}}^{n})\leq\epsilon(\beta_{n})\}\geq 1-\beta_{n}$ for each $n$. Then using the result from Lemma VIII.1 we have $\mathbf{P}^{n}\{d_{W}(\mathbb{P},\tilde{\mathbb{P}}^{n})\leq d_{W}(\tilde{\mathbb{P}}^{n},\hat{\mathbb{P}}^{n})+d_{W}(\mathbb{P},\hat{\mathbb{P}}^{n})\leq\epsilon(\beta_{n})+\omega\}\geq 1-\beta_{n}$, i.e., $\mathbf{P}^{n}\{d_{W}(\mathbb{P},\tilde{\mathbb{P}}^{n})\leq\tilde{\epsilon}(\beta_{n})\}\geq 1-\beta_{n}$ for each $n$. The rest of the proof follows directly from that in Lemma 5 and Theorem IV.2.