{'fraction_non_alphanumeric': 0.04818870155093309, 'fraction_numerical': 0.01480564078025164, 'mean_word_length': 4.237800383311219, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 8, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Wasserstein Distributionally Robust Kalman Filtering

We study a distributionally robust mean square error estimation problem over a nonconvex Wasserstein ambiguity set containing only normal distributions. We show that the optimal estimator and the least favorable distribution form a Nash equilibrium. Despite the non-convex nature of the ambiguity set, we prove that the estimation problem is equivalent to a tractable convex program. We further devise a Frank-Wolfe algorithm for this convex program whose direction-searching subproblem can be solved in a quasi-closed form. Using these ingredients, we introduce a distributionally robust Kalman filter that hedges against model risk.


Introduction
The Kalman filter is the workhorse for the online tracking and estimation of a dynamical system's internal state based on indirect observations [1]. It has been applied with remarkable success in areas as diverse as automatic control, brain-computer interaction, macroeconomics, robotics, signal processing, weather forecasting and many more. The classical Kalman filter critically relies on the availability of an accurate state-space model and is therefore susceptible to model risk. This observation has led to several attempts to robustify the Kalman filter against modeling errors.
The H ∞ -filter targets situations in which the statistics of the noise process is uncertain and where one aims to minimize the worst case instead of the variance of the estimation error [3,29]. This filter bounds the H ∞ -norm of the transfer function that maps the disturbances to the estimation errors. However, in transient operation, the desired H ∞ -performance is lost, and the filter may diverge unless some (typically restrictive) positivity condition holds in each iteration. In set-valued estimation the disturbance vectors are modeled through bounded sets such as ellipsoids [6,25]. In this framework, one attempts to construct the smallest ellipsoids around the state estimates that are consistent with the observations and the exogenous disturbance ellipsoids. However, the resulting robust filters ignore any distributional information and thus have a tendency to be over-conservative. A filter that is robust against more general forms of (set-based) model uncertainty was first studied in [22]. This filter iteratively minimizes the worst-case mean square error across all models in the vicinity of a nominal state space model. While performing well in the face of large uncertainties, this filter may be too conservative under small uncertainties. A generalized Kalman filter that addresses this shortcoming and strikes the balance between nominal and worst-case performance has been proposed in [28]. A risk-sensitive Kalman filter is obtained by minimizing the moment-generating function instead of the mean of the squared estimation error [27]. This risk-sensitive Kalman filter is equivalent to a distributionally robust filter proposed in [15], which minimizes the worst-case mean square error across all joint state-output distributions in a Kullback-Leibler (KL) ball around a nominal distribution. Extensions to more general τ -divergence balls are investigated in [30].
In this paper we use ideas from distributionally robust optimization to design a Kalman-type filter that is immunized against model risk. Specifically, we assume that the joint distribution of the states and outputs is uncertain but known to reside in a given ambiguity set that contains all distributions in the proximity of the nominal distribution generated by a nominal state-space model. The ambiguity set thus reflects our level of (dis)trust in the nominal model. We then construct the most accurate filter under the least favorable distribution in this set. The hope is that hedging against the worst-case distribution has a regularizing effect and will lead to a filter that performs well under the unknown true distribution. Distributionally robust filters of this type have been studied in [10,19] using uncertainty sets for the covariance matrix of the state vector and in [15,30] using ambiguity sets defined via information divergences. Inspired by recent progress in data-driven distributionally robust optimization [17], we construct here the ambiguity set as a ball around the nominal distribution with respect to the type-2 Wasserstein distance. The Wasserstein distance has seen widespread application in machine learning [2,9,21], and an intimate relation between regularization and Wasserstein distributional robustness has been discovered in [24,23,26,18]. Also, the Wasserstein distance is known to be more statistically robust than other information divergences [8].
We summarize our main contributions as follows: • We introduce a distributionally robust mean square estimation problem over a nonconvex Wasserstein ambiguity set containing normal distributions only, and we demonstrate that the optimal estimator and the least favorable distribution form a Nash equilibrium. • Leveraging modern reformulation techniques from [18], we prove that this problem is equivalent to a tractable convex program-despite the nonconvex nature of the underlying ambiguity set-and that the optimal estimator is an affine function of the observations. • We devise an efficient Frank-Wolfe-type first-order method inspired by [13] to solve the resulting convex program. We show that the direction-finding subproblem can be solved in quasi-closed form, and we derive the algorithm's convergence rate. • We introduce a Wasserstein distributionally robust Kalman filter that hedges against model risk.
The filter can be computed efficiently by solving a sequence of robust estimation problems via the proposed Frank-Wolfe algorithm. Its performance is validated on standard test instances.
All proofs are relegated to Appendix A, and additional numerical results are reported in Appendix B.
Notation: For any A ∈ R d×d we use Tr [A] to denote the trace and A to denote the spectral norm of A. By slight abuse of notation, the Euclidean norm of v ∈ R d is also denoted by v . Moreover, I d stands for the identity matrix in R d×d . For any A, B ∈ R d×d , we use A, B = Tr A B to denote the trace inner product. The space of all symmetric matrices in R d×d is denoted by S d . We use S d + (S d ++ ) to represent the cone of symmetric positive semidefinite (positive definite) matrices in S d . For any A, B ∈ S d , the relation . Finally, the set of all normal distribution on R d is denoted by N d .

Robust Estimation with Wasserstein Ambiguity Sets
Consider the problem of estimating a signal x ∈ R n from a potentially noisy observation y ∈ R m . In practice, the joint distribution of x and y is never directly observable and thus fundamentally uncertain. This distributional uncertainty should be taken into account in the estimation procedure. In this paper, we model distributional uncertainty through an ambiguity set P, that is, a family of distributions on R d , d = n + m, that are sufficiently likely to govern x and y in view of the available data or that are sufficiently close to a prescribed nominal distribution. We then seek a robust estimator that minimizes the worst-case mean square error across all distributions in the ambiguity set. In the following, we propose to use the Wasserstein distance in order to construct ambiguity sets. Definition 2.1 (Wasserstein distance). The type-2 Wasserstein distance between two distributions Q 1 and Q 2 on R d is defined as where Π(Q 1 , Q 2 ) is the set of all probability distributions on R d × R d with marginals Q 1 and Q 2 .
Proposition 2.2 ([12, Proposition 7]). The type-2 Wasserstein distance between two normal distribu- Consider now a d-dimensional random vector z = [x , y ] comprising the signal x ∈ R n and the observation y ∈ R m , where d = n + m. For a given ambiguity set P, the distributionally robust minimum mean square error estimator of x given y is a solution of the outer minimization problem in where L denotes the family of all measurable functions from R m to R n . Problem (2) can be viewed as a zero-sum game between a statistician choosing the estimator ψ and a fictitious adversary (or nature) choosing the distribution Q. By construction, the minimax estimator performs best under the worst possible distribution Q ∈ P. From now on we assume that P is the Wasserstein ambiguity set which can be interpreted as a ball of radius ρ ≥ 0 in the space of normal distributions. We will further assume that P is centered at a normal distribution P = N d (µ, Σ) with covariance matrix Σ 0.
Even though the Wasserstein ambiguity set P is nonconvex (as mixtures of normal distributions are generically not normal), we can prove a minimax theorem, which ensures that one may interchange the infimum and the supremum in (2) without affecting the problem's optimal value.
Remark 2.4 (Connection to Bayesian estimation). The optimal solutions ψ and Q of the two dual problems in (4) represent the minimax strategies of the statistician and nature, respectively. Theorem 2.3 implies that (ψ , Q ) forms a saddle point (and thus a Nash equilibrium) of the underlying zero-sum game. Hence, the robust estimator ψ is also the optimal Bayesian estimator for the prior Q . For this reason, Q is often referred to as the least favorable prior [14].
We now demonstrate that the minimax problem (2) is equivalent to a tractable convex program, whose solution allows us to recover both the optimal estimator ψ as well as the least favorable prior Q .
Theorem 2.5 (Tractable reformulation). The minimax problem (2) with the Wasserstein ambiguity set (3) centered at P = N d (µ, Σ),σ λ min (Σ) > 0, is equivalent to the finite convex program If S , S xx , S yy and S xy is optimal in (5) and µ = [µ x , µ y ] for some µ x ∈ R n and µ y ∈ R m , then the affine function ψ (y) = S xy (S yy ) −1 (y − µ y ) + µ x is the distributionally robust minimum mean square error estimator, and the normal distribution Q = N d (µ, S ) is the least favorable prior.
Theorem 2.5 provides a tractable procedure for constructing a Nash equilibrium (ψ , S ) for the statistician's game against nature. Note that if ρ = 0, then S = Σ is the unique solution to (5). In this case the estimator ψ reduces to the Bayesian estimator corresponding to the nominal distribution P = N d (µ, Σ). We emphasize that the choice of the Wasserstein radius ρ may have a significant impact on the resulting estimator. In fact, this is a key distinguishing feature of the Wasserstein ambiguity set (3) with respect to other popular divergence-based ambiguity sets.
Remark 2.6 (Divergence-based ambiguity sets). As a natural alternative, one could replace the Wasserstein distance in (3) with an information divergence. For example, ambiguity sets defined via τ -divergences, which encapsulate the popular KL divergence as a special case, have been studied in [15,30]. As shown in [15,Theorem 1] and [30, Theorem 2.1], the optimal estimator corresponding to any τ -divergence ambiguity set always coincides with the Bayesian estimator for the nominal distribution P = N d (µ, Σ) irrespective of ρ. Thus, in stark contrast to the setting considered here, the size of a τ -divergence ambiguity set has no impact on the corresponding optimal estimator. Moreover, the least favorable prior Q = N d (µ, S ) for a τ -divergence ambiguity set always satisfies Thus, in order to harm the statistician, nature only perturbs the second moments of the signal but sets all second moments of the observation as well as all cross moments to their nominal values. Example 2.7 (Impact of ρ on the Nash equilibrium). We illustrate the dependence of the saddle point (ψ , Q ) on the size ρ of the ambiguity set in a 2-dimensional example. Suppose that the nominal distribution P of [x, y] ∈ R 2 satisfies µ x = µ y = 0, Σ xx = Σ xy = 1 and Σ yy = 1.1, implying that the noise w y − x and the signal x are independent (E P [xw] = Σ xy − Σ xx = 0). Figure 1 visualizes the canonical 90% confidence ellipsoids of the the least favorable priors as well as the graphs of the optimal estimators for different sizes of the Wasserstein and KL ambiguity sets. As ρ increases, the least favorable prior for the Wasserstein ambiguity set displays the following interesting properties: (i) the signal variance S xx increases, (ii) the measurement variance S yy decreases, (iii) the signal-measurement covariance S xy decreases towards 0, and (iv) the noise variance decreases and is negative for all ρ > 0, and (vi) the optimal estimator ψ tends to the zero function. Note that the optimal estimator and the measurement variance remain constant in ρ when working with a KL ambiguity set.

Efficient Frank-Wolfe Algorithm
The finite convex optimization problem (5) is numerically challenging as it constitutes a nonlinear semi-definite program (SDP). In principle, it would be possible to eliminate all nonlinearities by using Schur complements and to reformulate (5) as a linear SDP, which is formally tractable. However, it is folklore knowledge that general-purpose SDP solvers are yet to be developed that can reliably solve large-scale problem instances. We thus propose a tailored first-order method to solve the nonlinear SDP (5) directly, which exploits a covert structural property of the problem's objective function Definition 3.1 (Unit total elasticity 1 ). We say that a function ϕ : S d + → R + has unit total elasticity if It is clear that every linear function has unit total elasticity. Maybe surprisingly, however, the objective function f (S) of problem (5) also enjoys unit total elasticity because Moreover, as will be explained below, it turns out problem (5) can be solved highly efficiently if its objective function is replaced with a linear approximation. These observations motivate us to solve (5) with a Frank-Wolfe algorithm [11], which starts at S (0) = Σ and constructs iterates where α k represents a judiciously chosen step-size, while the oracle mapping F : S + → S + returns the unique solution of the direction-finding subproblem In each iteration, the Frank-Wolfe algorithm thus minimizes a linearized objective function over the original feasible set. In contrast to other commonly used first-order methods, the Frank-Wolfe algorithm thus obviates the need for a potentially expensive projection step to recover feasibility. It is easy to convince oneself that any solution of the nonlinear SDP (5) is indeed a fixed point of the operator F . To make the Frank-Wolfe algorithm (7) work in practice, however, one needs (i) an efficient routine for solving the direction-finding subproblem (7b); (ii) a step-size rule that offers rigorous guarantees on the algorithm's convergence rate.
In the following, we propose an efficient bisection algorithm to address (i). As for (ii), we show that the convergence analysis portrayed in [13] applies to the problem at hand. The procedure for solving (7b) is outlined in Algorithm 1, which involves an auxiliary function h : Theorem 3.2 (Direction-finding subproblem). For any fixed inputs ρ, ε ∈ R ++ , Σ ∈ S d ++ and S ∈ S d + , Algorithm 1 outputs a feasible and ε-suboptimal solution to (7b).
We emphasize that the most expensive operation in Algorithm 1 is the matrix inversion (γI d − D) −1 , which needs to be evaluated repeatedly for different values of γ. These computations can be accelerated by diagonalizing D only once at the beginning. The repeat loop in Algorithm 1 carries out the actual bisection algorithm, and a suitable initial bisection interval is determined by a pair of a priori bounds LB and U B, which are available in closed form (see Appendix A).
The overall structure of the proposed Frank-Wolfe method is summarized in Algorithm 2. We borrow the step-size rule suggested in [13] to establish rigorous convergence guarantees. This is accomplished by showing that the objective function f has a bounded curvature constant. Our convergence result is formalized in the next theorem. Theorem 3.3 (Convergence analysis). If Σ 0, ρ > 0, δ > 0 and α k = 2/(2 + k) for any k ∈ N, then the k-th iterate S (k) computed by Algorithm 2 is feasible in (5) and satisfies where S is an optimal solution of (5),σ is the smallest eigenvalue of Σ, andσ (ρ + Tr [Σ]) 2 .
1 Our terminology is inspired by the definition of the elasticity of a univariate function ϕ(s) as dϕ(s) ds s ϕ(s) .

The Wasserstein Distributionally Robust Kalman Filter
Consider a discrete-time dynamical system whose (unobservable) state x t ∈ R n and (observable) output y t ∈ R m evolve randomly over time. At any time t ∈ N, we aim to estimate the current state x t based on the output history Y t (y 1 , . . . , y t ). We assume that the joint state-output process z t = [x t , y t ] , t ∈ N, is governed by an unknown Gaussian distribution Q in the neighborhood of a known nominal distribution P . The distribution P is determined through the linear state-space model where A t , B t , C t , and D t are given matrices of appropriate dimensions, while v t ∈ R d , t ∈ N, denotes a Gaussian white noise process independent of x 0 ∼ N n (x 0 , V 0 ). Thus, v t ∼ N d (0, I d ) for all t, while v t and v t are independent for all t = t . Note that we may restrict the dimension of v t to the dimension d = n + m of z t without loss of generality. Otherwise, all linearly dependent columns of [B t , D t ] and the corresponding components of v t can be eliminated systematically.
By the law of total probability and the Markovian nature of the state-space model (9), the nominal distribution P is uniquely determined by the marginal distribution P x0 = N n (x 0 , V 0 ) of the initial state x 0 and the conditional distributions Unlike P , the true distribution Q governing z t , t ∈ N, is unknown, and thus the estimation problem at hand is not well-defined. We will therefore estimate the conditional meanx t and covariance matrix V t of x t given Y t under some worst-case distribution Q to be constructed recursively. First, we assume that the marginal distribution Q x0 of x 0 under Q equals P x0 , that is, Q x0 = N n (x 0 , V 0 ). Next, fix any t ∈ N and assume that the conditional distribution Q xt−1|Yt−1 of x t−1 given Y t−1 under Q has already been computed as Q xt−1|Yt−1 = N n (x t−1 , V t−1 ). The construction of Q xt|Yt is then split into a prediction step and an update step. The prediction step combines the previous state estimate Q xt−1|Yt−1 with the nominal transition kernel P zt|xt−1 to generate a pseudo-nominal Algorithm 3 Robust Kalman filter at time t Input: Covariance matrix V t−1 0 State estimatex t−1 Wasserstein radius ρ t > 0 Tolerance δ > 0 Prediction: Form the pseudo-nominal distribution P zt|Yt−1 = N d (µ t , Σ t ) using (10) Observation: Observe the output y t Update: Use Algorithm 2 to solve (11) x t = S t,xy (S t,yy ) −1 (y t − µ t,y ) + µ t,x Figure 2: Wasserstein ball in the space S 2 + of covariance matrices centered at I 2 with radius 1.
for every Borel set B ⊆ R d and observation history Y t−1 ∈ R m×(t−1) . The well-known formula for the convolution of two multivariate Gaussians reveals that Note that the construction of P zt|Yt−1 resembles the prediction step of the classical Kalman filter but uses the least favorable distribution Q xt−1|Yt−1 instead of the nominal distribution P xt−1|Yt−1 .
In the update step, the pseudo-nominal a priori estimate P zt|Yt−1 is updated by the measurement y t and robustified against model uncertainty to yield a refined a posteriori estimate Q xt|Yt . This a posteriori estimate is found by solving the minimax problem equipped with the Wasserstein ambiguity set Note that the Wasserstein radius ρ t quantifies our distrust in the pseudo-nominal a priori estimate and can therefore be interpreted as a measure of model uncertainty. Practically, we reformulate (11) as an equivalent finite convex program of the form (5), which is amenable to efficient computational solution via the Frank-Wolfe algorithm detailed in Section 3. By Theorem 2.5, the optimal solution S t of problem (5) yields the least favorable conditional distribution Q zt|Yt−1 = N d (µ t , S t ) of z t given Y t−1 . By using the well-known formulas for conditional normal distributions (see, e.g., [20, page 522]), we then obtain the least favorable conditional distribution The distributionally robust Kalman filtering approach is summarized in Algorithm 3. Note that the robust update step outlined above reduces to the usual update step of the classical Kalman filter for ρ ↓ 0.

Numerical Results
We showcase the performance of the proposed Frank-Wolfe algorithm and the distributionally robust Kalman filter in a suite of synthetic experiments. All optimization problems are implemented in MATLAB and run on an Intel XEON CPU with 3.40GHz clock speed and 16GB of RAM, and the corresponding codes are made publicly available at https://github.com/sorooshafiee/WKF.

Distributionally Robust Minimum Mean Square Error Estimation
We first assess the distributionally robust minimum mean square error (robust MMSE) estimator, which is obtained by solving (2), against the classical Bayesian MMSE estimator, which can be viewed as the solution of problem (2) over a singleton ambiguity set that contains only the nominal distribution. Recall from Remark 2.6 that the optimal estimator corresponding to a KL or τ -divergence ambiguity set of the type studied in [15,30] coincides with the Bayesian MMSE estimator irrespective of ρ. Thus, we may restrict attention to Wasserstein ambiguity sets. In order to develop a geometric intuition, Figure 2 visualizes the set of all bivariate normal distributions with zero mean that have a Wasserstein distance of at most 1 from the standard normal distribution-projected to the space of covariance matrices.
In the first experiment we aim to predict a signal x ∈ R 4d/5 from an observation y ∈ R d/5 , where the random vector z = [x , y ] follows a d-variate Gaussian distribution with d ∈ {10, 50, 100}. The experiment comprises 10 4 simulation runs. In each run we randomly generate two covariance matrices Σ and Σ as follows. First, we draw two matrices A and A from the standard normal distribution on R d×d , and we denote by R and R the orthogonal matrices whose columns correspond to the orthonormal eigenvectors of A + (A ) and A + A , respectively. Then, we define ∆ = R Λ (R ) and Σ = RΛR , where Λ and Λ are diagonal matrices whose main diagonals are sampled uniformly from [0, 1] d and [0.1, 10] d , respectively. Finally, we set Σ = (Σ 1 2 + (∆ ) 1 2 ) 2 and define the normal distributions P = N d (0, Σ ) and P = N d (0, Σ). By construction, we have where · F stands for the Frobenius norm, and the first inequality follows from [16,Proposition 3]. We assume that P is the true distribution and P our nominal prior. The robust MMSE estimator is obtained by solving (5) for ρ = √ d via the Frank-Wolfe algorithm from Section 3, while the Bayesian MMSE estimator under P is calculated analytically. In order to provide a meaningful comparison between these two approaches, we also compute the Bayesian MMSE estimator under the true distribution P (denoted by MMSE ), which is indeed the best possible estimator. Figure 3 visualizes the distribution of the difference between the mean square errors under P of the robust MMSE (Bayesian MMSE) and MMSE estimators. We observe that the robust MMSE estimator produces better results consistently across all experiments, and the effect is more pronounced for larger dimensions d. Figures 4(a) and 4(b) report the execution time and the iteration complexity of the Frank-Wolfe algorithm for d ∈ {10, . . . , 100} when the algorithm is stopped as soon as the relative duality gap F (S k ) − S k , ∇f (S k ) /f (S k ) drops below 0.01%. Note that the execution time grows polynomially due to the matrix inversion in the bisection algorithm. Figure 4(c) shows the relative duality gap of the current solution as a function of the iteration count.

Wasserstein Distributionally Robust Kalman Filtering
We assess the performance of the proposed Wasserstein distributionally robust Kalman filter against that of the classical Kalman filter and the Kalman filter with the KL ambiguity set from [15]. To this end, we borrow the standard test instance from [22,28,15] with n = 2 and m = 1. In all numerical experiments we simulate the different filters over 1000 periods starting fromx 0 = 0 and V 0 = I 2 . Figure 5 shows the empirical mean square error 1 500 500 j=1 x j t −x j t 2 across 500 independent simulation runs, wherex j t denotes the state estimate at time t in the j th run. We distinguish four different scenarios: time-invariant uncertainty (∆ j t = ∆ j sampled uniformly from [−∆,∆] for each j) versus time-varying uncertainty (∆ j t sampled uniformly from [−∆,∆] for each t and j), and small uncertainty (∆ = 1) versus large uncertainty (∆ = 10). All results are reported in decibel units (10 log 10 (·)). As for the filter design, the Wasserstein and KL radii are selected from the search grids {a · 10 −1 : a ∈ {1, 1.1, · · · , 2}} and {a · 10 −4 : a ∈ {1, 1.1, · · · , 2}}, respectively. Figure 5 reports the results with minimum steady state error across all candidate radii.
Under small time-invariant uncertainty ( Figure 5(a)), the Wasserstein and KL distributionally robust filters display a similar steady-state performance but outperform the classical Kalman filter. Note that the KL distributionally robust filter starts from a different initial point as we use the delayed implementation from [15]. Under small time-varying uncertainty ( Figure 5(b)), both distributionally robust filters display a similar performance as the classical Kalman filter. Figures 5(c) and (d) corresponding to the case of large uncertainty are similar to Figures 5(a) and (b), respectively. However, the Wasserstein distributionally robust filter now significantly outperforms the classical Kalman filter and, to a lesser extent, the KL distributionally robust filter. Moreover, the Wasserstein distributionally robust filter exhibits the best transient behavior.

Appendix A Proofs
A.1 Proof of Theorem 2.3 The proof of Theorem 2.3 requires the following preparatory lemma, which we borrow from [18].
Moreover, if γI d D, the unique optimal solution of the above maximization problem is given by Proof of Theorem 2.3. The optimal value of the minimax problem (2) satisfies where (A.1a) follows from the max-min inequality, and (A.1b) holds because the inner minimization problem over ψ is solved by the conditional expectation function ψ (y) = E Q [x|y], which is affine in y for every fixed Gaussian distribution Q ∈ P, see, e.g., [20, page 522]. Without loss of generality, one can thus restrict the set of measurable functions L to the set of affine functions parametrized by a sensitivity matrix G ∈ R n×m and an intercept vector g ∈ R n . Recalling the definition of the Wasserstein ambiguity set P in (3) and encoding each normal distribution Q ∈ P by its mean vector c ∈ R d and covariance matrix S ∈ S d + , we can use Proposition 2.2 to reformulate (A.1b) as sup inf Solving the inner minimization problem over g analytically and substituting the optimal solution g = c x − Gc y back into the objective function shows that (A.2a) is equivalent to (A.2b) The minimization over G may now be interchanged with the maximization over c and S by using the classical minimax theorem [5,Proposition 5.5.4], which applies because c and S range over a compact feasible set. The inner maximization problem over c is then solved by c = µ, which maximizes the slack of the Wasserstein constraint. Thus, the minimax problem (A.2b) simplifies to Assigning a Lagrange multiplier γ ≥ 0 to the Wasserstein constraint and dualizing the inner maximization problem yields Strong duality holds because S = Σ 0 represents a Slater point for the primal maximization problem. Finally, by using Lemma A.1, problem (A.2d) can be reformulated as By construction, the optimal value of (A.3) provides a lower bound on that of the minimax problem (2). Next, we construct an upper bound by restricting L to the class of affine estimators.
As P is non-convex, we cannot simply use Sion's minimax theorem to show that the right-hand side of (A.4) equals (A.1b). Instead, we need a more involved argument. Recalling the definition of P in (3) and encoding each normal distribution Q ∈ P by its mean vector c ∈ R d and covariance matrix S ∈ S d + , we can use Proposition 2.2 to reformulate the right-hand side of (A.4) as inf (A.5a) Next, we introduce the set C {c ∈ R d : c − µ ≤ ρ} as well as the auxiliary functions We emphasize that the constraint c ∈ C is redundant in (A.5b) but will facilitate further simplifications below. Note also that D(G) 0, and thus the minimax problem (A.5b) involves a cumbersome convex maximization problem over c. By employing a penalty formulation of the Wasserstein constraint, the inner maximization problem over c and S in (A.5b) can be re-expressed as Here, the minimization over γ and the maximization over S may be interchanged by strong duality, which holds because S = Σ 0 constitutes a Slater point for the primal problem, see, e.g., [5,Proposition 5.3.1]. We note that when c − µ = ρ, the feasible set of S reduces to a singleton, and thus strong duality holds trivially. The emerging inner maximization problem over S can then be solved analytically by using Lemma A.1. In summary, the minimax problem (A.5b) is equivalent to inf Observe now that the optimal value function of the innermost minimization problem over γ in (A.5c) is convex in g and, thanks to the constraint γI d − D(G) 0, concave in c for every fixed G. By the classical minimax theorem [5,Proposition 5.5.4], which applies because c ranges over the compact set C, we may thus interchange the infimum over g with the supremum over c. After replacing D(G) and b(G, g) with their definitions, it becomes clear that the innermost minimization problem over g admits the analytical solution g = µ x − Gµ y . Thus, problem (A.5c) is equivalent to (A.5d) By invoking the minimax theorem [5,Proposition 5.5.4] once again, the inner infimum over γ can be interchanged with the supremum over c. As the resulting inner maximization problem over c is solved by c = µ, problem (A.5d) is thus equivalent to (A.3). In summary, we have shown that (A.3) provides both an upper bound on the left-hand side of (A.1) as well as a lower bound on the right-hand side of (A.1). Thus, the inequality in (A.1) is in fact an equality.

A.2 Proof of Theorem 2.5
The proof of Theorem 2.5 relies on the following lemma, which extends a similar result from [18].
Lemma A.2 (Analytical solution of direction-finding subproblem). For any fixed Σ ∈ S d ++ and D ∈ S d + \{0}, the optimization problem where γ is the unique solution with γ I d D of the algebraic equation Moreover, we have S σ I d , whereσ λ min (Σ).
Proof of Lemma A.2. The optimality of S follows immediately from [18,Theorem 5.1]. Moreover, the spectral norm of (S ) −1 obeys the following estimate.
As the largest eigenvalue of (S ) −1 is bounded byσ −1 , we may conclude that S σ I d .
Proof of Theorem 2.5. The proof of Theorem 2.3 has shown that the original infinite-dimensional minimax problem (2) is equivalent to the finite-dimensional minimax problem (A.2c). By Lemma A.2, the solution of the inner maximization problem in (A.2c) satisfies S σ I d . Thus, one may append the redundant constraint S σ I d to this inner problem without sacrificing optimality. By interchanging the minimization over G with the maximization over S, which is allowed by [5,Proposition 5.5.4], problem (A.2c) can thus be reformulated as (A.6) Recall thatσ > 0, which implies that S 0. Hence, the unconstrained quadratic minimization problem over G in (A.6) has a unique solution G , which can be obtained analytically by solving the problem's first-order optimality condition. Specifically, we have Substituting G into (A.6) yields the desired maximization problem (5). By construction, this convex program is equivalent to nature's decision problem on the right-hand side of (4), and thus it is easy to see that the least favorable prior is given by Q = N d (µ, S ). Next, we solve the Bayesian estimation problem inf An elementary analytical calculation reveals that this problem is solved by ψ (y) = S xy (S yy ) −1 (y − µ y ) + µ x . Moreover, this solution is unique because S σ I d , which implies that the objective function is strictly convex. By Theorem 2.3 and [7, Section 5.5.5], we may then conclude that ψ is also optimal in (2). This observation completes the proof.

A.3 Proof of Theorem 3.2
The following lemma suggests upper and lower bounds on the (unique) root γ of the function h(γ) defined in (8). Note that this root is computed approximately using bisection in Algorithm 1. the scalar λ 1 is the largest eigenvalue of D ∇f (S), and v 1 is a corresponding eigenvector.
The function h can be equivalently rewritten as where the summation admits the following bounds: Equating the two bounds to ρ 2 yields γ min and γ max , respectively.
Proof of Theorem 3.2. The proof of Lemma A.2 implies that L(γ) γ 2 (γI d − D) −1 Σ(γI d − D) −1 is feasible in (7b) for every γ with γI d D and h(γ) > 0. Moreover, L(γ ) is optimal in (7b) if γ I d D and h(γ ) = 0. Algorithm 1 uses a bisection procedure to compute an approximation γ of γ such that L(γ) is feasible and ε-suboptimal in (7b). The degree of suboptimality of L(γ) equals L(γ ) − L(γ), D . The true optimal value L(γ ), D is inaccessible but can be estimated above by the objective value of γ in the Lagrangian dual of (7b), which can be expressed as The proof of Theorem 3.3 widely parallels that of [13,Theorem 1]. The key ingredient is to prove that the curvature constant of the problem's (negative) objective function −f is bounded. Definition A.4 (Curvature constant). The curvature constant C g of the convex function g with respect to a compact domain S is defined as In order to bound the curvature constant of −f , we need several preparatory lemmas. Proof of Lemma A.6. We seek an upper bound on the maximum eigenvalue of S uniformly across all covariance matrices S feasible in (5), that is, we seek an upper bound on the optimal value of sup S 0 S s. t. Tr S + Σ − 2 Σ Proof of Lemma A.7. We first expand the negative objective function −f at S ∈ S d + . By Lemma A.5, for any symmetric perturbation matrix ∆ with a characteristic block structure of the form ∆ = ∆ xx ∆ xy ∆ xy ∆ yy ∈ S d , the negative objective function −f (S + ∆) can be expressed as Tr −S xx − ∆ xx + (S xy + ∆ xy )(S yy + ∆ yy ) −1 (S yx + ∆ xy ) = Tr [−S xx − ∆ xx ] + Tr (S xy + ∆ xy )S −1 yy (I m − ∆ yy S −1 yy + (∆ yy S −1 yy ) 2 + O( ∆ yy 3 ))(S yx + ∆ xy ) = Tr −S xx + S xy S −1 yy S yx − Tr ∆ xx − ∆ xy S −1 yy S yx + S xy S −1 yy ∆ yy S −1 yy S yx − S xy S −1 yy ∆ xy − Tr ∆ xy S −1 yy ∆ yy S −1 yy S yx − ∆ xy S −1 yy ∆ xy + S xy S −1 yy ∆ yy S −1 yy ∆ xy − S xy S −1 yy (∆ yy S −1 yy ) 2 S yx + O( The diameter of the feasible set S with respect to the Frobenius norm satisfies diam(S) = sup S1,S2∈S where the first inequality holds due to [4,Equation (9.2.16)], and the last inequality follows from the proof of Lemma A.6. Therefore, by [13,Lemma 7], the curvature constant C −f admits the estimate C −f ≤ (diam(S)) 2 Lip(∇f ) ≤ 2σ 4 σ 3 . This observation completes the proof.
Proof of Theorem 3.3. By Lemma A.7, the curvature constant C −f is bounded, and thus one can directly apply [13, Theorem 1] to find the convergence rate.