# KernelSOS for Global Sampling-Based Optimal Control and Estimation via Semidefinite Programming

Antoine Groudiev, Fabian Schramm, Eloïse Berthier, Justin Carpentier,
Frederike Dümbgen

<figure><img src="image_1.png" /></figure>

## To cite this version:

Antoine Groudiev, Fabian Schramm, Eloïse Berthier, Justin Carpentier, Frederike Dümbgen. KernelSOS for Global Sampling-Based Optimal Control and Estimation via Semidefinite Programming. 2025. hal-05184030

HAL Id: hal-05184030

https://inria.hal.science/hal-05184030v1

Preprint submitted on 24 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# KernelSOS for Global Sampling-Based Optimal Control and Estimation
via Semidefinite Programming

Antoine Groudiev<fnref n="1" />, Fabian Schramm<fnref n="1" />,<fnref n="2" />, Éloïse Berthier<fnref n="3" />, Justin Carpentier<fnref n="1" />,<fnref n="2" />, and Frederike Dümbgen<fnref n="1" />,<fnref n="2" />

**Abstract**— Global optimization has gained attraction over the past decades, thanks to the development of both theoretical foundations and efficient numerical routines to cope with optimization problems of various complexities. Among recent methods, Kernel Sum of Squares (KernelSOS) appears as a powerful framework, leveraging the potential of sum of squares methods from the polynomial optimization community with the expressivity of kernel methods widely used in machine learning. This paper applies the kernel sum of squares framework for solving control and estimation problems, which exhibit poor local minima. We demonstrate that KernelSOS performs well on a selection of problems from both domains. In particular, we show that KernelSOS is competitive with other sum of squares approaches on estimation problems, while being applicable to non-polynomial and non-parametric formulations. The sample-based nature of KernelSOS allows us to apply it to trajectory optimization problems with an integrated simulator treated as a black box, both as a standalone method and as a powerful initialization method for local solvers, facilitating the discovery of better solutions.

## I. INTRODUCTION

Many tasks in modern robotics and control applications are phrased as optimization problems, including estimation, calibration, inverse kinematics, inverse dynamics, trajectory optimization, frictional simulation, and policy learning. For many applications, the optimization problems are non-convex, and commonly used local solvers are likely to converge to local, rather than global, optima. Local optima may lead to inefficiencies when they prohibit the discovery of good solutions (e.g., in trajectory optimization and policy learning), or they may represent inaccurate solutions and thus pose safety problems (e.g., in estimation and calibration).

To ensure the quality of solutions for *safe and efficient* operation, certifiably optimal methods have been adopted in robotics over the past decades. These methods employ convex relaxations of the original non-tractable optimization problems to identify certifiably optimal solutions. Two dual views have received traction: the moment hierarchy view [1], which can be seen as a generalization of Shor's relaxation [2], and the Sums of Squares (SOS) view [3].

Both views suffer from limited generalizability as they are restricted to optimization problems in parametric, polynomial form. This restriction is problematic in at least two ways. First, many problems of practical interest are not inherently polynomial. As a result, various workaround techniques

<figure><img src="image_3.png" /><figcaption>Fig. 1: This paper demonstrates the effectiveness of KernelSOS for various optimization problems in robotics. Based on function evaluations only, KernelSOS can handle a larger class of problems than classic sums-of-squares methods, particularly problems with non-polynomial and black-box elements. Being based on semidefinite relaxations, one iteration of KernelSOS only requires the solution of a single semidefinite program.</figcaption></figure>

have been proposed to force a polynomial structure, such as reformulations for distance-based optimization [4], [5] or alternatives for exponential maps in problems involving rotations [6]. While often effective in practice, such reformulations may slightly deviate from the desired cost function and may require extensive domain expertise and creativity to discover. Second, scalability issues may arise when a high degree of the moment hierarchy is required to achieve sufficient *relaxation tightness*, a prerequisite to retrieve certifiable solutions [7], [8].

This paper highlights the potential of a recent method, *KernelSOS* [9], to overcome these limitations. The high-level idea is illustrated in Fig. 1. Intuitively, *KernelSOS* reasons in terms of kernels rather than feature functions. While standard SOS methods typically use finite-dimensional monomial basis functions as features, in *KernelSOS*, the feature functions are implicitly defined by a positive-definite kernel, evaluated at a chosen set of samples. For many choices of kernels (e.g., Gaussian, Laplace, Matern), the obtained feature functions are infinite dimensional and encompass a larger class of functions than finite-degree polynomials, thus leading to generally tighter formulations [9].

Because *KernelSOS* relies solely on function evaluations, it belongs to the class of zeroth-order optimization algorithms. This characteristic enables it to address not only non-polynomial objectives but even problems with no explicit parametric form, making it amenable to data-driven rather than model-based settings. In contrast to most black-box optimizers—which typically proceed through iteratively fitting a surrogate model and minimizing it—one step of *KernelSOS* fits and minimizes the surrogate in one shot by solving a single Semidefinite Program (SDP).

We demonstrate the broad applicability of the method through two prototypical examples. First, we treat range-only (RO) localization, where classical SOS approaches require

<footnote n="1">¹École Normale Supérieure, PSL Research University, Paris, France. firstname.lastname@ens.psl.eu</footnote>

<footnote n="2">²Inria, Département d'informatique de l'ENS, CNRS, PSL Research University, Paris, France. firstname.lastname@inria.fr</footnote>

<footnote n="3">³U2IS, ENSTA, Institut Polytechnique de Paris, Palaiseau, France. firstname.lastname@ensta.fr</footnote>
heuristics [4], [11], or a significant increase in problem dimension [5], to be effective. In contrast, *KernelSOS* matches their performance without relying on such reformulations. Second, in trajectory optimization, we show how easily *KernelSOS* integrates black-box components, such as high-fidelity dynamics simulators, into the solution process, and that it outperforms Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [12] in low-sample regimes. Furthermore, we demonstrate that a local solver initialized using *KernelSOS* uncovers qualitatively superior solutions compared with random initializations.

This paper is organized as follows. After reviewing polynomial and black-box optimization methods for robotics in Sec. II, we formally state in Sec. III the two robotics problems on which we will illustrate the performance of the *KernelSOS* method. We detail the method in Sec. IV and present implementation details in Sec. V. We provide numerical experiments in Sec. VI, and conclude in Sec. VII.

## II. RELATED WORK

### A. Sums of Squares methods in robotics

SOS and moment-hierarchy methods have been used for a variety of tasks in robotics. Within control, this includes region-of-attraction estimation [13], control synthesis using funnel methods [14], globally solving optimal control problems via linear programming (LP) formulations [15], or inverse kinematics [16]. In estimation, Shor's relaxation or the moment hierarchy [1] have received particular attention for problems such as robust estimation [17], rotation averaging [18], or simultaneous localization and mapping [19], [20].

*KernelSOS* can be seen as a generalization of SOS methods (and, by duality, moment hierarchy methods) beyond polynomial basis vectors.<fnref n="1" /> In particular, it extends naturally to infinite-dimensional basis vectors. While still in its early stages, it was already applied to optimal control problems [22], to learn value functions, using an LP formulation similar to that in [15]. Here, we explore a broader and still understudied application of the *KernelSOS* method in robotics: its potential to solve non-polynomial and black-box optimization problems.

### B. Non-parametric optimization in robotics

Non-parametric methods (sometimes called *black-box methods*) do not require models in parametric form but only require feasible samples of an optimization problem. *KernelSOS* thus also belongs to this category of methods, along with more established approaches including, for example, evolutionary algorithms such as CMA-ES, and sequential approaches such as Bayesian Optimization (BO).

The goal of CMA-ES [12] is to find one or multiple samples of sufficiently small function value by evolving a *population* over multiple *generations*. A multivariate normal distribution characterizes the population, and the iterative

optimization procedure alternates between sampling from the current distribution and then updating the distribution based on a subset of the new samples. While it may better explore the space than classic gradient-based methods, CMA-ES is still a fundamentally *local* search method.

BO [23], on the other hand, iterates between sampling an (unknown) function optimally (possibly in batches) by maximizing an acquisition function, and estimating then minimizing a surrogate of the unknown function. While the acquisition function is generally designed to be easier to optimize than the original function, finding its maximizer is still often a non-convex optimization problem that may exhibit local maxima. *KernelSOS*, on the other hand, fits and globally minimizes the chosen surrogate model class in one shot by solving an SDP.

As we will see, *KernelSOS* with restarts is similar to CMA-ES in that it also evolves a population of samples over multiple iterations, but with a global minimization performed per iteration. Furthermore, more complex surrogate models than Gaussian distributions are considered. Indeed, the classes of surrogate models that are fit implicitly within *KernelSOS* are more similar to those used in BO. In this paper, we compare *KernelSOS* to CMA-ES. A fair comparison with BO would require integrating an active sampling strategy into *KernelSOS*, which we leave for future work.

## III. PROBLEM STATEMENTS

We are interested in optimization problems of the form

$$ \min_{x \in \Omega} f(x), \qquad (1) $$

where *f* is either unknown or of a parametric form that is too complex to be treated by standard solvers, and *Ω* is a compact and convex set<fnref n="2" /> on which we can sample *f*. We take the black-box optimization approach, where all information available to the solver is evaluations of the function *f* at a finite number of points *x*ᵢ, *i* ∈ [*N*], with shorthand notation [*N*] := {1, ..., *N*}. In particular, we consider the following two example problems in this paper.

### A. Range-only localization

First, we consider a classic state estimation problem: RO localization. It can be formulated as

$$ \min_{x} \underbrace{\sum_{i=1}^{m} \frac{1}{\sigma_{i}^{2}} (d_{i} - \|x - a_{i}\|_{2})^{2}}_{=: f_{RO-non-sq}(x)}, \qquad (2) $$

where *x* is a robot's position, *d*ᵢ are distance measurements to known anchors *a*ᵢ, *i* = 1, ..., *m* (e.g., from time-of-flight measurements), and the *σ*ᵢ<fnref n="2" /> are chosen weights. The solution corresponds to the maximum-likelihood estimator under the assumption that the distance measurements are Gaussian distributed with the Euclidean distances of mean ||*x* − *a*ᵢ||₂ and variance *σ*ᵢ<fnref n="2" />. Note that *f*<sub>RO</sub> is not polynomial, so the

<footnote n="1">¹Although an often overlooked fact is that the theory of SOS is also not necessarily limited to polynomials [21].</footnote>

<footnote n="2">²While the original paper [9] has more general assumptions, we restrict ourselves to such sets for the implementation; all search spaces considered in the following satisfy these assumptions.</footnote>
standard SOS method does not apply. Two known remedies exist.

First, we can square the measurements $d_i$, in which case the function $f_{RO-sq}$ becomes quartic:

$$ \min_{\boldsymbol{x}} \underbrace{\sum_{i=1}^{m} \frac{1}{\sigma_i^2} (d_i^2 - \|\boldsymbol{x} - \boldsymbol{a}_i\|_2^2)^2}_{=:f_{RO-sq}(\boldsymbol{x})}. \quad (3) $$

In this case, known global solutions exist [4], [24]. However, as we will show, this may lead to skewed results as it does not correspond to the assumed noise model.

Another approach introduces new unit-vector variables $n_i$ and reformulates (2)to [5]

$$ \min_{\boldsymbol{x}, \boldsymbol{n}_i \in S^d} \sum_{i=1}^{m} \frac{1}{\sigma_i^2} \|\boldsymbol{n}_i d_i - (\boldsymbol{x} - \boldsymbol{a}_i)\|_2^2, \quad (4) $$

where $S^d$ denotes the unit sphere. One can show that the optimal $\hat{n}_i$ are always aligned with $\boldsymbol{x} - \boldsymbol{a}_i$; the two formulations are thus equivalent. We will show in Section VI-B that using this heuristic outperforms squaring the measurements in terms of accuracy, under the Gaussian noise model on distances. However, we need to introduce $md$ new variables, and turn the problem into a non-convex quadratically constrained problem. Thus, this approach scales poorly with problem dimension.

## B. Trajectory optimization

Secondly, we consider the following instance of Trajectory Optimization (TO):

$$ \begin{aligned} & \min_{\boldsymbol{u}_{1:T}} \underbrace{\|\boldsymbol{x}_{T+1}(\boldsymbol{u}_{1:T})\|^2 + \rho \sum_{t=1}^{T} \|\boldsymbol{u}_t\|^2}_{=:f_{TO}(\boldsymbol{u}_{1:T} | \boldsymbol{x}_{\text{start}})} \\ & \text{s.t. } \boldsymbol{x}_{t+1} = g(\boldsymbol{x}_t, \boldsymbol{u}_t), \quad \boldsymbol{x}_1 = \boldsymbol{x}_{\text{start}}, \end{aligned} \quad (5) $$

where the desired goal state is set to the origin without loss of generality, and $\boldsymbol{u}_{1:T} := \{\boldsymbol{u}_1, \dots, \boldsymbol{u}_T\}$ represent the inputs over the time window $T$. The initial state is $\boldsymbol{x}_{\text{start}}$ and the final state $\boldsymbol{x}_{T+1}$ is obtained through roll-outs of a given dynamics function $g$, placing us in a single-shooting setting. In this paper, the roll-outs happen in a simulator, but they could be equivalently carried out on a real-world system. Note that we use a quadratic cost, which is a common choice in the literature, but this is not a requirement for KernelSOS.

## C. Initializing local solvers

While KernelSOS is able to discover solutions across the entire space, it can be inefficient compared to local solvers. In (5), as it only samples a rollout function, it is not able to exploit the structure of the problem beyond the overall regularity of the optimized function. In contrast, first-order methods such as Iterative Linear Quadratic Regulator (iLQR) [25] or Feasibility-driven Differential Dynamic Programming (FDDP) [26] have direct access to the dynamics function and its derivatives, but they are prone to getting stuck in local minima.

We can use KernelSOS to provide a good initial guess for a local solver. This can be done by running KernelSOS for a limited number of restarts, which will yield a good approximation of the solution in the neighborhood of the global minimum. Then, switching to a local solver will refine the solution. This is particularly useful in cases where the problem is too complex for KernelSOS to solve in a reasonable time, but where local solvers are highly dependent on the initial guess. We will demonstrate that a naive approach of initializing the local solver with a sample of the lowest cost from a batch of samples is insufficient, and that the ability of KernelSOS to interpolate between samples is crucial for obtaining a good initial guess.

# IV. KERNELSOS FOR GLOBAL OPTIMIZATION

This section gives a high-level overview of the KernelSOS method, as introduced in [9]. A visualization of the different steps is shown in Fig. 2.

## A. SOS formulation

In the case where $f$ is a polynomial, the standard SOS formulation is a powerful tool to find its global minimizer. Assuming $f$ attains its global minimizer, the following problem is equivalent to (1):

$$ c^* = \max_{c \in \mathbb{R}} c \quad \text{s.t.} \quad \forall \boldsymbol{x} \in \Omega : f(\boldsymbol{x}) - c \ge 0. \quad (6) $$

Problem (6) finds the largest lower bound $c$ of the function $f$. The problem is convex but has infinitely many constraints. Evaluating if a function (here, $f(\boldsymbol{x}) - c$) is non-negative is, in general, NP-hard. However, a sufficient (but generally non-necessary) condition for non-negativity of a polynomial is that it can be written as a sum of other polynomials squared; such polynomials are called SOS. Using this restriction, we can derive the following problem that lower-bounds (6):

$$ c_{\text{d-SOS}}^* = \max_{c} c \quad \text{s.t.} \quad \forall \boldsymbol{x} \in \Omega : f(\boldsymbol{x}) - c \in \Sigma_{2d}, \quad (7) $$

where $\Sigma_{2d}$ is the set of sums-of-squares polynomials of degree $2d$. This set has the favorable property that it is SDP-representable, meaning that we can reformulate (7) as the following equivalent SDP problem:

$$ \begin{aligned} & c_{\text{d-SOS}}^* = \max_{c, B \succeq 0} c \\ & \text{s.t.} \quad \forall \boldsymbol{x} \in \Omega : f(\boldsymbol{x}) - c = \boldsymbol{v}_d(\boldsymbol{x})^\top B \boldsymbol{v}_d(\boldsymbol{x}), \end{aligned} \quad (8) $$

where $\boldsymbol{v}_d(\boldsymbol{x})$ is a chosen basis vector, e.g., the vector of monomials of $\boldsymbol{x}$ up to degree $d$. Because we have replaced Problem (6) by the SDP-representable restriction (8), we have $c_{\text{d-SOS}}^* \le c^*$. Since some non-negative polynomials cannot be expressed as a sum of squares, $c_{\text{d-SOS}}^* = c^*$ does not hold in general. The difference between the optima is also known as the *relaxation gap*, and may be too large for the solution of (8) to be of any use. When it is zero, we have a *tight* relaxation. Under mild assumptions, the gap is guaranteed to approach zero as we increase the degree $d$ [3], [27]. Sometimes, tightness is achieved at a low degree<fnref n="3" />; however, in

<footnote n="3">³More specifically, finite convergence is generic for polynomial optimization problems [27]. However, knowing whether finite convergence will appear for a given instance is NP-hard [28].</footnote>
many scenarios, the problem becomes prohibitively expensive before a sufficient level of tightness is attained.

## B. KernelSOS formulation

The KernelSOS formulation approaches the problem from a new angle. It replaces the handcrafted basis functions, such as the often used monomials $v_d(x)$, by possibly infinite-dimensional functions in a chosen Reproducing Kernel Hilbert Space (RKHS). The class of functions is specified by choosing a specific positive definite kernel (which uniquely defines the RKHS [29]). Choosing an RKHS can be significantly easier than providing the functions in parametric form, and is in particular more amenable to include black-box elements, such as real-world evaluations in a data-driven setting.<fnref n="4" />

Using this kernelized form allows us to formulate the problem only in terms of samples $x_i$, their pairwise kernel evaluations $k(x_i, x_j)$, and the cost evaluations $f(x_i)$. We define the kernel matrix $K = (k(x_i, x_j))_{i,j \in [N]}$, where $k$ is a chosen positive-definite kernel. We denote its Cholesky decomposition by $K = R^\top R$, with $R$ an upper-triangular matrix. The KernelSOS SDP is defined as follows:

$$c_{\text{k-SOS}}^* = \max_{c, \mathbf{B} \succeq 0} c - \lambda \text{Tr}(\mathbf{B}) \\ \text{s.t. } f(x_i) - c = \mathbf{\Phi}_i^\top \mathbf{B} \mathbf{\Phi}_i, \quad i \in [N], \qquad (9)$$

where the $\mathbf{\Phi}_i$ are the columns of $R$, Tr is the trace and $\lambda > 0$ is a regularization parameter related to the finite number of samples $N$ and the smoothness of $f$ with respect to the kernel $k$. In particular, setting $\lambda$ to 0, the optimal solution is $\min_i f(x_i)$, which does not exploit any regularity of $f$.

To make sense of (9), we visualize in Fig. 2 a univariate polynomial example problem. Along with the function to minimize $f$, we plot what we call the surrogate function,<fnref n="5" /> which corresponds to the optimal function in the chosen RKHS that passes through the given samples. For a well-chosen kernel and calibrated hyperparameters, this surrogate function lower-bounds $f$.

The advantage of (9) with respect to (8) is that it exhibits a more principled way to generate basis functions that reflect the chosen function class, and that it allows for a larger class of functions. It also comes with favorable guarantees: as shown in [9], the (9) problem satisfies $c_{\text{k-SOS}}^* \le c^*$ if $f(x) - c$ lies in the chosen RKHS. In that case, the relaxation gap is guaranteed to asymptotically approach zero as $N \to \infty$. The convergence as we increase the degree of the hierarchy, as in SOS methods, is thus replaced by convergence as we increase the number of samples $N$.

## C. Finding the minimizer

Following [9], and denoting by $\alpha_i$ for $i \in [N]$ the dual variables associated to the equality constraints in (9), we can recover a good candidate of $x^*$, the minimizer of $f$, using

<figure><img src="image_4.png" /><figcaption>Fig. 2: One iteration of KernelSOS explained visually. Given samples (green) of an arbitrary (possibly unknown) function (gray), KernelSOS minimizes a kernel-defined surrogate function (red), by solving a simple SDP. See also Fig. 7 for a visualization in two dimensions.</figcaption></figure>

the following formula:

$$x^* = \sum_{i=1}^{N} \alpha_i^* x_i, \qquad (10)$$

where the $\alpha_i^*$ are the optimal dual variables corresponding to the equality constraints in (9). This candidate is exact when using finite polynomial kernels and when there is no duality gap, which can be seen by considering the dual problem, as explained in Appendix I. We found that it gave satisfactory performance, while the alternative method proposed in [9, Sec. 7.1] exhibited poor conditioning.<fnref n="6" />

## V. IMPLEMENTATION DETAILS

**Custom interior-point solver.** Off-the-shelf solvers such as SCS [30] or MOSEK [31] can be used to solve (9), but can lead to a prohibitive cost for problems with a large number of samples (e.g., $N \ge 100$) as they do not exploit the rank-one structure of the constraints.<fnref n="7" /> Therefore, we implement the custom damped Newton algorithm detailed in [9], which only involves iterations on the $N$ dual variables $\alpha_i$. For range-only experiments, we use MOSEK for its efficiency and robustness on small-scale SDP instances. In large-scale trajectory optimization experiments, we solve the SDP using the custom solver, with precision $\epsilon = 10^{-3}$, and a fixed number of 100 Newton steps.

**Restarting strategy.** Increasing the number of samples $N$ usually improves the quality of the surrogate function, leading to a minimizer closer to the global minimum. However, the runtime of a Newton iteration grows cubically with $N$, making it prohibitive to solve the problem for high values of $N$. Furthermore, samples far from the minimizer are not informative.

To improve scalability, we implement the restarting procedure outlined in [9], which iteratively shrinks the region

<footnote n="4">⁴Note that choosing a polynomial kernel of finite degree recovers the classical SOS formulation [9].</footnote>

<footnote n="5">⁵This is called the model in the original work and can be computed using [9, Equation (37)].</footnote>

<footnote n="6">⁶The alternative essentially replaces the maximization of a constant lower-bound $c$, in which the locality of the solution gets lost, by the maximization of a lower-bounding parabola, which, for sufficiently small curvature, will attain its minimum at the minimizer of the original function.</footnote>

<footnote n="7">⁷Each constraint in (9) is of the form $b_i = \text{Tr}(\mathbf{B} \mathbf{A}_i)$ with $\mathbf{A}_i$ a rank-one matrix.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 3: Illustration of the restarting procedure in 2D. The algorithm starts with a large region $B(z_0, r_0)$, and iteratively shrinks it down to $B(z_w, r_w)$.</figcaption></figure>

to sample from, as shown in Fig. 3. The algorithm is given an initial center $z_0$ and a radius $r_0$. At each iteration, the algorithm solves the optimization problem with $N$ samples in the region $B(z_i, r_i)$, where $B(z_i, r_i)$ is the hypercube of radius $r_i$ centered at $z_i$. The algorithm then updates the center and the radius: $z_{i+1}$ is set to the minimizer found at iteration $i$, and $r_{i+1}$ is set to $\gamma r_i$, where $\gamma \in (0, 1)$ is a *decay factor*. This process is repeated a fixed number of times; we denote the number of restarts by $w$ in the following. Note that samples from previous iterations are discarded, and the algorithm starts from scratch at each iteration.

**Choice of kernels.** An appropriate choice of kernel is crucial for the performance of *KernelSOS*, as it should be in correspondence with the expected smoothness of the function to optimize. For the RO problem, the cost function is expected to be very smooth, and we thus choose the Gaussian kernel, defined by $k(x, x') = \exp(-\|x - x'\|^2 / \sigma^2)$ of scale factor $\sigma$. For TO, since the function to approximate is not as smooth, we choose the *Laplace kernel*, also known as *exponential kernel*, with $k(x, x') = \exp(-\|x - x'\| / \sigma)$. As the associated RKHS space of both kernels is infinite-dimensional (unlike polynomial as for standard SOS), this choice allows to leverage the expressivity of *KernelSOS*.

# VI. EXPERIMENTS

## A. Choice of hyperparameters

Hyperparameters of the *KernelSOS* method include the number of samples per restart $N$, the number of restarts $w$, the kernel scale $\sigma$ and the penalization factor $\lambda$. Initial center $z_0$ and radius $r_0$ are also provided to the algorithm to define the initial search space depending on the properties of the problem; usually, $z_0$ is set to 0 and $r_0$ is chosen in accordance with the problem's scale.

In Fig. 4, we vary the parameters $N$ and $w$ to illustrate the trade-off between algorithm performance and computational cost. The results for multiple ($N, w$) couples are aggregated and presented w.r.t. the total number of samples. We observe that both the performance and the stability of the algorithm increase with the number of samples, but the improvement is less significant after a certain number of samples. Likewise,

<figure><img src="image_6.png" /><figcaption>Fig. 4: Analysis of the effect of the number of restarts and the number of samples on the performance of trajectory optimization (TO).</figcaption></figure>

increasing the number of restarts significantly improves the reward but is subject to a diminishing returns phenomenon.

The choice of kernel parameters $\sigma$ and $\lambda$ is problem-dependent, and we use a grid search on a calibration dataset to find the best values for a given problem. The choice of kernel parameters depends on the smoothness of the function being modeled. As the search space is decreased in size, the function appears smoother. For this reason, we use the following simple heuristics: $\sigma_{r+1}$, the kernel scale for the restart $r+1$, is always set to $\gamma\sigma_r$. A more sophisticated approach would also update $\lambda$, which could help improve the algorithm's performance. We did not find this necessary for the considered examples.

## B. Range-only localization

We compare the following formulations and methods for RO localization. The algorithm *SampleSOS* solves the problem using feasible samples, as suggested by [32] for the unconstrained case (3), and by [33] for the constrained case (4).<fnref n="8" /> By *EquationSOS*, we denote the approach to solving the problem outlined, for example, in [4], [11].<fnref n="9" /> We also compare to a naive baseline, which samples the search space on a uniform grid, and chooses the *best sample* as the solution, and a local solver initialized at the ground truth. We introduce the labels *sq* or *non-sq* to denote the cost with squared and non-squared distance measurements, respectively. Note that *KernelSOS* can sample from the original formulation (2) rather than the reformulation (4), which is considerably simpler.

We plot the performance of the different methods, measured by the localization accuracy, over 10 randomly sampled instances, as a function of the standard deviation of Gaussian noise added to the ground-truth distances. Five known anchors

<footnote n="8">⁸Note that although sample-based, these algorithms assume a polynomial parametric form of known degree.</footnote>

<footnote n="9">⁹In these works, the SDP relaxation from the Shor's relaxation of a QCQP formulation of the problem, which can be seen as the dual of the *EquationSOS* view.</footnote>
<figure><img src="image_7.png" /><figcaption>Fig. 5: Distance to ground truth as a function of the noise level in range-only (RO) localization. We compare a local solver initialized at ground truth (local), the classic SOS method (EquationSOS), the (polynomial) sampling-based SOS method by [32], [33] (SampleSOS), our implementation of KernelSOS, and the best-sample baseline. Note that for almost all noise levels, the performance of the (parametric) local, EquationSOS, and SampleSOS solvers is indistinguishable. Notably, the non-parametric KernelSOS method achieves similar performance for a wide range of noise levels, and significantly outperforms the best-sample baseline.</figcaption></figure>

are sampled uniformly at random and rescaled such that their bounding box is 2 × 2 units. The unknown target is sampled uniformly at random within this box. We use the Gaussian kernel with $\sigma = 1.4$ and $\lambda = 10^{-3}$ for the squared cost function, and $\sigma = 1.0$ for the non-squared cost function. These parameters are determined using a grid search on the localization accuracy for a calibration dataset. We use $N = 36$ points throughout all experiments.

The results are shown in Fig. 5. We make the following observations. First, the non-squared formulation results in slightly better accuracy, particularly at high noise levels. This is expected, as it reflects the chosen noise model. KernelSOS, for the chosen number of samples $N$, struggles to realize this better accuracy on the non-squared formulation. This may be a reflection of the relatively localized cost function at the lower noise levels, as shown in the left-most plot of Fig. 7. Such a sparse cost function makes it more challenging to accurately find the solution from only finite samples. However, it is notable that KernelSOS clearly outperforms the best-sample baseline, showing that the interpolation based on the chosen kernel helps retrieve a better solution.

Secondly, we observe that for the higher noise levels of the non-squared formulation and for all noise levels of the squared formulation, there is almost no difference in accuracy between the equation-free KernelSOS method and the other solvers. This is notable, as KernelSOS, conversely to the other methods, does not use any parametric form of the cost, but only a smoothness prior, incorporated by the choice of kernel and parameters $\lambda$ and $\sigma$.

Coming back to Fig. 7, we can observe qualitatively that KernelSOS implicitly fits a relatively accurate surrogate

<figure><img src="image_8.png" /><figcaption>Fig. 6: Runtime analysis for the range-only (RO) problem as a function of the number of samples. We use MOSEK in this problem instance due to ill-conditioning, which leads to higher runtimes compared with the rank-exploiting solver used in the trajectory optimization experiments.</figcaption></figure>

function (middle plots), obtained from evaluating the quadratic form with the optimal matrix $B^*$ post-hoc on a grid, as explained in [9, Eq. 37]. For the non-squared setting on the left, the first step of KernelSOS identifies the correct region to zoom in on, while the second step is necessary to achieve higher accuracy. For the squared setting on the right, the cost is less localized, and a good solution is already found with the first iteration.

The runtime of the different methods is provided in Fig. 6. We confirm that KernelSOS has the same empirical complexity as SampleSOS, since the bottleneck for both methods is the solution of an SDP of the same structure, while the cost for EquationSOS and the well-initialized local solver is significantly lower. We empirically observed that the Hessian becomes singular for the RO localization problems; thus, the low-rank exploiting solver suggested by [9] encounters numerical issues. We resort to MOSEK here, but are working on an extension of the customized solver to deal with ill conditioning. The expectation is that with such a solver, KernelSOS and SampleSOS can become more competitive with EquationSOS. However, we also note that the main computational advantage is expected when EquationSOS needs a high degree in the moment hierarchy, which is not the case for this particular example problem.

### C. Trajectory optimization

We evaluated the performance of KernelSOS for TO on simulations of simple and double pendulum systems, simulated using the Pinocchio library [34], [35]. We consider the swing-up task, hence the desired goal state of $f_{TO}$ in (5) is the upright position with no velocity.

In Fig. 8, we compare the normalized cost of trajectories found by KernelSOS and CMA-ES, as a function of the number of calls to the simulator’s dynamics function. We
<figure><img src="image_9.png" /><figcaption>Fig. 7: Illustration of the surrogate cost and the benefit of restarts, on the range-only problem, using non-squared distances (left half) or squared distances (right half). In each block, we have on the left: real cost; center: surrogate function found by the first *KernelSOS* step; right: surrogate function found by the first restart (second step). The surrogate cost is sampled within the bounding box of the sampling domain for the current step, as indicated by squares of the corresponding color for the real-cost plots. Black crosses represent the known landmarks, black x-marks represent the used samples. Orange and green x-marks denote the minimizers found by *KernelSOS* initially and after restarting, respectively.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8: Comparison of performance of *KernelSOS* and CMA-ES on TO tasks. *KernelSOS* results are averaged over 10 runs of the algorithm to evaluate the influence of the random samples. The shaded region represents the standard deviation for the randomness of the method.</figcaption></figure>

can see that *KernelSOS* outperforms CMA-ES for a low number of samples, while they approximately match for a higher number of samples. This is expected, as the global nature of *KernelSOS* allows it to find good trajectories from the entire space from the start, while CMA-ES tends to explore the search space more incrementally. However, CMA-ES handles different scales of the search space better than *KernelSOS*, which might require a more sophisticated automated adaptation of the kernel parameters than the simple rescaling of σ described in Section VI-A.

<figure><img src="image_11.png" /><figcaption>Fig. 9: Visualization of the trajectories found by *KernelSOS* (left, cost: 0.119) and CMA-ES (right, cost: 0.090) for the double pendulum swing-up task.</figcaption></figure>

We plot in Fig. 9 trajectories found by both methods for

the double pendulum swing-up task. Both methods, despite being fully black-box, can find trajectories that swing up the pendulum. While the final state of the pendulum is further from the goal state than the one found using first-order local methods, as shown in Fig. 11, this shows how using a global optimization method such as *KernelSOS* can find good trajectories despite the unstable nature of the problem. The ability to find relatively good solutions with a small number of samples motivates the use of *KernelSOS* as an initialization method for local optimization methods, as discussed next.

### D. Initializing local solvers

<figure><img src="image_12.png" /><figcaption>Fig. 10: Comparison of the performance of FDDP depending on the initialization. Using *KernelSOS* as an initializer reduces the cost by 12% on average compared to a random initialization. “Iterations” refers to the number of iterations needed by FDDP to converge after initialization. We observe the emergence of three discrete modes, where the best one is only found by *KernelSOS* initialization.</figcaption></figure>

We compare the performance of FDDP [26], implemented in the Aligator [36], [37] library, for different initialization mechanisms. Because it exhibits numerous local minima, we used the double pendulum swing-up task as a benchmark, and compared the following initializations:

* random: a single random input sequence is sampled from the search space;

* best-sample: 100 random input sequences are sampled from the search space, and the one with the lowest cost
<figure><img src="image_13.png" /><figcaption>Trajectory found by random initialization (cost: $6.2 \times 10^{-4}$)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Trajectory found by best sample initialization (cost: $5.5 \times 10^{-4}$)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Trajectory found by *KernelSOS* initialization (cost: $5.0 \times 10^{-4}$)</figcaption></figure>

<caption>Fig. 11: Visualization of the trajectories found by FDDP after initialization with different methods. On the left is the only trajectory found by random initialization; in the middle, the best trajectory found by best-sample initialization; on the right, the best trajectory found by *KernelSOS* initialization.</caption>

is used as an initialization;

* *KernelSOS*: the solver is initialized using *KernelSOS* applied to (5), with 100 samples and no restarts.

Cost performance results of 30 trajectories per initialization method are shown in Fig. 10. We observe that using *KernelSOS* as an initialization reduces the cost by 12% on average compared to a random initialization, and by 10% compared to the one using the best sample. Due to the high dimensionality of the problem, the best-sample initialization only provides a small improvement over the random initialization, whereas *KernelSOS* significantly improves the initialization. Indeed, while the problem has only two degrees of freedom, planning a trajectory requires optimizing the control inputs at each time step, resulting in a high-dimensional search space (here, of size 400). In such a setup, the ability of *KernelSOS* to interpolate the cost function based on a limited number of samples is crucial to find a good initialization.

It is interesting to note that the trajectories found after running FDDP always converge to one of three distinct costs, which correspond to three local minima of the problem. These three modes are visualized in Fig. 11. Notably, *KernelSOS* is the only initialization method that manages to find the lowest-cost solution among the three.

A reasonable question to ask is whether the additional time taken by the initialization step is worth the improvement in performance. Counterintuitively, we observe that the total time taken when initializing with *KernelSOS* is similar to that taken when randomly initializing, despite *KernelSOS* being considerably more costly than random sampling. This is because *KernelSOS* initialization reduces the number of iterations needed by FDDP to converge by 77% on average, as shown in Fig. 12. As the initial guess provided by *KernelSOS* is closer to a local minimum or in a region with better conditioning, FDDP converges faster and requires fewer iterations to converge. The best-sample initialization reduces the number of iterations needed to converge only slightly, while requiring a significant time to evaluate the cost of 100

<figure><img src="image_16.png" /><figcaption>Fig. 12: Comparison of the number of FDDP iterations needed after initialization, and time taken by the full optimization (including initialization). *KernelSOS* initialization reduces the number of iterations by 77% on average, leading to a similar total time as random initialization despite the additional initialization time.</figcaption></figure>

trajectories; as a result, the best-sample method is the slowest of the three. In summary, *KernelSOS* discovers higher-quality solutions in the same amount of time as random sampling methods.

## VII. CONCLUSION AND FUTURE WORK

We presented *KernelSOS*, a strong new contender in the optimization field, which is a sampling-based optimization method that generalizes the SOS framework to non-polynomial, non-parametric problems. Its black-box nature allows it to be applied to a wide range of problems in estimation, control, and beyond. While it can be used as a standalone solver, its demonstrated use as an initialization method enables the combination of the global search capabilities of *KernelSOS* with the efficiency and accuracy of local optimization methods.
We believe that *KernelSOS* has considerable potential for future work, both in terms of theoretical analysis and practical applications.

The method itself could be further extended by improving the sampling strategy, leading to more efficient exploration of the search space. Active sampling methods, as used, for example, in Bayesian optimization, should be investigated. Furthermore, adding the ability to handle constraints could be a promising direction to generalize the method to an even broader class of problems.

On the application side, the effectiveness of *KernelSOS* as an initialization method could be further explored for policy optimization, where it could be used as an initializer or for improved search direction updates in reinforcement learning, potentially leading to faster convergence, better solutions, and enhanced sample efficiency.

## ACKNOWLEDGMENTS

We warmly thank Alessandro Rudi, Francis Bach, and Corbinian Schlosser for many insightful discussions, as well as Lorenzo Shaikewitz for his helpful comments on the manuscript. This work has received support from the French government, managed by the National Research Agency, through the NIMBLE project (ANR-22-CE33-0008) and under the France 2030 program with the references Organic Robotics Program (PEPR O2R) and “PR[A]RIE-PSAI” (ANR-23-IACL-0008). This work was also supported by the European Union through the AGIMUS project (GA no.101070165), and Paris Île-de-France Région has also supported this work in the framework of DIM AI4IDF. Views and opinions expressed are those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

## REFERENCES

[1] J. B. Lasserre, “Global Optimization with Polynomials and the Problem of Moments,” *SIAM Journal on Optimization*, vol. 11, no. 3, pp. 796–817, 2001.

[2] N. Z. Shor, “Quadratic Optimization Problems,” *Soviet Journal of Computer and Systems Sciences*, vol. 25, pp. 1–11, 1987.

[3] P. A. Parrilo, “Semidefinite programming relaxations for semialgebraic problems,” *Mathematical Programming*, vol. 96, no. 2, pp. 293–320, 2003.

[4] A. Beck, P. Stoica, and J. Li, “Exact and Approximate Solutions of Source Localization Problems,” *IEEE Transactions on Signal Processing*, vol. 56, no. 5, pp. 1770–1778, 2008.

[5] T. Halsted and M. Schwager, “The Riemannian Elevator for Certifiable Distance-based Localization,” *Preprint*, 2022.

[6] T. D. Barfoot, C. Holmes, and F. Dümbgen, “Certifiably Optimal Rotation and Pose Estimation Based on the Cayley Map,” *International Journal of Robotics Research*, 2024.

[7] H. Yang, J. Shi, and L. Carlone, “TEASER: Fast and certifiable point cloud registration,” *IEEE Transactions on Robotics*, vol. 32, no. 2, pp. 314–333, 2020.

[8] B. E. Khadir, J. Bernard Lasserre, and V. Sindhwani, “Piecewise-Linear Motion Planning amidst Static, Moving, or Morphing Obstacles,” in *IEEE International Conference on Robotics and Automation (ICRA)*, 2021, pp. 7802–7808.

[9] A. Rudi, U. Marteau-Ferey, and F. Bach, “Finding global minima via kernel approximations,” *Mathematical Programming*, vol. 209, no. 1, pp. 703–784, 2025.

[10] S. Caron, N. Perrin-Gilbert, V. Ledoux, et al., *Upkie open source wheeled biped robot*, https://github.com/upkie/upkie, 2024.

[11] F. Dümbgen, C. Holmes, and T. D. Barfoot, “Safe and Smooth: Certified Continuous-Time Range-Only Localization,” *IEEE Robotics and Automation Letters*, vol. 8, no. 2, pp. 1117–1124, 2023.

[12] N. Hansen, “The CMA Evolution Strategy: A Tutorial,” *arXiv:1604.00772*, 2023.

[13] D. Henrion and M. Korda, “Convex computation of the region of attraction of polynomial control systems,” *IEEE Transactions on Automatic Control*, vol. 59, no. 2, pp. 297–312, 2013.

[14] R. Tedrake, I. R. Manchester, M. Tobenkin, and J. W. Roberts, “LQR-trees: Feedback motion planning via sums-of-squares verification,” *The International Journal of Robotics Research*, vol. 29, no. 8, pp. 1038–1052, 2010.

[15] J. B. Lasserre, D. Henrion, C. Prieur, and E. Trélat, “Nonlinear optimal control via occupation measures and LMI-relaxations,” *SIAM Journal on Control and Optimization*, vol. 47, no. 4, pp. 1643–1666, 2008.

[16] M. Giamou, F. Marić, D. M. Rosen, et al., “Convex Iteration for Distance-Geometric Inverse Kinematics,” *IEEE Robotics and Automation Letters*, vol. 7, no. 2, pp. 1952–1959, 2022.

[17] H. Yang and L. Carlone, “Certifiably Optimal Outlier-Robust Geometric Perception: Semidefinite Relaxations and Scalable Global Optimization,” *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 45, no. 3, pp. 2816–2834, 2023.

[18] A. Eriksson, C. Olsson, F. Kahl, and T.-J. Chin, “Rotation Averaging and Strong Duality,” in *IEEE/CVF Conference on Computer Vision and Pattern Recognition*, 2018, pp. 127–135.

[19] D. M. Rosen, L. Carlone, A. S. Bandeira, and J. J. Leonard, “SE-sync: A certifiably correct algorithm for synchronization over the special Euclidean group,” *International Journal of Robotics Research*, vol. 38, no. 2-3, pp. 95–125, 2019.

[20] L. Carlone and F. Dellaert, “Duality-based verification techniques for 2D SLAM,” in *IEEE International Conference on Robotcis and Automation (ICRA)*, 2015, pp. 4589–4596.

[21] F. Bach and A. Rudi, “Exponential Convergence of Sum-of-Squares Hierarchies for Trigonometric Polyno-
mials,” *SIAM Journal on Optimization*, vol. 33, no. 3, pp. 2137–2159, 2023.

[22] E. Berthier, J. Carpentier, A. Rudi, and F. Bach, “Infinite-dimensional sums-of-squares for optimal control,” in *IEEE Conference on Decision and Control (CDC)*, 2022, pp. 577–582.

[23] H. J. Kushner, “A New Method of Locating the Maximum Point of an Arbitrary Multipeak Curve in the Presence of Noise,” *Journal of Basic Engineering*, vol. 86, no. 1, pp. 97–106, 1964.

[24] J. Fredriksson, V. Larsson, C. Olsson, and F. Kahl, “Optimal relative pose with unknown correspondences,” in *IEEE/CVF Conference on Computer Vision and Pattern Recognition*, 2016.

[25] W. Li and E. Todorov, “Iterative linear quadratic regulator design for nonlinear biological movement systems,” in *First International Conference on Informatics in Control, Automation and Robotics*, SciTePress, 2004.

[26] C. Mastalli, W. Merkt, J. Marti-Saumell, et al., “A feasibility-driven approach to control-limited DDP,” *Auton. Robots*, vol. 46, no. 8, pp. 985–1005, 2022.

[27] J. Nie, “Optimality conditions and finite convergence of Lasserre’s hierarchy,” *Mathematical Programming*, vol. 146, no. 1, pp. 97–121, 2014.

[28] L. F. Vargas, “On the hardness of deciding the finite convergence of Lasserre hierarchies,” *AIMS Numerical Algebra, Control and Optimization*, Jan. 2024.

[29] J.-P. Vert, K. Tsuda, and B. Schölkopf, “A primer on kernel methods,” *Kernel Methods in Computational Biology*, vol. 47, pp. 35–70, 2004.

[30] B. O’Donoghue, E. Chu, N. Parikh, and S. Boyd, “Conic optimization via operator splitting and homogeneous self-dual embedding,” *Journal of Optimization Theory and Applications*, vol. 169, no. 3, pp. 1042–1068, 2016.

[31] MOSEK ApS, MOSEK Optimization Toolbox for Python, Version 10.0, https://www.mosek.com, 2022.

[32] J. Lofberg and P. Parrilo, “From coefficients to samples: A new approach to SOS optimization,” in *IEEE Conference on Decision and Control (CDC)*, 2004, pp. 3154–3159.

[33] D. Cifuentes and P. A. Parrilo, “Sampling algebraic varieties for sum of squares programs,” *SIAM Journal on Optimization*, vol. 27, no. 4, pp. 2381–2404, 2017.

[34] J. Carpentier, F. Valenza, N. Mansard, et al., *Pinocchio: fast forward and inverse dynamics for poly-articulated systems*, https://stack-of-tasks.github.io/pinocchio, 2015–2021.

[35] J. Carpentier, G. Saurel, G. Buondonno, et al., “The Pinocchio C++ library: A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives,” in *IEEE/SICE International Symposium on System Integration (SII)*, 2019, pp. 614–619.

[36] W. Jallet, A. Bambade, S. El Kazdadi, C. Justin, and M. Nicolas, *Aligator*, https://github.com/Simple-Robotics/aligator.

[37] W. Jallet, A. Bambade, E. Arlaud, S. El-Kazdadi, N. Mansard, and J. Carpentier, “ProxDDP: Proximal constrained trajectory optimization,” *IEEE Transactions on Robotics*, vol. 41, pp. 2605–2624, 2025.

# APPENDIX I

## FINDING THE MINIMIZER VIA THE DUAL PROBLEM

In this section, we explain where the solution candidate given by (10) comes from. The dual problem of (9) is:

$$ \begin{array}{ll} \min_{\alpha, M} & \sum_i f(\mathbf{x}_i)\alpha_i \\ \text{s.t.} & \left\{ \begin{array}{l} \sum_i \alpha_i = 1 \\ M := \sum_i \alpha_i v_d(\mathbf{x}_i) v_d(\mathbf{x}_i)^\top \ge 0, \end{array} \right. \end{array} \quad (11) $$

where $M$ is the dual variable associated with $B \succeq 0$, sometimes called the pseudo-moment matrix, and $\alpha \in \mathbb{R}^N$ contains the dual variables associated with the equality constraints. Here, we assume $v_d$ is the monomial basis vector of degree $d \ge 1$. When we have strong duality and a unique minimizer, the dual problem must have a rank-one optimal solution of the form $M^* = v_d(\mathbf{x}^*)v_d(\mathbf{x}^*)^\top$, from which the elements of $\mathbf{x}^*$ can be extracted. The identity (10) emerges since we also have $M^* = \sum_i \alpha_i^* v_d(\mathbf{x}_i^*) v_d(\mathbf{x}_i^*)^\top$ from (11). Note that when there are multiple global minima, the estimate may, in theory, be poor. We did not find this to be an issue for the considered example problems.