A Robust Unscented Transformation for Uncertain Moments

TL;DR

This paper introduces a robust version of the unscented transform that accounts for uncertain moments within intervals, using polynomial optimization and relaxation techniques to improve robustness in one-dimensional stochastic analysis.

Contribution

It develops a robust UT framework based on Chebychev centers for uncertain moments, reformulating the problem as a polynomial optimization with proposed relaxation algorithms.

Findings

Provides a robust UT method for uncertain moments.

Formulates the problem as a polynomial optimization.

Proposes algorithms to relax NP-hard problems.

Abstract

This paper proposes a robust version of the unscented transform (UT) for one-dimensional random variables. It is assumed that the moments are not exactly known, but are known to lie in intervals. In this scenario, the moment matching equations are reformulated as a system of polynomial equations and inequalities, and it is proposed to use the Chebychev center of the solution set as a robust UT. This method yields a parametrized polynomial optimization problem, which in spite of being NP-Hard, can be relaxed by some algorithms that are proposed in this paper.

Equations

$$\text{E}\left\{f(X)\right\}={\int }_{\mathbb{R}}f(\xi )p_X(\xi )d\xi .$$

$$E\left\{f(Z)\right\}=\sum _{i=1}^{m+1}f(z_i)p_Z(z_i).$$

$$E\{X^k\}$$

$$=\sum _{i=1}^{m+1}{z}_i^k{p}_Z(z_i),k=1,\dots ,m,$$

$$E\left\{g_k(Z)\right\}=E\left\{g_k(X)\right\},k=1,\dots ,m,$$

$$E\{g_k(X)\}=\sum _{i=1}^{m+1}{\theta }_i{g}_k(z_i),k=1,\dots ,m$$

$$\sum _{i=1}^{m+1}{\theta }_i=1.$$

$$p_Z(k)=\{\begin{array}{ll}{\theta }_i,& \text{if }k=z_i,\\ 0,& \text{otherwise,}\end{array}$$

$$\frac{f^{(2m)}(\xi )}{(2m)!}{\int }_a^b{p}_X(x)h_m^2(x)dx,$$

$$\begin{array}[]{ccc}z_{1}&=&\mu-\sqrt{3V},\\ z_{2}&=&\mu,\\ z_{3}&=&\mu+\sqrt{3V},\end{array}$$

$$z_1=$$

$$z_3=$$

$$z_1=$$

$$z_3=$$

$$S:\qquad\begin{array}[]{cccc}\omega_{j}&\geq&0,&j\in\{1,\ldots,m\},\\ \sum_{i=1}^{m+1}\omega_{i}&=&1,\\ \sum_{i=1}^{m+1}z_{i}^{k_{1}}\omega_{i}&=&E\{X^{k_{1}}\},&k_{1}\in\mathcal{K}_{1},\\ \sum_{i=1}^{m+1}z_{i}^{k_{2}}\omega_{i}&\leq&u_{k_{2}},&k_{2}\in\mathcal{K}_{2},\\ \sum_{i=1}^{m+1}z_{i}^{k_{2}}\omega_{i}&\geq&\ell_{k_{2}},&k_{2}\in\mathcal{K}_{2},\end{array}$$

$$(z_1^C,\dots ,z_{m+1}^C,{\omega }_1^C,\dots ,{\omega }_{m+1}^C)=\underset{\hat{x}\in {\mathbb{R}}^{2m+2}}{\mathrm{argmin}}\underset{x\in S}{\max }\Vert x-\hat{x}{\Vert }^2,$$

$$x=(z_1,\dots ,z_{m+1},{\omega }_1,\dots ,{\omega }_{m+1})$$

$$S_{\varepsilon}:\qquad\begin{array}[]{cccc}\omega_{j}&\geq&\varepsilon,&j\in\{1,\ldots,m\},\\ \sum_{i=1}^{m+1}\omega_{i}&=&1,\\ \sum_{i=1}^{m+1}z_{i}^{k_{1}}\omega_{i}&=&E\{X^{k_{1}}\},&k_{1}\in\mathcal{K}_{1},\\ \sum_{i=1}^{m+1}z_{i}^{k_{2}}\omega_{i}&\leq&u_{k_{2}},&k_{2}\in\mathcal{K}_{2},\\ \sum_{i=1}^{m+1}z_{i}^{k_{2}}\omega_{i}&\geq&\ell_{k_{2}},&k_{2}\in\mathcal{K}_{2}.\end{array}$$

$$(z_1^C,\dots ,z_{m+1}^C,{\omega }_1^C,\dots ,{\omega }_{m+1}^C)=\underset{\hat{x}\in {\mathbb{R}}^{2m+2}}{\mathrm{argmin}}\underset{x\in S_{\epsilon }}{\max }\Vert x-\hat{x}{\Vert }^2,$$

$$\sum _{i=1}^{m+1}{z}_i^2{\omega }_i\le u_2,$$

$${\underline{z}}_i$$

$${\underline{\omega }}_i$$

$${\overline{z}}_i$$

$${\overline{\omega }}_i$$

$$B=\{(z_1,\dots ,z_{m+1},{\omega }_1,\dots ,{\omega }_{m+1})\in {\mathbb{R}}^{2m+2}:{\underline{z}}_i\le z_i\le {\overline{z}}_i,$$

$${\underline{\omega }}_i\le {\omega }_i\le {\overline{\omega }}_i,i=1,\dots ,m+1\}$$

$$d\le \frac{\text{diam}(B)}{2},$$

$$(z_1^C,\dots ,z_{m+1}^C,{\omega }_1^C,\dots ,{\omega }_{m+1}^C)=\underset{\hat{x}\in B}{\mathrm{argmin}}\underset{x\in S_{\epsilon }}{\max }\Vert x-\hat{x}{\Vert }^2.$$

$$J(\hat{x}):=\underset{x\in S_{\epsilon }}{\max }\Vert x-\hat{x}{\Vert }^2,$$

$$(z_1^C,\dots ,z_{m+1}^C,{\omega }_1^C,\dots ,{\omega }_{m+1}^C)=\underset{\hat{x}\in B}{\mathrm{argmin}}{\tilde{J}}_{\tau }(\hat{x}),$$

$$\begin{array}[]{cc}\min_{\hat{x},r}r&\text{ s.t. }\|x-\hat{x}\|^{2}\leq r,\\ &x\in S_{\varepsilon},\,\hat{x}\in B,\end{array}$$

$$U{T}_f(z_1,\dots ,z_{m+1},w_1,\dots ,w_{m+1}):=\sum _{i=1}^{m+1}{w}_{i}f(z_i).$$

$$\min D\text{s.t. }$$

$$\le (1+D)\Vert x-y\Vert ,$$

Full text

t1t1footnotetext: ©2019. This manuscript version is made available under the CC BY-NC-ND 4.0 license http://creativecommons.org/licenses/by-nc-nd/4.0/. This article has been accepted for publication in a future issue of the Journal of Franklin Institute, but has not been fully edited. Content may change prior to final publication. The final version of record is available at https://doi.org/10.1016/j.jfranklin.2019.02.018.

1 Introduction

In numerous problems of statistics and stochastic filtering, one is often interested in calculating the posterior expectation of a continuous random variable $X$ that undergoes a nonlinear transform $f$, viz.:

[TABLE]

It is not always possible to have a closed-form expression for this integral in terms of elementary functions: if this integral does not satisfy the hypothesis of Liouville’s theorem (see for instance, Section 12.4 of [1]), then the antiderivative of this integral cannot be expressed in terms of elementary functions. Thus, instead of using analytical methods to calculate (1), in many situations numerical methods must be employed.

A common way to numerically calculate (1) is by using the Monte Carlo integration method [2], which is a stochastic sampling technique: by taking a sufficiently large number of samples of the random variable $X$, one can approximate the probability density function $p_{X}$ and obtain an estimate for (1). However, this method can be very demanding computationally, since it frequently employs several thousands of simulations to obtain the statistics of the final result.

Another way to calculate (1) with less computational burden than Monte Carlo integration method is the technique of Unscented Transform (UT). Originally proposed for the problem of extending the Kalman filter for nonlinear dynamical systems [3], this method has also been applied in several problems of engineering, such as in the analysis of the sensitivity of antennas [4] and in circuit optimization [5].

Different from the Monte Carlo integration method, the UT is a deterministic sampling technique: in place of choosing random points to approximate $p_{X}$, points known as sigma points are deterministically selected to capture the statistics of $p_{X}$. This is accomplished by generating a discrete distribution $p^{\prime}_{X}$ having the same first and second (and possibly higher) moments of $p_{X}$. The mean, covariance and higher moments of the transformed ensemble of sigma points can then be computed as the estimate of the nonlinear transformation $f$ of the original distribution [3], [6]. Thus, at least the first moments of $p_{X}$ must exist and be exactly known. In several circumstances, however, this is not valid: for instance, it may be the case that the distribution does not even have the first moment (one example is the Cauchy distribution [7]). In the case that the moments can be assumed to exist, it is usual that they are not precisely known, but it is still possible to have upper and lower bounds for the exact value of moments from statistical experiments or by the use, for instance, of Chebychev’s inequality [2].

In this paper it is designed a technique for generating sigma points when the exact value of the moments are not known, but upper and lower bounds for the unknown moments are known. In this case, the moment matching equations of UT are no longer just a system of polynomial equations but a system of polynomial equations with polynomial inequalities. Furthermore, since this system can have more than one solution, it is possible to choose a set of sigma points which minimizes a given cost function by formulating the problem as a polynomial optimization problem (POP). Although the solution to this problem is in general computationally infeasible [8], by using Lasserre’s hierarchy of semidefinite programming relaxations one can approximate the solution of the original POP problem by the solution of a computationally feasible convex optimization problem[9, 10].

The main contribution of this paper is the introduction of the concept of UT robustness in the sense of exploiting the upper and lower bounds for moments. A robust UT is proposed in [11] but robustness has different meaning and it is achieved by matching precisely known high order moments. It is important to note that Lasserre’s hierarchy of relaxations was applied in earlier investigations of the moment matching problem [12], but its use was limited to polynomial equations while here polynomial inequalities are also taken into account.

This paper is organized as follow. Some preliminaries for the robust UT are presented in Section 2. The robust UT itself is detailed in Section 3. Details about the computation of robust sigma points is presented in Section 4. The computation of UT transform is presented in Section 5. Finally, the conclusion is presented in Section 7.

2 Unscented Transform

The rationale behind the unscented transformation is that it is easier to calculate a moment for a discrete distribution than a continuous one. In fact, to compute the posterior distribution of a random scalar variable $X$ with distribution $p_{X}$ by a function $f$, one must calculate the integral (1) while, in other hand, if $Z$ is a discrete random variable with $m+1$ atoms $z_{i}$ and distribution $p_{Z}$, one only needs to calculate the sum

[TABLE]

Henceforth, if it is possible to choose the atoms and their weights of an adequate $p_{Z}$ to approximate $p_{X}$, then the value $E\{f(Z)\}$ would be a good approximation to $E\{f(X)\}$, but with (2) being easier to compute than (1). Thus, the principle behind the UT is to approximate the continuous distribution $p_{X}$ by the discrete distribution $p_{Z}$ by equating the first $m$ moments of these distributions. In other words, the following equations must be satisfied:

[TABLE]

where is supposed the access to $E\{X^{k}\}$ for $k=1,\ldots,m$.

Given $E\{X^{k}\}$ it is always possible to find $z_{i}$ and $p_{Z}(z_{i})$, $i=1,\ldots,m+1$ (which will be called sigma point and its weight, respectively) satisfying (4). In fact, as it will be stated next, there is at least one solution for the equations

[TABLE]

where $g_{k}:\mathbb{R}\rightarrow\mathbb{R}$ are continuous functions, $g_{k}\not\equiv g_{j}$ for $k\neq j$ and from which (4) is a particular case with $g_{k}(x)=x^{k}$. That (5) has at a least one solution is stated next in Theorem 1.

Theorem 1**.**

Consider that the $m$ moments $E\{g_{k}(X)\}$, $k=1,\ldots,m$, are given. The system (5) in terms of variables $z_{i}$ and $p_{Z}(z_{i})$ has at least one solution with at most $m+1$ sigma points.111The proof of the theorem is known in the mathematical literature in the context of the Caratheodory’s theorem. Since it is less known in the context of UT literature, the proof is presented here for easy reference. See, e.g., [13].

Proof.

Since the point $P=\left(E\left\{g_{1}(X)\right\},\ldots,E\left\{g_{m}(X)\right\}\right)$ belongs to the convex hull of the set $G=\left\{\left(g_{1}(x),\ldots,g_{m}(x)\right)\mid x\in\mathbb{R}\right\}$, then Caratheodory’s theorem [14, Theorem 1.3.6] gives that $P$ can be written as a convex combination of at most $m+1$ points in $G$. Thus, one has that there exist $\theta_{i}\geq 0$ and $z_{i}\in\mathbb{R}$ such that

[TABLE]

and

[TABLE]

Taking $Z$ to be the discrete random variable with probability distribution given by

[TABLE]

one has that equations (6) and (7) are exactly the equations given by (5). ∎

It is important to note that Theorem 1 only states the existence of sigma points satisfying (5), but it does not give any conclusion about uniqueness. In fact, many solutions are possible, and the choice of an adequate set of sigma points has been investigated thoroughly in the literature [6, 15].

Remark 1*.*

Specifically for the first moments of $X$ focused in this paper, i.e. the case that $g_{k}(x)=x^{k}$, $k=0,\ldots,m+1$, one can also see (4) as a Gaussian quadrature integration scheme with $p_{X}$ being the weighting function and $z_{i}$ and $\omega_{i}:=p_{Z}(z_{i})$, $i=1,\ldots,m+1$, being respectively the nodes and their weights in the quadrature formula. Thus, depending on the probability density function $p_{X}$, the sigma points can be readily calculated as the roots of an orthogonal polynomial [16, Section 4.6]. For instance, if $p_{X}$ is a normal distribution, then the sigma points are the roots of a Hermite polynomial.

One can note that if the only intention were to find some discrete random variable such that $E\{f(X)\}=E\{f(Z)\}$, it is always possible to find a $Z$ variable with two sigma points. In fact, in Theorem 1, consider $m=1$ and $f=g_{1}$. However, this would imply the knowledge of the function $f$. In the UT reasoning, it is sought a greater number of sigma points in order that the approximation $E\{f(X)\}=E\{f(Z)\}$ be valid for a greater number of functions $f$. In [17] it is shown that the approximation is good for any function $f$ which can be well approximated by its $m$-order Taylor representation.

On top of that, the larger is $m$, the more precise is the approximation for $E\{f(X)\}$: an estimate for the error $\int_{a}^{b}p_{X}(x)f(x)\,dx-\sum_{i=1}^{m}\omega_{i}f(z_{i})$ is given by

[TABLE]

where $\xi\in\left(a,b\right)\subset\mathbb{R}$ and $h_{m}$ is the associated monic orthogonal polynomial of degree $m$ associated to $p_{X}$ [18, Theorem 3.6.24]. As a matter of fact, the Gauss quadrature computed integral for $E\{f(X)\}$ is exact for all $f(x)$ that are polynomials of degree less or equal than $2m-1$ [18, pp. 172-175].

Since Theorem 1 takes in account generic functions $g_{k}$, one can analyze very general moment settings (for instance, the case of fractional moments is worked in [19]).

Summing up, one can note that the basic assumption for the UT theory until now is that the moments are precisely known. However, this assumption can be too strong in practical situations where the moments are estimated from experiments and then, their values are only known to be in some intervals. In this work, we deal with the scenario of unknown moments and propose to calculate a robust set of sigma-points in the sense of minimizing the worst possible error between this robust choice of sigma-points and the sigma-points computed by using the real values of the moments.

3 Robust UT

To motivate the proposed robust UT consider a normally distributed random variable $X$ with mean $\mu$ and variance $V$. In the case where $\mu$ and $V$ are given and $m=2$, it is known that a set of sigma points are given by [6]

[TABLE]

with its weights given respectively by $\omega_{1}=\frac{1}{6}$, $\omega_{2}=\frac{2}{3}$ and $\omega_{3}=\frac{1}{6}$.

Suppose now that the value of $V$ is not exactly known, but an upper bound $\overline{V}$ and a lower bound $\underline{V}$ for $V$ is known, i.e. $V\in\left[\underline{V},\overline{V}\right]$. A naive approach for choosing the sigma points would be to use the mean value between $\overline{V}$ and $\underline{V}$ for $V$ in (8). In this case, the sigma points $z_{1}$ and $z_{3}$ would be given by

[TABLE]

and $z_{2}$, and the weights $\omega_{1}$, $\omega_{2}$ and $\omega_{3}$ unchanged.

As Figure 1 illustrates, this can be a pessimistic choice for the sigma points, since if, for example, the real value of $V$ is in fact nearer of $\bar{V}$, a better choice of sigma points would be nearer to the point $\left(z_{1},z_{3}\right)=\left(\mu-\sqrt{3\overline{V}},\mu+\sqrt{3\overline{V}}\right)$. In fact, as can be seen in Figure 1, a preferable choice for the sigma points would be at the center of the region:

[TABLE]

This choice is precisely the Chebychev center222There are two non-equivalent definitions of Chebychev center of a bounded set with non-empty interior: the first definition is the center of the minimal-radius ball enclosing this set, and the second is the center of the largest incribed ball in this set [20]. In this paper, only the first definition will be used. of the set of possibles sigma points given that $V\in\left[\underline{V},\overline{V}\right]$. It has the property of having the minimum worst possible error between the chosen sigma point set and the sigma point set corresponding to the true value of $V$.

Consider now a general scenario where some of the moments are not known, but upper and lower bounds for them are. In this case, the set of possible sigma points and their weights would be in the semialgebraic set $S\subset\mathbb{R}^{2m+2}$ of elements $x=(z_{1},\ldots,z_{m+1},\omega_{1},\ldots,\omega_{m+1})$, defined as solution of the system

[TABLE]

where $\mathcal{K}_{1}$ (resp. $\mathcal{K}_{2}$) is the set of indexes $k$ for which the $k$-moment is known (resp. unknown), $\mathcal{K}_{1}\cup\mathcal{K}_{2}=\{1,\ldots,m\}$ and $u_{k_{2}}$ and $\ell_{k_{2}}$ are respectively known upper and lower bounds for the $k_{2}$-moment.

Although Theorem 1 guarantees a solution for (9), this solution may not be unique. On that account, there is more than one set of sigma points and corresponding weights, and it is natural to ask which is a good choice of sigma points and weights. Based on the previous example, a so-called robust choice would be the Chebychev center of the semialgebraic set defined by (9) since this choice minimizes the worst possible error between the chosen sigma point set and the sigma point set corresponding to the possible true values of the moments of $X$.

Different from the previous example, however, in this generic scenario it may be the case that an analytic formula to express the sigma points as the function of the moments of the random variable $X$ is not available, and the following optimization problem must be solved to find the Chebychev center:

[TABLE]

where $S$ is the solution set for (9) and $x=(z_{1},\ldots,z_{m+1},\omega_{1},\ldots,\omega_{m+1})$.

It is important to note that the optimization problem in (10) is not always guaranteed to have a solution, since the set $S$ may not be bounded (if, for example, $\omega_{1}$ is zero, then $z_{1}$ can take any value). Nevertheless, it still makes sense to use the Chebychev center of a large bounded subset of $S$ to try to choose a set of sigma points minimizing the worst possible estimation error. In fact, for any (bounded or unbounded) set $S$, one can take an $\varepsilon>0$ and consider the Chebychev center of the semi-algebraic set $S_{\varepsilon}\subseteq S$ defined as the set of points

[TABLE]

satisfying

[TABLE]

As $\varepsilon$ decreases, $S_{\varepsilon}$ covers a larger part of $S$. Fig. 2 illustrates how the set $S_{\varepsilon}$ approximates the set $S$ by choosing a sufficiently small $\varepsilon>0$.

While the Chebychev center of $S$ may not exist, the next theorem assures the existence of the Chebychev center of $S_{\varepsilon}$ for any $\varepsilon>0$.

Theorem 2**.**

If $m\geq 2$ and $\varepsilon>0$, then the optimization problem

[TABLE]

has a solution.

Proof.

It suffices to prove that the set $S_{\varepsilon}$ is bounded. The variables $\omega_{i}$, $i=1,\ldots,m+1$ are all inside the $m+1$ dimensional simplex and thus bounded. Since

[TABLE]

and $\omega_{i}>0$, it is impossible for any variable $z_{i}$ to grow without bound. ∎

Remark 2*.*

It is interesting to note that alternatively, instead of the compact (closed and bounded) set $S_{\varepsilon}$ one could consider the set $\hat{S}$ defined as the set $S$ up to the first inequalities replaced by a strict inequality (that is, $\omega_{j}\geq 0$ is replaced by $\omega_{j}>0$, $j\in\{1,\ldots,m\}$). This set is indeed bounded as $S_{\varepsilon}$. However, strict inequalities are in general not well-handled by numeric solvers.

The inner maximization problems of (10) and (12) are polynomial optimization problems (POP). Such problems are ubiquitous and are encountered in several fields [21], such as: finance [22, 23, 24], robust and nonlinear control [25, 26], signal processing [27, 28], quantum physics [29] and materials science [30]. It is known that this problem is NP-Hard in general [8], but despite this, it is possible to approximate a POP by convex optimization problems which are computationally feasible using Lasserre’s hierarchy of semidefinite programming relaxations [9, 10]. In fact, by using the GloptiPoly 3 package333Available for download at http://homepages.laas.fr/henrion/software/gloptipoly3/. [31] or the SparsePOP package444Available for download at http://sourceforge.net/projects/sparsepop/. [32] for MATLAB, or ncpol2sdpa library555Available for download at https://gitlab.com/peterwittek/ncpol2sdpa. [33] for Python, these relaxations can be easily implemented and enables the user to transparently construct an increasing sequence of convex LMI relaxations whose optima are guaranteed to converge monotonically to the global optimum of the original non-convex global optimization problem [26]. Moreover, it is possible to numerically certify the global optimality of the problem.

A direct approach to solve (10) would be to use the mentioned Lasserre’s hierarchy to solve the inner maximization problem for a fixed $\hat{x}=\hat{x}_{0}$ and a local optimization algorithm to search which $\hat{x}_{0}$ minimizes (10). This however, besides being a hard numeric problem, would not guarantee a good approximation to the real value of the Chebychev center. In the next section alternative ways to compute (10) with different trade-off between precision of the result and speed of algorithm are proposed.

4 Computation of robust sigma points

4.1 Computation of an outer box to approximate the Chebychev center

The first proposed approach to find an approximate solution to (12) is to compute an outer-bounding box $B$ of $S_{\varepsilon}$ and approximate the Chebychev center of $S_{\varepsilon}$ by the Chebychev center of $B$. By Theorem 2, for $m\geq 2$ the constraint set $S_{\varepsilon}$ is bounded. Thus, each one of the following polynomial optimization problems on $x=(x_{1},\ldots,x_{2m})\in S_{\varepsilon}$ has solution

[TABLE]

and can be used to construct an outer box

[TABLE]

such that $S_{\varepsilon}\subset B$. The next theorem gives an estimate of the error of the outer-bounding box approximation.

Theorem 3**.**

Let $c_{B}$ be the center of Chebychev of $B$ and let $c_{S}$ be the center of Chebychev of $S_{\varepsilon}$. If $d\coloneqq\|c_{B}-c_{S}\|$ is the defined as the distance between the centers of Chebychev, then

[TABLE]

where $\textrm{diam}(B)$ is the diameter of $B$, that is, the least upper bound of the set of all distances between pairs of points in $B$.

Proof.

Suppose that $d>\textrm{diam}(B)/2$. Thus $c_{S}\not\in B$, and since $c_{S}$ is the center of Chebychev of $S_{\varepsilon}$, one has that $c_{S}$ is in the convex hull of $S_{\varepsilon}$. But this implies that $c_{S}$ is in $B$, and the result follows by contradiction. ∎

4.2 Polynomial optimization program to compute minimum enclosing ball

Another approximate solution to problem (12) can be found by using the two-stage approach of [34]. Since the Chebychev center of $S_{\varepsilon}$ is always inside of the outer box $B$ computed in Section 4.1, problem (10) is equivalent to

[TABLE]

In the first stage, the inner optimization problem of (14), namely,

[TABLE]

is approximated by a polynomial function $\tilde{J}_{\tau}(\hat{x})$ of degree $2\tau$ using the semidefinite optimization program outlined in [34]. Then, in the second stage, the outer minimization problem is replaced with the polynomial optimization problem

[TABLE]

which can be solved using the Lasserre’s hierarchy. As the degree $2\tau$ of the approximation polynomial increases, the solution of problem (15) converges to the solution of problem (14) in the sense of [34]. Alternatively, problem (12) can be shown to be equivalent to the following polynomial optimization problem:

[TABLE]

where $\hat{x}$ is the Chebychev center of $S_{\varepsilon}$. In other words, the Chebychev center is the center of the radius of the minimum volume ball that encloses $S_{\varepsilon}$.

Remark 3*.*

It is interesting to note that if the gap between $\ell_{i}$ and $u_{i}$ is not large enough, the resulting semi-definite programs from relaxing the polynomial optimization problems (13) and (16) can be numerically unstable. In this case, however, one may use the naive approach of computing sigma points by using the arithmetic mean of the moments without loss, since the difference between the real Chebychev center and the point computed by this method would be not so great due to the small difference between the upper and lower bounds of the moments.

5 Computation of UT transform

Suppose that $x_{CB}\coloneqq(z_{1}^{C},\ldots,z_{m+1}^{C},\omega_{1}^{C},\ldots,\omega_{m+1}^{C})$ is the Chebychev center of $S$ computed by one of the methods of Section 4. Based on (2), define the function

[TABLE]

As discussed above, $UT_{f}$ is a good approximation for $E\{f(X)\}$ for a sufficiently large $m$. While it is true that $x_{CB}$ approximates of the center of Chebychev of $S$, one also wish to know if $UT_{f}(x_{CB})$ is near to the Chebychev center of $UT_{f}(S)$ (that is, the image of the set $S$ by the function $UT_{f}$). A class of functions $f$ such that $UT_{f}(x_{CB})$ is near the Chebychev center of $UT_{f}(S)$ is given by the functions such that the solution of the following optimization problem

[TABLE]

is sufficiently small. If these values are sufficiently small, the function $UT_{f}$ is a low-distortion geometric embedding [35], and this implies that $UT_{f}(x_{CB})$ is near the Chebychev center of $UT_{f}(S)$. Finally, to compute (17), one can estimate a solution for (17) by uniformly sampling random values $(x_{i},y_{i})$ of $S$ and computing $\max\,D_{i}$, where

[TABLE]

6 Numerical experiments

In this example, the computation of Chebychev center using the methods proposed in this paper will be illustrated. Consider that one desires to compute $2$ sigma points in a scenario where the values of the first and second moment are not precisely known, but it is known that $E\{X\}\in[0,1]$ and that $E\{X^{2}\}\in[0,1]$. In this case, one has (9) with $m=1$, $\mathcal{K}_{1}=\emptyset$, $\mathcal{K}_{2}=\{1,2\}$, and $\ell_{1}=-3$ and $u_{1}=4$, and $\ell_{2}=0$ and $u_{2}=5$.

Naive method: For fixed values of $E\{X\}=\mu$ and $E\{X^{2}\}=V$, one can compute a point inside the set $S$ by using the canonical formula

[TABLE]

A naive choice of sigma points would be the use of the arithmetic mean of the lower and upper bounds for the moments in (19). Using (19) with $\mu=\frac{1}{2}(\ell_{1}+u_{1})=0.5$ and $V=\frac{1}{2}(\ell_{2}+u_{2})=2.5$ results in the sigma points $z_{1}=2$, $z_{2}=-1$ with weights given by $w_{1}=0.5$ and $w_{2}=0.5$.

**Method 1: **By using the method proposed in Section 4.1 it is possible to compute an outer-bounding box approximation to the semi-algebraic set $S_{\varepsilon}$. Using $\varepsilon=0.01$, the computation of an approximate point of the Chebychev center using outer-bounding box approximation results in the sigma points $z_{1}=0$, $z_{2}=0$ with weights respectively given by $w_{1}=0.5$ and $w_{2}=0.5$.

Method 2: Finally, by computing a point using the method proposed in Section 4.2 with $\varepsilon=0.01$ results in the sigma points $z_{1}=-0.0001$, $z_{2}=0.0001$ with weights respectively given by $w_{1}=0.1$ and $w_{2}=0.9$.

In both Method 1 and 2, the MATLAB toolbox Gloptipoly 3 was used to relax the polynomial optimization problems (13) and (16). A relaxation order of $2$ was used and the global optimality of the problems was numerically certified by Gloptipoly3. Moreover, higher relaxation orders could not be used, since using bigger relaxation orders results in semi-definite programs with a large number of variables and constraints, and thus yields an unstable numeric problem.

Figure 3 illustrates the semi-algebraic variety $S_{\varepsilon}$ and the coordinates of the sigma points calculated by the three methods. It can be seen by Figure 3 that the point computed by using the outer-bounding box approximation is the best approximation to the real Chebychev center of the variety. Moreover, it is important to note that the point calculated by Method 2 is very far from the real Chebychev center of the set $S_{\varepsilon}$. This is due to the semi-definite relaxation of (16) being far from the exact solution. In principle, one can increase the relaxation, but this can lead to an unstable numeric problem, which global optimality cannot be more certified.

Finally, to illustrate the robustness of the computed point, $100$ random samples from $E\{X\}$ and $E\{X^{2}\}$ are respectively draw from the uniform distributions $U(\ell_{1},u_{1})$ and $U(\ell_{2},u_{2})$. Solving (18) for $f(x)=\sin(x)$ with $500$ samples gives an estimate for $D$ of $0.9956$, and the posterior expectation in (1) is approximated as

[TABLE]

where $\omega_{1}$, $\omega_{2}$, $z_{1}$ and $z_{2}$ are computed according to the above methods. The mean error $\left|E\{\sin(X)\}-\sum_{i=1}^{2}\omega_{i}\sin(z_{i})\right|$ by the naive method is $0.0339$, while by Method 1 is $3.6932\cdot 10^{-6}$ and by Method 2 is $1.2085\cdot 10^{-4}$.

7 Conclusion

In this paper, it was proposed a way to devise a robust unscented transform by the computation of the Chebychev center of the semialgebraic set defined by the possible choices of sigma points and its weights. Although, in general, this problem is NP-Hard, some methods are proposed in this paper to approximate the solution of the original problem by convex optimization problems. Further works intend to generalize the present work to higher dimensions, as this could enable novel filter designs for multivariate dynamical systems. Another possible extension to this work is to consider adjustments of the proposed algorithm to work with confidence bounds for the moments instead of absolute upper and lower bounds, since the former is more common to be used on interval estimation techniques, such as the bootstrapping method.

Finally, one limitation of the method is that a great number of sigma points can lead to a numerically unstable relaxation of the polynomial optimization problems. Nevertheless, recent advances in polynomial optimization techniques such as novel relaxation hierarchies based on linear programming [36] could help to scale the presented approaches to a higher number of sigma points.

Acknowledgments

The authors would like to thanks Dario Piga and Henrique M. Menegaz for the worthy suggestions and discussions relevant to this paper. We would also like to thanks the Brazilian agencies CNPq and CAPES which partially supported this work.