Support Sobol sequences sampling, periodic kernel, and improve CI

.github/workflows/python-package-conda.yml

@@ -1,13 +1,10 @@
 name: Python Package using Conda
 
-on: [push]
+on: [push, pull_request]
 
 jobs:
-  build-linux:
+  lint:
     runs-on: ubuntu-latest
-    strategy:
-      max-parallel: 5
-
     steps:
     - uses: actions/checkout@v4
     - name: Set up Python 3.12
@@ -34,3 +31,104 @@ jobs:
       run: |
         conda install pytest
         pytest
+
+  linux:
+    runs-on: ${{ matrix.platform.runner }}
+    strategy:
+      matrix:
+        platform:
+          - runner: ubuntu-22.04
+            target: x86_64
+          - runner: ubuntu-22.04
+            target: x86
+          - runner: ubuntu-22.04
+            target: aarch64
+          - runner: ubuntu-22.04
+            target: armv7
+          - runner: ubuntu-22.04
+            target: s390x
+          - runner: ubuntu-22.04
+            target: ppc64le
+
+    steps:
+    - uses: actions/checkout@v4
+    - name: Set up Python 3.12
+      uses: actions/setup-python@v3
+      with:
+        python-version: '3.12'
+    - name: Add conda to system path
+      run: |
+        # $CONDA is an environment variable pointing to the root of the miniconda directory
+        echo $CONDA/bin >> $GITHUB_PATH
+    - name: Install dependencies
+      run: |
+        conda env create --file ksos_env.yml --name ksos_env
+    - name: Test with pytest
+      env:
+          MOSEKLM_LICENSE_FILE: ${{ secrets.MOSEK_LICENSE }}
+      run: |
+        source $CONDA/etc/profile.d/conda.sh
+        conda activate ksos_env
+        pytest
+
+  windows:
+    runs-on: ${{ matrix.platform.runner }}
+    strategy:
+      matrix:
+        platform:
+          - runner: windows-latest
+            target: x64
+          - runner: windows-latest
+            target: x86
+
+    steps:
+    - uses: actions/checkout@v4
+    - name: Set up Python 3.12
+      uses: actions/setup-python@v3
+      with:
+        python-version: '3.12'
+        architecture: ${{ matrix.platform.target }}
+    - uses: conda-incubator/setup-miniconda@v3
+      with:
+        auto-update-conda: true
+        python-version: '3.12'
+    - name: Install dependencies
+      run: |
+        conda env create --file ksos_env.yml --name ksos_env
+    - name: Test with pytest
+      env:
+          MOSEKLM_LICENSE_FILE: ${{ secrets.MOSEK_LICENSE }}
+      run: |
+        conda activate ksos_env
+        python -m pytest
+
+  macos:
+    runs-on: ${{ matrix.platform.runner }}
+    strategy:
+      matrix:
+        platform:
+          - runner: macos-14
+            target: x86_64
+          - runner: macos-14
+            target: aarch64
+
+    steps:
+    - uses: actions/checkout@v4
+    - name: Set up Python 3.12
+      uses: actions/setup-python@v3
+      with:
+        python-version: '3.12'
+    - uses: conda-incubator/setup-miniconda@v3
+      with:
+        auto-update-conda: true
+        python-version: '3.12'
+    - name: Install dependencies
+      run: |
+        conda env create --file ksos_env.yml --name ksos_env
+    - name: Test with pytest
+      env:
+          MOSEKLM_LICENSE_FILE: ${{ secrets.MOSEK_LICENSE }}
+      run: |
+        source $CONDA/etc/profile.d/conda.sh
+        conda activate ksos_env
+        python -m pytest

CHANGELOG.md

@@ -0,0 +1,17 @@
+# Changelog
+
+This file is used to track changes made to the project over time.
+
+## [Unreleased]
+### Additions
+- Support for Sobol sequences in polynomial problem sampling
+- Support for periodic kernel
+
+## [0.2.2] - 2025-11-14
+### Additions
+- Add penalty parameter for soft-constrained version used with polynomial kernels
+- Create helper functions for monomial feature function creation for polynomial kernel
+
+### Improvements
+- Improve timing calculations
+- Improve efficiency of surrogate function calculation

README.md

@@ -1,12 +1,12 @@
 ![](overview.png)
-# KernelSOS Toolbox
+# Kernel Sum-of-Squares (KernelSOS) Toolbox
 Implementation of the Kernel Sum-of-Squares optimization framework and various related solvers.
 
 The kernelSOS algorithm solves the following general optimization problem:
 ```math
 \min_{x\in\Omega} f(x)
 ```
-where $f$ is any function that can be sampled over the convex bouded domain $\Omega$. Its generality allows it to be applied to a wide range of problems.
+where $f$ is any function that can be sampled over the convex bounded domain $\Omega$. Its generality allows it to be applied to a wide range of problems.
 
 ## Usage
 The main solver is implemented as `ksos_tools.solvers.ksos.solve`. A minimal example of how to use it is as follows (see `example.py`):
@@ -51,6 +51,32 @@ print(f"True solution:  x={-np.pi/2:.4f}, f=-1")
 
 A large variety of options can be passed to the solver, and can be found in the [`ksos.solve`](ksos_tools/solvers/ksos.py) function documentation. Many examples are showcased in the `tests` directory.
 
+## Citation
+If you use this code in your research, please use the following citation:
+```bibtex
+@software{ksos_tools,
+  author = {Groudiev, Antoine and Dümbgen, Frederike},
+  title = {{ksos-tools}: Implementation of the {Kernel Sum-of-Squares} optimization framework and various related solvers},
+  url = {https://github.com/Simple-Robotics/ksos-tools},
+  version = {0.2.2},
+  date = {2025-11-14},
+}
+```
+
+You can also cite the following related paper:
+
+```bibtex
+@misc{groudiev2025ksos,
+    title={Sampling-Based Global Optimal Control and Estimation via Semidefinite Programming}, 
+    author={Groudiev, Antoine and Schramm, Fabian and Berthier, Éloïse and Carpentier, Justin and Dümbgen, Frederike},
+    year={2025},
+    eprint={2507.17572},
+    archivePrefix={arXiv},
+    primaryClass={cs.RO},
+    url={https://arxiv.org/abs/2507.17572}, 
+}
+```
+
 ## Acknowledgments
 
 This project received funding from the French government, managed by the National Research Agency (ANR-22-CE33-0008,PEPR O2R,ANR-23-IACL-0008), and by the European Union (GA no.101207106), the ARTIFACT project (GA no.101165695, GA no.101070165). Paris Île-de-France Région supported this work in the framework of DIM AI4IDF.

example.py

@@ -2,24 +2,24 @@
 import numpy as np
 
 # Standard parameters
-center = [0.0] # center for sampling
-radius = np.pi # radius of sampling
-sampling = "linspace" # use a uniform grid for sampling
-n_samples = 10 # how many samples to use
+center = [0.0]  # center for sampling
+radius = np.pi  # radius of sampling
+sampling = "linspace"  # use a uniform grid for sampling
+n_samples = 10  # how many samples to use
 
-# Which solver to use: can be MOSEK, newton (simple log-barrier Newton 
+# Which solver to use: can be MOSEK, newton (simple log-barrier Newton
 # solver), newton-features or newton-kernel (more advanced solvers using
 # feature or kernel matrices, respectively; adding e.g. a linesearch option
-# and more advanced diagnostics to the original solver. 
-solver = "newton" 
+# and more advanced diagnostics to the original solver.
+solver = "newton"
 
-# Which kernel to use: use Gauss for smooth, Laplace for less smooth, 
-# or provide a kernel of your choice. 
-kernel = "Gauss" 
+# Which kernel to use: use Gauss for smooth, Laplace for less smooth,
+# or provide a kernel of your choice.
+kernel = "Gauss"
 
 # Kernel parameter: good to use roughly the minimum distance between
-# samples. 
-sigma = 2 * np.pi / n_samples 
+# samples.
+sigma = 2 * np.pi / n_samples
 
 solution, info = ksos.solve(
     f=np.sin,
@@ -29,7 +29,7 @@
     center=center,
     radius=radius,
     solver=solver,
-    sigma=sigma
+    sigma=sigma,
 )
 print(f"Found solution: x={solution[0]:.4f}, f={info['cost']:.4f}")
 print(f"True solution:  x={-np.pi/2:.4f}, f=-1")

ksos_env.yml

@@ -12,6 +12,7 @@ dependencies:
   - matplotlib
   - numpy
   - scipy
+  - pytest
   - pip:
     - mosek
     - -e .

ksos_tools/examples/polynomial.py

@@ -1,5 +1,3 @@
-import warnings
-
 import matplotlib.pylab as plt
 import numpy as np
 
@@ -36,8 +34,6 @@ def __repr__(self):
 
 
 if __name__ == "__main__":
-    import matplotlib.pylab as plt
-
     poly = Polynomial()
     poly.plot()
     plt.show(block=False)

ksos_tools/solvers/helpers.py

@@ -1,7 +1,6 @@
 import cvxpy as cp
-import numpy as np
 
-from ksos_tools.utils import duplication_matrix, hvec
+from ksos_tools.utils import hvec
 
 
 def find_feasible_B(alpha, c, problem, soft=False):

ksos_tools/solvers/ksos.py

@@ -20,6 +20,7 @@ def solve(
     radius: np.ndarray | float | None = None,
     n_samples: int | None = None,
     samples: np.ndarray | None = None,
+    f_samples: np.ndarray | None = None,
     sampling: str = "uniform",
     sampling_function: Callable[[], np.ndarray] | None = None,
     lambd: float = 1e-3,
@@ -101,7 +102,12 @@ def solve(
         radius = np.array([radius] * dim)
     assert n_samples is None or n_samples >= 1
     assert lambd >= 0
-    assert sigma > 0
+    assert (
+        kernel == "Periodic"
+        and isinstance(sigma, tuple)
+        and sigma[0] > 0
+        and sigma[1] > 0
+    ) or sigma > 0
     assert warm_iterations >= 1
     if samples is not None:
         assert np.ndim(samples) >= 2
@@ -158,7 +164,10 @@ def solve(
         # generate samples and kernel matrix
 
         if samples is not None:
-            problem.register_fixed_samples(samples, f)
+            if f_samples is not None:
+                problem.register_fixed_samples(samples, None, f_samples)
+            else:
+                problem.register_fixed_samples(samples, f, None)
 
             if solver != "naive":
                 success = problem.initialize_kernel(
@@ -279,7 +288,7 @@ def solve(
             radius = radius * d
             if verbose:
                 print(f"Sobolev norm: {norm}  |  Decay: {d}  |  New radius: {radius}")
-        else:
+        elif iteration < warm_iterations - 1:
             assert radius is not None
             assert sigma is not None
             assert isinstance(decay, float)
@@ -393,7 +402,7 @@ def get_surrogate(
 
         for i, (fi_interp, fi) in enumerate(zip(values_samples, f_samples_min_c)):
             if abs((fi_interp - fi) / fi) > 1e-2:
-                msg = f"Warning: at sample {i}, surrogate function not passing directly through f: {float(fi_interp):.4f}, {float(fi):.4f}"
+                msg = f"Warning: at sample {i}, surrogate function not passing directly through f: {float(fi_interp.item()):.4f}, {float(fi.item()):.4f}"
                 if errors == "print":
                     print(msg)
                 elif errors == "raise":

ksos_tools/solvers/newton.py

@@ -1,19 +1,15 @@
 import warnings
 
-import eigenpy
 import numpy as np
 import scipy
 import scipy.linalg
+from ksos_tools.solvers.problem import Problem
 
 DEBUG = False
 
 TOL_EIG_MIN = -1e-8  # values below this are considered negative.
 
 
-from ksos_tools.solvers.helpers import find_feasible_B
-from ksos_tools.solvers.problem import Problem
-
-
 def armijo_linesearch(cost_before, alpha, stepsize, delta, grad_H, problem: Problem):
     beta = 0.25
     gamma = 1e-4
@@ -175,7 +171,7 @@ def solve_newton_system(alpha):
             # solve = lambda x: llt.solve(x)
             a, b = scipy.linalg.cho_factor(H_pp)
             solve = lambda x: scipy.linalg.cho_solve((a, b), x)
-        except:
+        except scipy.linalg.LinAlgError:
             warnings.warn("Hessian is not positive definite?")
             a, b = scipy.linalg.lu_factor(H_pp)
             solve = lambda x: scipy.linalg.lu_solve((a, b), x)
@@ -302,7 +298,7 @@ def solve_newton_system(alpha, grad_H, hess_H):
             solve = lambda x: scipy.linalg.cho_solve((a, b), x)
             # llt = eigenpy.LLT(hess_H)
             # solve = lambda x: llt.solve(x)
-        except:
+        except scipy.linalg.LinAlgError:
             a, b = scipy.linalg.lu_factor(hess_H)
             solve = lambda x: scipy.linalg.lu_solve((a, b), x)
         g1 = solve(grad_H).flatten()

ksos_tools/solvers/problem.py

@@ -43,6 +43,9 @@ def kernel_function(x, y, sigma, kernel):
         return np.exp(-np.linalg.norm(x - y) ** 2 / (2 * sigma**2))
     elif kernel == "Polynomial":
         return (1 + np.inner(x, y)) ** int(sigma)
+    elif kernel == "Periodic":
+        p, l = sigma
+        return np.prod(np.exp(-2 * (np.sin(np.pi * (x - y) / p) ** 2) / l**2))
     else:
         raise ValueError(f"unknown kernel function {kernel}")
 
@@ -63,7 +66,7 @@ def decompose(K, method):
     elif method == "numpy":
         try:
             R = np.linalg.cholesky(K).T
-        except:
+        except np.linalg.LinAlgError:
             raise ValueError("Kernel matrix not positive!")
     elif method == "eigenpy":
         llt = eigenpy.LLT(K)
@@ -254,7 +257,7 @@ def generate_new_samples(
         ):
             radius = np.array([radius] * dim)
         assert n_samples >= 1
-        assert sampling in ["linspace", "uniform"]
+        assert sampling in ["linspace", "uniform", "sobol"]
 
         # Generate samples
         if sampling_function is not None:

ksos_tools/solvers/sos.py

@@ -100,7 +100,7 @@ def solve_from_samples(
             return None, {"cost": None, "ttot": ttot}
         else:
             ttot = time.time() - t1
-            if not "optimal" in prob.status:
+            if "optimal" not in prob.status:
                 print("No solution found:", prob.status)
                 return None, {"cost": None, "ttot": ttot}