Fitting from values alone

Use fit_point_cloud when each observation is a point on a curve or surface and its domain coordinate is unknown. The existing fit solves a different problem: regression from known parameter/value pairs.

The value-only methods implement Algorithms 2 and 3 of Kobayashi et al. [KHS+19]. Both currently support a single simplex. fit_point_cloud returns a PointCloudFitResult; fit_skeleton returns a SkeletonFitResult. Each result holds an ordinary BezierModel in result.model for prediction, plotting and serialization. Inspect result.success, result.status and result.message for numerical termination, or call result.require_success() to raise on failure.

All-at-once fitting

The unknowns are both the control points and a barycentric coordinate for every observation. fit_point_cloud alternates two steps:

Hold the model fixed and minimize each point’s squared distance to the surface using analytic first and second derivatives and Newton steps.
Hold the inferred coordinates fixed and solve for the controls using the existing float64 SVD least-squares implementation.

For example, fit a two-objective Schaffer front without providing its parameter coordinates to the fitting routine. The three examples on this page share the named point-cloud-demo session: run them in page order; the later examples reuse torch, space and values from this one.

import torch
from torch_bsf import BezierSpace, fit_point_cloud

u = torch.linspace(0, 1, 21, dtype=torch.float64)
values = torch.stack((u**2, (1 - u)**2), dim=1)
space = BezierSpace.simplex(dimension=1, degree=2)
result = fit_point_cloud(space, values).require_success()

assert result.success and result.status == "converged"
assert result.parameters.shape == (len(values), 2)
assert torch.all(result.parameters >= 0)
torch.testing.assert_close(
    result.parameters.sum(1), torch.ones(len(values), dtype=torch.float64),
    rtol=0, atol=8 * torch.finfo(torch.float64).eps,
)
fitted_values = result.model(result.parameters)
# Row i reconstructs values[i], rather than a reordered cloud.
row_errors = (fitted_values - values).square().sum(1)
assert row_errors.max() < 1e-12

result.parameters contains the inferred barycentric rows in the original observation order, with shape (n_observations, dimension + 1) for a positive-dimensional simplex. Here the 21 observations have two barycentric coordinates each; result.model(result.parameters)[i] approximates values[i]. A zero-dimensional simplex uses empty coordinate rows, consistent with BezierSpace. objective_history records mean squared Euclidean distances, starting with the initial model and coordinates; control_result retains the last linear solve’s diagnostics. Controls supplied through fixed=FixedControls(...) remain exact throughout fitting.

Automatic initialization places vertex i at the sample minimizing output i, then initializes the remaining controls on the affine grid between vertices. This follows the paper and assumes minimization objectives, with one output per vertex. Use initial=BezierModel(...) to specify another orientation or fit a simplex embedded in a different output dimension. The initial model is copied.

Inductive skeleton fitting

fit_skeleton fits vertices first, then edges, then higher faces. Each face keeps every previously estimated boundary control fixed exactly. Only faces with at most min(degree, number_of_vertices) vertices are needed: a degree-three simplex has no controls supported on four or more vertices, even when its domain dimension is larger.

from torch_bsf import fit_skeleton

skeleton = fit_skeleton(space, values).require_success()
model = skeleton.model
assert skeleton.success and skeleton.status == "converged"
assert set(skeleton.faces) == {(0,), (1,), (0, 1)}
edge_result = skeleton.faces[(0, 1)]
assert edge_result.success
assert edge_result.parameters.shape == (len(values), 2)
# The assembled edge retains the fitted vertex controls exactly.
for vertex, powers in [(0, (2, 0)), (1, (0, 2))]:
    vertex_result = skeleton.faces[(vertex,)]
    assert vertex_result.success
    control = model.coefficients[space.control_index(powers)]
    assert torch.equal(control, vertex_result.model.coefficients[0])

skeleton.faces maps sorted vertex-ID tuples to PointCloudFitResult objects. In this segment example, (0,) and (1,) identify the vertices and (0, 1) identifies the whole edge. Each face has its own numerical status and locally inferred coordinates; the assembled predictor is skeleton.model.

By default, a face’s samples are the points nondominated under the objectives corresponding to its vertices, using the paper’s definition of X_J. All output columns are retained. This extraction assumes that output columns are minimization objectives and that the Pareto front has the stated simplex structure. It takes quadratic time in sample count per face and linear auxiliary memory. For larger data or observations already separated by subproblem, supply the samples explicitly:

faces = {
    (0,): values[:1],
    (1,): values[-1:],
    (0, 1): values,
}
explicit = fit_skeleton(space, values, face_values=faces).require_success()
assert all(face.success for face in explicit.faces.values())
for vertex, powers in [(0, (2, 0)), (1, (0, 2))]:
    control = explicit.model.coefficients[space.control_index(powers)]
    assert torch.equal(control, faces[(vertex,)][0])
edge = explicit.faces[(0, 1)]
assert edge.parameters.shape == (len(values), 2)
row_errors = (edge.model(edge.parameters) - faces[(0, 1)]).square().sum(1)
assert row_errors.max() < 1e-12

Face keys are sorted tuples of zero-based vertex IDs. Supply exactly all required faces and full output vectors for each sample. No parameter coordinates are needed. Explicit face_values and an initial model also allow fitting arbitrary embeddings. A vertex with repeated or noisy observations is fitted to their mean. Each face result contains coordinates local to that face, in its supplied or extracted sample order. A failed face stops assembly and leaves an explicitly unsuccessful partial result.

Numerical interpretation

The parameter updates enforce nonnegative coordinates summing to one. Projection, backtracking and a gradient fallback safeguard Newton updates at boundaries and singular or indefinite Hessians. These safeguards extend the paper’s plain Newton description. The algorithm is a local nonlinear method: different initializations can give different surfaces and domain parameterizations are generally not identifiable from values alone.

tolerance bounds the absolute change in root mean squared distance; parameter_tolerance bounds the projected gradient residual. The code uses an absolute change rather than the signed stopping expression printed in the paper. Success requires both stopping checks. Iteration exhaustion returns success=False, status="iteration_limit" and emits ConvergenceWarning, preserving the finite model. Use require_success() when a failed fit must halt a workflow. Convergence is a numerical stopping statement, not a certificate of global optimality or reconstruction accuracy.

Distances and tolerances use the original output units. Normalize outputs explicitly when their scales should contribute equally. These routines do not automatically filter the input cloud to a Pareto front or infer face memberships for arbitrary embeddings. Choose degree and validate surface coverage separately from numerical convergence.