{'fraction_non_alphanumeric': 0.041681530060642645, 'fraction_numerical': 0.009000356721455424, 'mean_word_length': 4.422407379854189, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Variational Inference MPC using Tsallis Divergence

In this paper, we provide a generalized framework for Variational Inference-Stochastic Optimal Control by using thenon-extensive Tsallis divergence. By incorporating the deformed exponential function into the optimality likelihood function, a novel Tsallis Variational Inference-Model Predictive Control algorithm is derived, which includes prior works such as Variational Inference-Model Predictive Control, Model Predictive PathIntegral Control, Cross Entropy Method, and Stein VariationalInference Model Predictive Control as special cases. The proposed algorithm allows for effective control of the cost/reward transform and is characterized by superior performance in terms of mean and variance reduction of the associated cost. The aforementioned features are supported by a theoretical and numerical analysis on the level of risk sensitivity of the proposed algorithm as well as simulation experiments on 5 different robotic systems with 3 different policy parameterizations.


I. INTRODUCTION
Variational Inference (VI) is a powerful tool for approximating the posterior distribution of the unobserved random variables [1]. VI recasts the approximation problem as an optimization problem. Instead of directly approximating the target distribution p(z|x) of the latent variable z, VI minimizes the Kullback-Leibler (KL) divergence between a tractable variational distribution q(z) and the target distribution. Due to its faster convergence and comparable performance to Markov Chain Monte Carlo sampling methods, VI has received increasing attention in machine learning and robotics [2,3,4].
VI has been applied to Stochastic Optimal Control (SOC) problems recently. In Okada and Taniguchi [5], the authors formulated the SOC problem as a VI problem by setting the desired policy distribution as the target distribution. The VI-SOC framework works directly in the space of policy distributions instead of specific policy parameterizations in most SOC and Reinforcement Learning (RL) frameworks. This gives rise to a unified derivation for a variety of parametric policy distributions, such as unimodal Gaussian and Gaussian mixture in Okada and Taniguchi [5]. Lambert et al. [6] derived the VI-SOC algorithm for non-parametric policy distribution using Stein variational gradient descent. Apart from the policy distribution, the VI-SOC framework is characterized by two components, the optimality likelihood function and distributional distance metric. The optimality likelihood function measures the likelihood of trajectory samples to be optimal and defines the cost/reward transform to allow for different algorithmic developments. The distributional distance metric is the measure of distance between the variational distribution and target distribution.
In existing VI-SOC works, the KL divergence is used as the distributional distance metric due to its simplicity. On the other hand, recent advances in VI research involve extending the framework to other statistical divergences, such as the α-divergence [7] and χ-divergence [8]. In Wang et al. [9], Regli and Silva [10], the authors proposed variants of the α-divergence to improve the performance and robustness of the inference algorithm. Wan et al. [11] further extended the VI framework to f -divergence, which is a broad statistical divergence family that recovers the KL, α and χ-divergence as special cases.
Tsallis divergence is another generalized divergence rooted in non-extensive statistical mechanics [12]. Tsallis divergence and Tsallis entropy are generalizations of the KL divergence and Shannon entropy respectively. The Tsallis entropy is nonadditive (hence non-extensive) in the sense that for a system composed of (probabilistically) independent subsystems, the total entropy differs from the sum of the entropies of the subsystems [13]. Over the last two decades, an increasing number of complex natural, artificial and social systems have verified the predictions and consequences derived from the Tsallis entropy and divergence [14,15,16,17]. Lee et al. [17] demonstrated that Tsallis entropy regularization leads to improved performance and faster convergence in RL applications.
In this paper, we provide a generalized formulation of the VI-SOC framework using the Tsallis divergence and introduce a novel Model Predictive Control (MPC) algorithm. The main contribution of our work is threefold: • We propose the Tsallis VI-MPC algorithm, which allows for additional control of the shape of the cost transform compared to previous VI-MPC algorithms using KL divergence. • We provide a holistic view of connections between Tsallis VI-SOC and the state-of-the-art Model Predictive Path Integral (MPPI) control, Cross Entropy Method (CEM) and Stochastic Search (SS) methods. • We show, both analytically and numerically, that the proposed Tsallis VI-SOC framework achieves lower cost variance than MPPI and CEM. We further demonstrate the superior performance of Tsallis VI-SOC in mean cost and variance minimization on simulated systems from control theory and robotics for 3 different choices of policy distributions. The rest of this paper is organized as follows: in Section II, we review KL VI-SOC and derive the Tsallis VI-SOC framework. In Section III, we discuss the connections between Tsallis VI-SOC and related works. A reparameterization and analysis of the framework is included in Section IV, and we propose the novel Tsallis VI-MPC algorithm in Section V. We showcase the performance of the proposed algorithm against related works in Section VI and conclude the paper in Section VII.

II. TSALLIS VARIATIONAL INFERENCE-STOCHASTIC OPTIMAL CONTROL
In this section, we review the VI-SOC formulation using KL divergence [5] and derive the novel Tsallis divergence VI-SOC framework.

A. SOC Problem Formulation
In discrete time SOC, we work with trajectories τ := {X, U } of state, X := {x 0 , x 1 , . . . , x T } ∈ R nx×T +1 , and control, U := {u 0 , u 1 , . . . , u T −1 } ∈ R nu×T over a finite time horizon T > 0. The goal in SOC problems is to minimize the expected cost defined by an arbitrary cost function J : R nx×T × R nu×T −1 → R + with initial state distribution p(x 0 ) and state transition probability p(x t+1 |x t , u t ) corresponding to stochastic dynamics x t+1 = F (x t , u t , t ) where t ∈ R n is the system stochasticity.

B. KL VI-SOC
We can formulate the SOC problem as an inference problem and apply methods from VI. To apply VI to SOC, we introduce a dummy optimality variable o ∈ {0, 1} with o = 1 indicating that the trajectory τ = {X, U } is optimal. Table I compares the differences in notation between conventional VI methods and SOC formulated as a VI problem.
The objective of VI-SOC is to sample from the posterior distribution Here p(x t+1 |x t , u t ) is the state transition probability, p(x 0 ) is the initial state distribution and p(u t ) is some prior control distribution (e.g. zero mean Gaussian or uniform distribution). For simplicity, we use o to indicate o = 1 and o for o = 0 from here on. We can now formulate the VI-SOC objective as minimizing the distance between a controlled distribution q(τ ) and the target distribution p(τ |o) where the constant p(o) is dropped. The first term in the objective measures the likelihood of a trajectory being optimal while the second term serves as regularization. Splitting the expectation in the first term, we get where p(U ) and q(U ) represent T −1 t=0 p(u t ) and T −1 t=0 q(u t ) respectively due to independence and Optimality Likelihood: The optimality likelihood in Eq. (3) can be parameterized by a non-increasing function of the trajectory cost J(X, U ) as p(o|τ ) := f (J(X, U )). The monotonicity requirement ensures that trajectories incurring higher costs are always less likely to be optimal. Common choices of f include f (x) = exp(−x) and f (x) = 1 {x≤γ} . In this paper, we choose f (x) = exp(−x) such that log p(o|τ ) = −J(X, U ). To avoid excessive notation, we abuse the notation to define J(U ) := E p(X|U ) [J(X, U )] and the N -sample empirical approximation of the expectation J = which can be interpreted as an application of the famous maximum entropy principle.

C. Tsallis VI-SOC
In this subsection, we use the Tsallis divergence as the regularization function and derive the Tsallis-VI-SOC algorithm. First, we define the deformed logarithm and exponential as where (·) + := max(0, ·) and r > 0. Using log r and exp r we can define the Tsallis entropy and the corresponding Tsallis divergence [12] as and D r (q(z) p(z)) :=E q log r q(z) p(z) Note that as r → 1, log r → log, exp r → exp, D r → D KL and S r → S, where S is the Shannon entropy, recovering the KL VI-SOC framework.
Versions of Tsallis Statistics: In our definition of the deformed logarithm and exponential in Eqs. (5) and (6), we use the variable r instead of the q used in most literature to avoid assigning multiple meanings to q. Also, our definition of log r and exp r differ from its original definitions [18], but the original can be recovered with r = 2 − r, where r corresponding to the value used in [18]. Additionally, there are multiple formulations of the Tsallis entropy differing in their definitions of the internal energy and how expectations are taken [12,19,20]. We have chosen to use the formulation from [12] due to its simplicity in computing the expectation. However, as shown in [21], they are equivalent and can be recovered from each other via a change of variables.
We can now define a new objective by replacing the KL divergence regularizer in Eq. (4) with the Tsallis divergence and introducing a parameter λ that multiplies the optimality likelihood to make the regularization strength tunable: The controlled distribution has to satisfy the additional constraint of q(U ) dU = 1. The optimal policy distribution q * can be explicitly solved for with the form whereλ = α(r − 1)λ with α being the Lagrange multiplier for the constraint that q * integrates to 1. A detailed derivation can be found in Section SM-1 of the supplementary materials. We can now use Eq. (10) to obtain the optimal control distribution by transforming the prior distribution.

D. Update Laws
In general, it is computationally inefficient to sample from q * (U ) via Eq. (10) directly. Instead, we can approximate q * (U ) by some policy π(U ) lying in a class Π of tractable distributions and solve for an iterative update law by minimizing the KL divergence between π(U ) and q * (U ): However, because one can only evaluate q * (U ) at a finite number of points {U (n) } N n=1 , we instead approximate q * (U ) by the empirical distributionq * (U ) with weights w (n) : .
We now solve Eq. (11) for 3 different policy classes: unimodal Gaussian, Gaussian mixture and a nonparametric policy corresponding to Stein Variational Gradient Descent (SVGD) [6,22]. The full derivations for each can be found in the supplementary material in Section SM-2.
Unimodal Gaussian: For a unimodal Gaussian policy distribution with parameters Θ := {(µ t , Σ t )} T −1 t=0 , the update laws for the k + 1th iteration take the form of Gaussian Mixture: Alternatively, the policy distribution can be an L-mode mixture of Gaussian distribution with parameters where φ l is the mixture weight for the lth component. Although it is not possible to directly solve Eq. (11) in this case, we draw from Expectation Maximization (EM) to derive an iterative update scheme for the k + 1th iteration with the form where Stein Variational Policy: The policy can also be a nonparametric distribution approximated by a set of particles Θ := {θ l } L l=1 for some parametrized policyπ(U ; θ). In [6],π is taken to be a unimodal Gaussian with fixed variance, where θ ∈ R nx×(T −1) corresponds to the mean. The update law of each Stein particle for the k + 1th iteration has the form where the N is chosen such that N = LS,k is a kernel function, and w (l,s) = w (m+L(l−1)) , where S denotes the number of rollouts for each of the L particle. As noted in [23], SVGD becomes less effective as the dimensionality of the particles increases due to the inverse relationship between the repulsion force in the update law and the dimensionality. Hence, we follow [6] in choosing a sum of local kernel functions as our choice ofk.

III. CONNECTIONS TO RELATED WORKS
In this section, we compare VI-SOC with three samplingbased methods from optimization, thermodynamics, and information theory that have been applied to stochastic control problems. A comparison of problem formulations and update laws for different approaches with a unimodal Gaussian policy is in Table II.

A. Stochastic Search
SS [24] is a stochastic optimization scheme and has also been applied to the SOC setting [25,26]. SS-SOC parameterizes the control policy with a distribution from the exponential family during problem formulation and optimizes with respect to policy parameters whereas VI-SOC performs optimization at the distribution level. The SS-SOC framework formulates the problem as where S(·) is a monotonically decreasing shape function. Note that in Eq. (24), the expectation taken with respect to p(X|U ) is inside the shape function S as opposed to outside as in Eq. (3). For convex shape/optimality likelihood functions, the objective in VI-SOC corresponds to an upper bound of that in SS-SOC, i.e., A detailed comparison of the two formulations can be found in [5].
The parameter update of SS-SOC has the form where T (·) denotes the sufficient statistic for the corresponding parameter and β is the step size. In the case of a unimodal Gaussian policy, the update in Eq. (26) is equivalent to Eqs. (14) and (15) with S(J) = exp r (−λ −1 J).  28)). Left: r is varied with fixed γ = 1. It is clear that exp r → CEM as r → ∞. Right: γ is adjusted accordingly such that exp r → MPPI as r → 1.

B. Cross Entropy Method
CEM [27] is a widely used algorithm in reinforcement learning and optimal control problems [28,29]. The objective function of CEM minimizes the expected cost E[J]. The policy update law for CEM corresponds to that of SS-SOC in Eq. (26) with shape function S(J) = 1 {J≤γ} where γ is the elite threshold. As will be shown in Section IV-A, CEM is also a special case of the reparameterized Tsallis VI-SOC with r → ∞.

C. Model Predictive Path Integral Control
MPPI is another approach closely related to VI-SOC and SS-SOC [30,31]. The framework solves for the controls by minimizing the KL divergence between a controlled distribution and the optimal control distribution The optimal distribution achieves the free energy lower bound, , which has been shown to be the solution of the Hamilton-Jacobi-Bellman equation [32]. The optimal distribution here is equivalent to Eq. (10) when r → 1. The corresponding update law can also be obtained from the SS-SOC framework with S(J) = exp(−λ −1 J).

IV. ANALYSIS A. Effect of exp r
To facilitate the analysis, we focus on the exp r term in Eq. (10). Reparametrizingλ = (r − 1)γ, we get where γ is now the threshold beyond which the optimality weight is set to 0. The reparameterization adjusts the original parameterλ at every iteration to maintain the same threshold γ, which is a more intuitive parameter. In addition, we have observed that the reparameterized framework is easier to tune and achieves better performance than the original formulation. Therefore, we focus our analysis and simulations on only the reparameterized version hereon after. Note that in practice, it is easier to define an elite fraction, which adjusts γ based on the scale of the costs, instead of using γ for easier tuning. Figure 1 illustrates the shapes of the function corresponding to different r and γ values. For J < γ, as r → ∞, which tends to 0 as r → 1. Hence, exp 1 r−1 log 1 − J γ converges to 1 {J=0} as r → 1.

B. Variance Reduction
With the connection to SS, the effect of Tsallis divergence can be analyzed through the equivalent problem formulation in optimization. In Section III-A, it is shown that Tsallis VI-SOC corresponds to an upper bound of SS-SOC with objective The variance reduction effect can be analyzed through the coefficient of Absolute Risk Aversion (ARA) [33] for the optimization problem where Tsallis VI-SOC corresponds to S(J) = exp r (−λ −1 J)).
The ARA coefficient measures a scaled ratio of terms corresponding to the mean and variance terms in the Taylor series expansion of the objective. Negative value of the coefficient corresponds to a risk-averse objective, and positive value of the coefficient corresponds to risk-seeking behavior. 1 The ARA coefficient of Tsallis VI-SOC is .
The ARA coefficients for MPPI and CEM are Since A MPPI is a positive constant, it corresponds to a riskseeking algorithm. For a cost below the elite threshold, J < γ, A CEM = −∞ and A Tsallis ≶ 0 for r ≶ 2. This leads to the Tsallis VI-SOC framework achieving lower variance than MPPI. For the same elite threshold γ and a properly selected r, we hypothesize that Tsallis VI-SOC results in lower mean cost than CEM since CEM penalizes variance infinitely harder than the mean and assigns equal weights to all elite samples. A more detailed analysis and the derivation of ARA coefficient are included in Section SM-3 of the supplementary material.

V. MODEL PREDICTIVE CONTROL ALGORITHM
With the update laws in Section II-D and a choice of the optimality likelihood function f , we propose the novel Tsallis VI-MPC algorithm, summarized in Algorithm 1.
Given the initial state distribution p(x 0 ) and a prior policy distribution, N × M initial states are sampled. Given the initial states, N control trajectories are sampled from the policy distribution. For each control trajectory sample, M state trajectories are propagated for a total of N × M rollouts. The state transitions at each timestep are sampled from the stochastic dynamics F (x, u, ) where corresponds to system stochasticity (zero mean Gaussian ∼ N (0, σ 2 I) used in this paper). The cost of each trajectory is normalized to [0, 1] for numerical stability and easier tuning. Depending on the choice of policy distribution class, the policy parameters are updated based on the update laws in Section II-D. Execute q K (u t ,0 ) 16:  and performing K optimization iterations in the ensuing MPC steps.
Control Selection: With the optimized policy distribution from the VI-SOC framework, the control to be executed on the real system is selected differently based on the choice of policy class. For the unimodal Gaussian policy, the mean of the distribution is used. For the Gaussian mixture policy, the mean of the model with the highest mixture weight is executed. In terms of the Stein policy, the Stein particle with the highest weight is used.
Receding Horizon: After the control execution, the policy distribution is shifted to warm start the optimization at the next MPC timestep. Let θ be the next iteration's starting sequence and θ be the current iteration's sequence and set θ t = θ t+1 for t = 0, ..., T −2. Finally, set the last item in the new sequence as θ T −1 = θ T −1 . This is known as the receding horizon technique in MPC. Note that in Algorithm 1, Θ = {θ 0 , · · · , θ T −1 }.

VI. SIMULATIONS
In this section, we compare the proposed Tsallis VI-MPC algorithm against MPPI and CEM, which represent state-ofthe-art sampling-based SOC algorithms. Since in practice, the shape functions used in SS-SOC corresponds to the ones which result in MPPI and CEM, we have chosen to only compare to MPPI and CEM. We first validate our analysis on a simple numerical experiment. We then showcase the scalability and performance of the Tsallis VI-MPC algorithm on 2D point mass, quadcopter, ant, manipulator and humanoid systems in simulation under the 3 aforementioned policy distributions.

A. Numerical Analysis
We verify the analysis in Section IV-B that Tsallis VI-SOC results in greater variance reduction by comparing the exp r cost transform with the CEM and MPPI cost transforms on the following simple single-stage stochastic optimization problem: where ξ ∼ N (0, 1), erf is the corresponding error function, the constants λ and σ are chosen to be λ = 0.2, σ = 2.5, and c and d are chosen such that the noiseless objective function is normalized between 0 and 1 for u ∈ [−5, 5]. Figure 3 plots the objective function near its minimum and the noisy objective values.
To ensure that only the weight computation of the three methods is tested, we sample 64 instances of u uniformly from [−5, 5] and use the same set of samples for all three methods. We use the update law corresponding to the unimodal Gaussian with fixed variance for each algorithm and compare the resulting mean and standard deviation of the cost obtained from each method after a single optimization iteration. A grid search is performed over 4096 different hyperparameters, and the optimization metric is computed as the cost averaged over 4096 random seeds. The results for the best performing hyperparameters are shown in Table III and agree with our intuition that the objective function of Tsallis VI-SOC should result in a lower variance compared to MPPI and a lower mean compared to CEM. In addition, in Fig. 2, we observe that mean cost increases slower as the elite fraction decreases from the optimal value than when increases, while the opposite is true for cost standard deviation.

B. Controls and Robotics Systems
The dynamics for the planar navigation and quadcopter tasks are solved via an Euler discretization. Details of the dynamics can be found in Section SM-4. For the manipulator, ant, and humanoid tasks, we use the GPU-accelerated Isaac-Gym [34] to sample trajectory rollouts in parallel. Additional system stochasticity is injected to each system through the controls channel such that F (x t , u t , t ) = F (x t , u t + t ).
The hyperparameters and system configurations for all simulations are included in Section SM-4. To ensure fair comparisons, all hyperparameters for each method are tuned using a combination of the TPE algorithm [35] from the Neural Network Intelligence (NNI) AutoML framework and hand tuning.
1) Planar Navigation: We first test the different algorithms on a point-mass planar navigation problem. The task is for the point-mass with double-integrator stochastic dynamics to navigate through an obstacle field to reach the target location. The dynamics and obstacle field are set up the same way as [6]. If a crash occurs, a crash cost is incurred and no further movement is allowed.
2) Quadcopter: We also set up a quadcopter 3-D navigation task. The task is for the quadcopter to reach a target location while avoiding obstacles. The quadcopter dynamics are taken from [36]. The quadcopter task is similiar to the planar navigation task, where the system is expected to fly through a randomly generated forest to reach the target location. If a crash occurs, a crash cost is incurred and no further movement is allowed.
3) Franka Manipulator: We next test on the Franka manipulator modified to have 7-DOF by removing the last joint and fixing the fingers. The objective of this task is to move the end effector around the obstacles to the goal. An illustration of the task is in Fig. 4. It is worth noting that no crash cost is in place for the Franka manipulator. Instead, contact is included as a part of the simulation dynamics.

4) Ant:
We also consider the task of locomotion. We test on the ant system, which has 29 state dimensions, 8 control dimensions, and is a widely used testbed for RL algorithms. 5) Humanoid: Finally, we test our approach on the complex humanoid locomotion task. The humanoid system has 56 state dimensions, 21 control dimensions, and has very unstable dynamics.

C. Discussion
In Table IV, we can see the results of each experiment summarized between each robotic system, policy parameterization, IV: Comparisons of mean and standard deviation of cost against MPPI and CEM on different systems and policy classes. The policy classes are defined as Unimodal Gaussian (UG), Gaussian Mixture (GM), and Stein (S). Note that the negative of the reward is used for the locomotion tasks (ant and humanoid). The best mean cost and cost variance for each system-policy distribution pair is boldfaced. The mean and standard deviation reduction percentages are included to the right where positive values correspond to a reduction. First in the planar navigation case, we see that Tsallis VI-MPC outperforms the other algorithms in almost all policy parameterizations. During testing, there was high variability in the trajectories computed by each algorithm. The best performing algorithms would have a large increase in velocity towards the goal state, and would have enough variation in sampled trajectories to "see" obstacles via large costs in order to avoid them. When comparing the trajectories computed through the GM policy parameterization we tend to see a reduced ability to stabilize at the goal resulting in larger costs overall when compared to other parameterizations. For these systems, the Stein policy results in longer trajectories to the goal. Next in the Quadcopter experiments, we observe that the Tsallis VI-MPC outperforms both MPPI and CEM in mean cost for the unimodal Gaussian and Stein policies. When looking at the high dimensional simulation environments, we see the same trend of Tsallis VI-MPC reducing the mean and variance in comparison to CEM and MPPI holds even for complex systems with contact dynamics.
The performance of the Tsallis VI-MPC can be attributed to the ability of the algorithm to sample policies that are low cost (via the elite fraction), but then improve beyond the CEM by performing the cost weighted averaging similar to MPPI. These characteristics of Tsallis VI-MPC are additionally heavily related to the choice of hyperparameters. From hand tuning, we observe that reducing the elite fraction generally results in lower mean cost but higher standard deviation and vice versa. As the cost transform approaches that of CEM or MPPI, the mean and standard deviation approaches their corresponding value. This verifies that the Tsallis VI-SOC framework can be thought of as an interpolation between CEM and MPPI to a certain degree. The best performing configuration is usually somewhere in between these two extremes.
Finally note that even though the Tsallis VI-SOC framework is typically used over multiple iterations in optimization schemes where the parameters can evolve as the optimization progresses, we have shown the benefits even in MPC mode, where a single iteration of optimization is performed, reiterating the empirical result that a single step has higher reduction in mean cost and variance compared to traditional MPC methods.

VII. CONCLUSION
We present a generalized Variational Inference-Stochastic Optimal Control framework using Tsallis divergence, which allows for additional control of the cost/reward transform and results in lower cost/reward variance. We provide a unifying study of the connections between Tsallis VI-SOC, MPPI, CEM, and SS methods. The performance and variance reduction benefits of the proposed Tsallis VI-SOC framework is verified analytically and numerically. We further showcase advantages of the Tsallis VI-MPC algorithm against MPPI and CEM on 5 different systems with 3 different policy distributions. We leave 2 extensions of this work as future research directions: 1. Implementation of the proposed algorithm on real systems; 2. Comparison against VI-MPC algorithms using other generalized divergences.

SM-1. DERIVATION OF OPTIMAL DISTRIBUTION
To solve for the optimal distribution, we first formulate the Lagrangian as which can be solved by setting ∂L ∂q = 0 as (40) We now use the constraint that q(U ) dU = 1 by integrating both sides to get Plugging the normalizing constants back in we can get the optimal policy q * to be

SM-2. DERIVATION OF UPDATE LAWS
For each of the following update laws, we make use of the following equality: We can drop the first term as it doesn't relate to the arg min π(U )∈Π : w (n) log π(U (n) ). (50)

A. Unimodal Gaussian
For the case of a policy class of unimodal Gaussian distributions where the controls u t ∈ R nu of each timestep are independent, the p.d.f. for each timestep's control has the form where we have that Minimizing the KL divergence between π(U ) andq * (U ) and using Eq. (46), we have π k+1 (U ) = arg min Hence, we have that Using Eqs. (52) and (53) then results in the update laws

B. Mixture of Gaussian
For a Gaussian mixture model π(U ; θ) with L components with parameters θ := {θ l } L l=1 and θ l : Directly trying to minimize the KL divergence between π(U ; θ) andq * (U ) and using Eq. (46) gives us However, it is not straightforward to solve Eq. (62) directly. Instead, we can use the expectation maximization algorithm and introduce the latent variables {Z n } N n=1 , Z n ∈ [L] which form a categorical distribution q(Z) over each of the L mixtures and denotes which mixture each sample U (n) was sampled from. Then, we can write log π(U ; θ) as log π(U ; θ) = q(z) log π(U ; θ)dz (63) where F (q, θ) is the ELBO, and provides a lower bound for the log likelihood log p(U ; θ). Now, to optimize F (q, θ), we estimate q(z) via p(z|U ; θ k ), where θ k is the previous iteration's parameters. Then, the ELBO becomes Since the second term is a function of θ k and thus is not dependent on θ, it suffices to optimize Q(θ, θ k ). To do so, the EM algorithm separates this into two steps: 1) E-step: Compute p(z|U ; θ k ).
In our case, we have that Q(θ, θ k ) = p(z|U ; θ k ) log p(U, z; θ)dz (72) where we have defined η l (u (n) t ; θ k ) for simplicity. To compute this, note that To solve for the optimal θ * in the M-step, we apply the first order conditions to µ l,t and Σ l,t to obtain where we have defined N l,t := N n=0 η l (u (n) t )w (n) for convenience. Since φ * l has the constraint L l=1 φ * l = 1, we solve for φ * l by applying the Lagrangian multiplier λ, yielding the following Lagrangian: Applying first order conditions for φ l gives To solve for the Lagrange multiplier λ, we plug in the constraint to get Hence, we finally get that where we have defined N l := T −1 t=0 N l,t .

C. Non-parametric Policy via SVGD
We now derive the update law for the choice of a fully non-parametric policy via SVGD. SVGD solves the variational inference problem π * = arg min π∈Π {D KL (π(U ) q * (U ))}.
Hence, one can then compute the optimal distribution π * by iteratively applying φ * to some initial distribution p 0 . However, given that we only have an empirical approximationq * of the optimal distribution q * , the gradient of log q * cannot be easily computed. To solve this, following [6], instead of optimizing directly over the distribution of controls U , one can instead optimize over some distribution g of parameters θ of a parametrized policyπ(U ; θ), such that By doing so, we obtain the following new VI problem and steepest descent direction respectively: As a result, we can now take the gradient of log E U ∼π(·;θ) [q * (U )] as Discretizing the above since we only have a discrete approximationq * of q * , we get that Hence, by letting g be an empirical distribution g(θ) = L l=1 1 {θ l } , we obtain the following update rule for the particles

SM-3. VARIANCE REDUCTION ANALYSIS
The variance reduction analysis is conducted through the connection between Tsallis VI-SOC, MPPI, CEM, and SS-SOC. With the appropriate selection of the shape function, the update laws of Tsallis VI-SOC, MPPI and CEM can be derived from SS-SOC. Associated with the SS-SOC update laws are their corresponding problem formulation in stochastic optimization. For the reparameterized Tsallis VI-SOC framework, its associated problem formulation in optimization is Note that the scenario where J ≥ γ is ignored since the Taylor series expansion and analysis will be around 0 < J < γ. Similarly, the optimization problem formulation corresponding to MPPI is For CEM, the corresponding problem formulation is To analyze degree of risk aversion of each problem formulation, we can look at its Taylor series expansion. For a cost likelihood function c(J θ ) where J θ is a stochastic cost term parametrized by the parameter θ, the Taylor series expansion around E[J θ ] :=J θ has the form: Taking the expectation: On the other hand, if we define the certainty equivalent cost J θ CE such that c(J θ CE ) = E[c(J θ )], we have that Hence, Defining π A :=J θ − J θ CE as the Absolute Risk Premium (which is a negative quantity for cost minimization as opposed to a positive quantity for utility/reward maximization), we now have that where A(·) is the ARA coefficient, defined by The absolute risk premium measures the difference between the mean cost and average cost scaled by the cost transform as a function of the cost function variance. Therefore, the ARA coefficient indicates the degree of risk aversion with respect to the cost function variance. The ARA coefficients for Tsallis VI-SOC and MPPI can be easily computed from equation 107 as Since CEM's problem formulation includes the non-differentiable indicator function, we cannot directly apply the same analysis. However, by taking the indicator function as the limit of the sigmoid function σ we can state the CEM problem formulation as With this, the ARA coefficient can be computed as SM-4. SIMULATION DETAILS A table of the system parameters is shown in Table V. A. Planar Navigation  We consider the task of navigating a 2D double integrator through an obstacle field (see Fig. 5). The state consists of (x t ,ẋ t ), while the control is u t =ẍ t , where x t ,ẋ t ,ẍ t ∈ R 2 represent the position, velocity and acceleration of the system. For the initial position x 0 := ( −9 −9 T , 0 0 T ) and goal position x goal := ( 9 9 T , 0 0 T ), we defined the quadratic cost function: where c crash is the indicator variable for crashing into an obstacle, and (116)

B. Quadrotor
The quadrotor task includes navigating a 3D obstacle course with 35 randomly placed obstacles while trying to reach a target location and hover. The control of the quadrotor is done through angular rates and a thrust command, we assume a low level tracking controller that is common for flight systems. The state of the quadrotor x is composed of the 3 coordinates x, y, z angular position represented in quaternions and linear and angular velocities for a total of 13 states. The location of target is represented in cartersian coordinates by x target and is located at (25,25,5). The initial location x 0 of the quadrotor is at (0,0,5). For this task, we use the following quadratic cost function: where x is the cartesian coordinates of the quadrotor,ẋ is the linear velocity, ω is the angular rates, and c crash is the indicator function that is 1 if the following conditions are met, and 0 otherwise 1) If the vehicle gets within 0.75 of an obstacle.
2) If the vehicle goes above 10 meters or below 0 meters in the z axis.

C. Franka
The Franka manipulator is modified to have 7-DOF by removing the last joint and fixing the fingers in place. For state x := (e,ė) and control u :=q, where e,ė ∈ R 3 are the position and velocity of the end-effector, andq ∈ R 7 are the joint velocities, the objective is to reach the goal position x goal with obstacles between the starting position of the end effector and the goal. For this task, we use the following quadratic cost function:

D. Ant
For the Ant, we take inspiration from the cost function for the HalfCheetah from [5] and modify the cost function to reward forward velocity. For the instantaneous x velocity v t , torso height h t and angle to the vertical axis α, the cost function J is defined as

E. Humanoid
Similar to the Ant, we modify the cost function to reward forward velocity only if the height of the agent's root is not too close or far from the ground and the humanoid is oriented upwards. For the instantaneous x velocity v t , torso height h t and angle to the vertical axis α,

SM-5. ADDITIONAL RESULTS
In this section we include the results of different framework's sensitivity with respect to its hyperparpameters on the planar navigation example with fixed variance unimodal Gaussian policy in Table X. We also include the cumulative cost comparison plots of the simulator systems. All plots in this section show error bars of ±1 standard deviation of the cost.