What is Bézier simplex fitting?
You are probably familiar with Bézier curves (1-D) and Bézier triangles (2-D) from computer graphics and CAD software. A Bézier simplex is their natural generalization to any number of dimensions: the same elegant polynomial construction, extended to a hypersurface of arbitrary dimension defined over a standard simplex.
At its core, Bézier simplex fitting is a general-purpose regression technique. Just as a regular Bézier curve smoothly interpolates or approximates a 1-D point cloud using a small set of control points, a Bézier simplex can approximate any continuous map from a standard simplex to a high-dimensional Euclidean space. Given a point cloud dataset defined over simplex coordinates, it can fit a highly flexible and mathematically well-behaved parametric surface to the data.
This page introduces the formal definition of Bézier simplices, the least-squares fitting algorithm used by PyTorch-BSF, and its most prominent real-world applications. When the data are better described on a product of simplices (a prism, a square, a cube) or on several such cells glued together, PyTorch-BSF fits Bézier simploids instead (see Beyond Simplices: Bézier Simploids).
Bézier Simplex
Let \(D, M, N\) be non-negative integers, \(\mathbb N\) the set of non-negative integers (including zero), and \(\mathbb R^N\) the \(N\)-dimensional Euclidean space. We define the index set by
and the simplex by
An \((M-1)\)-dimensional Bézier simplex of degree \(D\) in \(\mathbb R^N\) is a polynomial map \(\mathbf b: \Delta^{M-1}\to\mathbb R^N\) defined by
where \(\mathbf t^{\mathbf d} = t_1^{d_1} t_2^{d_2}\cdots t_M^{d_M}\), \(\binom{D}{\mathbf d}=D! / (d_1!d_2!\cdots d_M!)\), and \(\mathbf p_{\mathbf d}\in\mathbb R^N\ (\mathbf d\in\mathbb N_D^M)\) are parameters called the control points.
Fig. 2 A 2-D Bézier simplex of degree 3 in \(\mathbb R^3\) and its control points. The shape of the simplex is determined by the control points \(\mathbf p_{(3,0,0)}, \mathbf p_{(2,1,0)}, \ldots, \mathbf p_{(0,0,3)}\).
Bézier Simplex Fitting
Assume we have a finite dataset \(B\subset\Delta^{M-1}\times\mathbb R^N\) and want to fit a Bézier simplex to the dataset. What we are trying can be formulated as a problem of finding the best vector of control points \(\mathbf p=(\mathbf p_{\mathbf d})_{\mathbf d\in\mathbb N_D^M}\) that minimizes the least squares error between the Bézier simplex and the dataset:
For fixed parameter locations, the predictions are linear in the control
points. PyTorch-BSF solves this linear least-squares problem with a dense
solver by default, with an optional sparse iterative solver for larger
systems. The solve uses CPU float64. This is separate from using
minimize to optimize a differentiable application objective; see
Objective minimization.
Fig. 3 A Bézier simplex fitted to a dataset. The control points are determined by the least squares fitting algorithm.
Approximation Theorem
Any continuous map from a simplex to a Euclidean space can be approximated by a Bézier simplex. More precisely, the following theorem holds [KHS+19]:
Theorem 1 (Universal Approximation Theorem)
For any continuous map \(\phi: \Delta^{M-1} \to \mathbb{R}^N\) and any \(\epsilon > 0\), there exists a degree \(D\) and control points \(\mathbf{p}\) such that the Bézier simplex \(\mathbf{b}(\mathbf{t} \mid \mathbf{p})\) satisfies \(\max_{\mathbf{t} \in \Delta^{M-1}} \| \phi(\mathbf{t}) - \mathbf{b}(\mathbf{t} \mid \mathbf{p}) \| < \epsilon\).
This guarantees that Bézier simplices are universal approximators for any continuous simplex-domain function.
Beyond Simplices: Bézier Simploids
A single simplex is not always the right parameter domain. A simploid of type \((a_1, \ldots, a_n)\) is the product \(\Delta^{a_1} \times \cdots \times \Delta^{a_n}\) of simplices, e.g., a prism \(\Delta^2 \times \Delta^1\) or a square \(\Delta^1 \times \Delta^1\), and a simplex is the simploid of type \((M-1)\). A Bézier simploid is the tensor product of Bézier simplices, one per factor, and is fitted by the same least-squares formulation as above with \(\Delta^{M-1}\) replaced by the simploid. Since the fitting problem is linear in the control points, torch_bsf.fit() solves it exactly by default, for Bézier simplices and simploids alike.
The approximation theorem carries over: polynomials are dense in the continuous maps on a compact product of simplices (Stone–Weierstrass), and every polynomial of degree at most \(d_i\) in the \(i\)-th factor is a Bézier simploid of degree \((d_1, \ldots, d_n)\). Simploidal sets go further and glue several simploids along common faces. See Single and assembled domains for details.
Build a geometric picture before choosing a domain
Start with a segment \(\Delta^1\). Choosing one point on it takes one number, such as \(u\in[0,1]\). Choosing points on two independent segments takes two numbers \((u,v)\): all their combinations fill a square. Likewise, choosing a point in a triangle and an independent height fills a triangular prism. This is what the word product means here. The type lists the dimensions of the factors, so a square has type \((1,1)\) and a prism has type \((2,1)\). These are dimensions, not polynomial degrees.
Fig. 4 Left: the image of a square-domain Bézier simploid. Center: the boundary images of a prism-domain Bézier simploid. Right: two different domain cells, a triangle and a square, fitted with a shared edge. Red points and dashed lines show control nets; the colored patches show model images.
Read the panels as three distinct constructions:
Left, two independent intervals. The domain is a flat square, but its image in three output coordinates is a curved saddle. Holding one interval coordinate fixed gives a Bézier curve as the other varies. Its biquadratic degree \((2,2)\) gives \(3\times3=9\) controls. Interior controls generally lie away from the fitted surface.
Center, a triangle and a height. At a fixed height, the domain slice is a triangle. Sweeping height creates a three-dimensional prism. Its boundary has two triangular and three square faces; these are what the renderer displays. The twisted output image does not change the domain type \((2,1)\). Degree \((2,2)\) gives \(6\times3=18\) controls.
Right, gluing rather than taking a product. A triangle and a square are placed side by side and their selected edges are identified. Each keeps its own local coordinates. Shared edge controls make predictions agree along the whole seam (C0 continuity). They do not require matching tangents across it. This is an assembled simploidal domain rather than one larger product cell.
The distinction matters for data. Three competing priorities use a triangle: their three nonnegative weights share one sum-to-one constraint. Two priorities and an independent operating condition use a square: the trade-off and the condition can vary independently. Three priorities and an operating condition use a prism. Joining neighboring experimental regions calls for gluing charts if the response must agree at their common boundary.
In the API these choices are BezierSpace.simplex(2, degree=...),
BezierSpace.simploid((1, 1), degree=...) and
BezierSpace.simploid((2, 1), degree=...). For the right panel, build a
SimploidalSet and compile it with BezierSpace.on_set. Each product
factor has its own barycentric block: a square row is
[1-u, u, 1-v, v] and sums to two, while a triangle row sums to one.
Reproduce the three panels with
python examples/generate_simploidal_plot.py. Continue to
Single and assembled domains for runnable fitting and gluing examples, or
Multi-Objective Model Predictive Control for an operating-condition
parameter that makes a product domain useful.
Relation to Multi-Objective Optimization
Data suitable for modeling with a Bézier simplex are point clouds distributed along a low-dimensional (e.g., 1 to 10 dimensions) curved simplex lying within a high-dimensional ambient space (e.g., tens to thousands of dimensions). Such data can be regarded as samples from the solution set of a specific class of multi-objective optimization problems.
Multi-Objective Optimization
In many real-world applications, from engineering design to machine learning, we often need to optimize multiple conflicting criteria simultaneously (e.g., maximizing performance while minimizing cost). Because these objectives naturally compete with one another, it is generally impossible to find a single perfect solution. Instead, multi-objective optimization seeks to find a set of optimal trade-offs, providing the mathematical foundation to explore and select the best compromise.
Definition 1 (Multi-Objective Optimization Problem)
Let \(X\) be a feasible decision space. A multi-objective optimization problem with \(M\) objectives aims to minimize a vector-valued function:
Since the objectives typically conflict with one another, there is rarely a single solution that minimizes all \(M\) objectives simultaneously. Instead, optimality is defined in terms of trade-offs.
Definition 2 (Pareto Set and Pareto Front)
A solution \(x \in X\) is said to dominate another solution \(x' \in X\) if \(f_m(x) \le f_m(x')\) for all \(m \in \{1, \ldots, M\}\) and \(f_j(x) < f_j(x')\) for at least one index \(j\).
A solution \(x^* \in X\) is Pareto optimal if no other solution within \(X\) dominates it.
The Pareto set is the set of all Pareto optimal solutions in the decision space \(X\).
The Pareto front is the image of the Pareto set in the objective space \(f(X) \subset \mathbb{R}^M\).
Weakly Simplicial Problems
In many real-world multi-objective optimization problems, the Pareto set and Pareto front exhibit a highly structured, continuous shape. Specifically, they often mirror the topological structure of a standard simplex (e.g., a curve for a two-objective problem, a curved triangle for a three-objective problem, and so on). The concept of a weakly simplicial problem formally captures this property, ensuring that the optimal trade-off surfaces are well-behaved and can be elegantly approximated by Bézier simplices.
Definition 3 (Weakly Simplicial Problem [MHI21])
Let \(f:X\to\mathbb R^M\) be a mapping, where \(X\) is a subset of \(\mathbb R^N\). The problem of minimizing \(f\) is \(C^r\)-simplicial if there exists a \(C^r\)-mapping \(\Phi:\Delta^{M-1} \to X^\star(f)\) such that both the mappings \(\Phi|_{\Delta_I}:\Delta_I\to X^\star(f_I)\) and \(f|_{X^\star(f_I)}:X^\star(f_I)\to f(X^\star(f_I))\) are \(C^r\)-diffeomorphisms for any nonempty subset \(I\) of \(\{1,\ldots,M\}\), where \(0\le r\le\infty\). The problem of minimizing \(f\) is \(C^r\)-weakly simplicial if there exists a \(C^r\)-mapping \(\phi:\Delta^{M-1}\to X^\star(f)\) such that \(\phi(\Delta_I)=X^\star(f_I)\) for any nonempty subset \(I\) of \(\{1,\ldots,M\}\), where \(0\le r\le\infty\).
Fig. 5 A simplicial problem: the Pareto set and Pareto front are homeomorphic to a simplex, i.e., they have no pinched topology.
Fig. 6 A weakly simplicial problem: the Pareto set and Pareto front are a continuous image of a simplex, i.e., they may have a pinched topology.
In weakly simplicial problems, there exists a continuous map from a simplex to the Pareto set and Pareto front, which is guaranteed to be approximable by the Universal Approximation Theorem.
Strongly Convex Problems
While the concept of weakly simplicial problems provides a powerful framework, verifying this property for an arbitrary problem can be challenging. Fortunately, a broad and highly practical class of optimization problems—strongly convex problems—are mathematically guaranteed to be weakly simplicial. This means that if an optimization problem is strongly convex, its Pareto front inherently possesses a simplex-like structure, making it an ideal candidate for Bézier simplex fitting. As outlined by Hamada et al. [HHI+20], this simplicial topology of the Pareto set is a fundamental, rigorous mathematical property of strongly convex multi-objective problems.
Strong convexity commonly arises in practice through four natural structural mechanisms: (1) L2 regularization, such as ridge penalties; (2) Positive-definite quadratic forms, like risk covariances or LQR stage costs; (3) Squared Euclidean distances, prevalent in facility location; and (4) Composite structures, where a strongly convex smooth term is paired with a convex nonsmooth penalty.
Definition 4 (Strongly Convex Problem)
A multi-objective optimization problem is strongly convex if its feasible decision space \(X\) is convex, and every objective function \(f_m\) (\(m=1,\ldots,M\)) is strongly convex. Formally, a function \(f_m\) is strongly convex with parameter \(\mu > 0\) if for all \(x, y \in X\) and \(t \in [0, 1]\):
The formal mathematical foundation for this connection is given by the following theorems, which prove that strongly convex problems are inherently weakly simplicial and that their optimal solutions can be continuously parameterized by a standard simplex.
Theorem 2 (Theorems 1 and 2 in [MHI21])
Let \(f: \mathbb{R}^n \to \mathbb{R}^m\) be a \(C^r\)-strongly convex mapping (\(0 \le r \le \infty\)). Then, the problem of minimizing \(f\) is \(C^{r-1}\)-weakly simplicial for \(r > 0\) and \(C^0\)-weakly simplicial for \(r = 0\).
Theorem 3 (Proposition 1 of [MHI21])
Let \(f: \mathbb{R}^n \to \mathbb{R}^m\) be a strongly convex mapping. Then, the mapping \(x^*: \Delta^{m-1} \to X^*(f)\) defined by
is surjective and continuous.
This guarantees that for strongly convex models, their Pareto fronts admit a simplex structure and can be efficiently reconstructed using Bézier simplex fitting.
Statistical Test for Simplicial Topology
When the underlying problem class is not analytically known in advance (e.g., in complex real-world designs using simulations), it is not immediately clear whether the Pareto set admits a simplex structure. In such cases, data-driven statistical tests based on persistent homology can determine whether this assumption is warranted before committing to a Bézier simplex model.
Earlier work introduced the formal concept of a simple problem:
Definition 5 (Simple Problem [HG18])
A problem is simple if it satisfies two conditions:
Theorem 4 (Equivalence of Simple and Simplicial Problems)
All \(C^r\)-simplicial problems are simple (\(0\le r \le \infty\)). A simple problem is \(C^0\)-simplicial if the interiors of Pareto sets for any two subproblems do not intersect.
Furthermore, it is conjectured that “simple” and “simplicial” are equivalent even without this extra condition. Therefore, testing whether a problem is simple serves as a practical test for whether a problem is simplicial. Notably, the non-intersection condition itself can also be tested statistically.
To strictly determine whether a problem is simple, two complementary statistical tests are required. In standard statistical testing, failing to reject a null hypothesis does not allow one to affirmatively adopt it. Therefore, we need both a test to reject simplicity and a test to affirmatively confirm it.
1. Testing that a problem is not simple
Hamada and Goto [HG18] introduced a data-driven method using persistent homology and its confidence sets to test for violations of the simplicity conditions. Specifically, we can reject the simplicity hypothesis if the estimated Betti numbers do not match those of a topological simplex:
Theorem 5 (Theorem 3.1 in [HG18] : Test for (S1) violation)
Let \(K(\mathcal{F})\) be a simplicial complex representing the approximated Pareto set. Assuming its geometric realization is homotopy equivalent to the true Pareto set, if one of the following homological conditions is satisfied:
\(H_0(K(\mathcal{F})) \not\cong \mathbb{Z}\)
\(\exists i > 0 : H_i(K(\mathcal{F})) \not\cong 0\)
then the problem does not satisfy the simplicity condition (S1).
2. Testing that a problem is simple
Building upon this, Hamada [Ham20] demonstrated that this homology-based approach can also be used to affirmatively test whether a problem is simple. By leveraging the h-cobordism theorem from algebraic topology, they established a mathematical characterization showing that if the Pareto set (modeled as a compact smooth manifold) and its boundary are simply connected, and its estimated homology matches that of a standard simplex, then the Pareto set is guaranteed to be homeomorphic to a simplex. This affirmative test complements the non-simpliciality test, allowing practitioners to strictly verify the simpliciality of a multi-objective optimization problem from sampled Pareto solutions before applying Bézier simplex fitting.
Theorem 6 (Theorem 3 in [Ham20])
Let \(M\) be an \(n\)-dimensional compact \(C^\infty\) manifold. If both \(M\) and its boundary \(\partial M\) are simply connected, and its homology groups satisfy:
then \(M\) is homeomorphic to the standard simplex \(\Delta^n\).