Unsupervised Segmentation and Alignment of Multi‑Demonstration Trajectories via Multi‑Feature Saliency and Duration‑Explicit HSMMs

1. Introduction

Learning from Demonstration (LfD) [1] seeks to transfer skills from a handful of expert executions to new task instances. In real deployments across aerial, driving, and manipulation domains, multiple demonstrations routinely exhibit irregular local time warping (hovering, waiting, backtracking) [2,3,4], sensor noise on high-dimensional signals [5,6], and alternating quasi-stationary and transition phases [7,8,9]. Our objective is to discover phase structure without labels, align multiple demonstrations onto a shared semantic time base, and provide calibrated, segment-wise probabilistic models that plug into constrained planners.

A central difficulty is that widely used segmentation and alignment tools either operate at the signal-threshold level or impose duration assumptions that do not reflect how humans actually perform tasks [10,11]. Threshold/template methods or simple peak detectors are brittle to sampling-rate and environment shifts; their decision surfaces move with scale changes, causing under- and over-segmentation in the presence of jitter and speed variations [12,13]. Bayesian online changepoint detection (BOCPD) relieves fixed thresholds but still hinges on hazard-rate and noise assumptions; in practice it often fragments long dwells into multiple short segments when features oscillate, and it provides no explicit duration model for downstream use [14,15,16]. Dynamic Time Warping (DTW) offers pairwise alignment but lacks a generative mechanism and cannot represent dwell distributions or uncertainty, which limits its utility for planning and simulation [17,18,19].

Latent-state models address some of these issues, yet common instantiations come with their own bias-variance trade-offs. Hidden Markov Models (HMMs) assume geometric state durations; this geometric-duration bias shortens non-geometric dwells and shifts boundaries when operators hesitate or loiter, leading to misaligned phases across trials [20]. Hidden semi-Markov Models (HSMMs) are a principled remedy because they explicitly model durations, but most prior pipelines still depend on hand-crafted feature stacks with fixed weights/hyper-parameters, trained per-demonstration or per-session [1,21,22]; robustness degrades under operator/platform shifts and tempo changes. Bayesian nonparametrics can adapt model complexity [23], but the computational overhead and stability concerns have limited their use in long-horizon, multi-demo settings.

Downstream generation has a parallel set of limitations. Dynamical Movement Primitives (DMPs) provide smooth, time-scalable execution but lack closed-form uncertainty; Gaussian Mixture Regression (GMR) and Probabilistic Movement Primitives (ProMPs) offer distributional predictions with covariances, yet both are sensitive to boundary errors—mis-segmentation inflates variances and biases means within segments [20]. More importantly, most pipelines treat segmentation/alignment and probabilistic encoding as loosely coupled stages; there is rarely a mechanism to self-calibrate features using multi-demo consistency or to learn an alignment that is simultaneously scale-aware and duration-aware [24]. Finally, while planners such as OMPL/MPC can exploit covariance for risk-aware control, few works deliver calibrated uncertainty on a shared semantic time base that is directly consumable by these planners [25].

To address these gaps, we propose a label-free pipeline that couples unsupervised segmentation, duration-explicit alignment, and probabilistic encoding into a single, self-calibrating loop. First, we compute a dimensionless multi-feature saliency by fusing velocity, acceleration, curvature, and direction-change rate; then we apply topology-aware keyframe extraction using persistent peak–valley pairs to keep only structurally significant extrema and non-maximum suppression to avoid clustered responses [26]. Second, we jointly train an HSMM across demonstrations with explicit duration distributions and an extended forward–backward recursion, producing a shared semantic time axis and phase-consistent boundaries [25]. Third, within each phase we fit probabilistic motion models—GMM/GMR or ProMP, optionally combined with DMP for execution—to obtain mean trajectories with calibrated covariances [20]. Crucially, we close the loop by learning the saliency weights without labels: a CMA-ES [27] search on the probability simplex minimizes cross-demonstration structural dispersion, automatically rebalancing features so that segmentation and alignment are mutually consistent. Compared with BOCPD, thresholding, or HMM-based pipelines, this design is explicitly duration-aware, scale/time-warp robust, and planner-ready with uncertainty that is calibrated and comparable across operators.

Contributions.

• Scale-/time-warp-robust saliency. A topology-aware, multi-feature saliency (persistence and non-maximum suppression) that stabilizes keyframes under noise and tempo variations, yielding sparser yet more stable anchors than signal-level detectors.
• Joint HSMM alignment with explicit durations. Cross-demo training with extended forward–backward/Viterbi recursions; model order selected by a joint criterion combining BIC and alignment error (AE) to balance fit and parsimony.
• Label-free feature-weight self-calibration. CMA-ES on the weight simplex to minimize cross-demo structural dispersion, eliminating hand-tuned fusion and improving phase consistency.
• Calibrated probabilistic encoding for planning. Segment-wise GMR/ProMP (optionally fused with DMP) returning means and covariances that integrate directly with OMPL/MPC for risk-aware execution.

Empirically, on UAV flight, autonomous-driving, and manipulation benchmarks, our method reduces phase-boundary dispersion by ≈31% on UAV-Sim and by 30–36% under monotone time warps, additive noise, and missing data compared with HMM variants; it improves the sparsity–fidelity–smoothness trade-off (higher time compression at comparable reconstruction error with lower jerk) and achieves nominal 2σ coverage (94–96%), indicating well-calibrated uncertainty. Section 2 details the pipeline; Section 3 reports datasets, metrics, baselines (including BOCPD and HMM), ablations, and robustness; Section 4 discusses positioning, limitations, and future work; Section 5 concludes.

2. Materials and Methods

Figure 1 sketches the end-to-end workflow that maps multiple demonstrations to executable trajectories under physical and safety constraints: (i) unsupervised segmentation from multi-feature saliency with topology-aware keyframes; (ii) duration-explicit alignment across demonstrations via a hidden semi-Markov model (HSMM) with shared parameters; and (iii) probabilistic encoding of each phase (GMM/GMR or ProMP, optionally combined with DMP) producing mean trajectories with calibrated covariances. Outputs include phase-consistent labels, segment-wise probabilistic models, and constraint-aware executable trajectories.

2.1. Inputs, Outputs, and Assumptions

• Inputs.

We observe M independent demonstrations (optionally multi-modal), ${\left\{P^{(m)}\right\}}_{m=1}^M,P^{(m)}={\left\{P^{(m)}(t)\right\}}_{t=1}^{T_m},$sampled with fixed period $\mathsf{\Delta }t$. When available, auxiliary channels (e.g., pose, force/tactile, depth) are concatenated into the observation vector used downstream.

• Outputs.

A shared set of semantic phases ${\left\{S_k\right\}}_{k=1}^{N^{\star }}$ and per-demo boundaries $\left\{{\tau }_k^{(m)}\right\}$;

For each phase, a segment-wise generative model—DMP, GMM/GMR, or ProMP, alone or in combination—returning mean trajectories and covariance estimates;

Under constraints $\mathcal{C}$ (terminal/waypoints, velocity/acceleration limits, etc.), an executable trajectory and associated risk measures computed from covariances.

• Assumptions.

• Demonstrations consist of alternating quasi-stationary and transition segments.
• Time deformation is order-preserving (the semantic phase order does not change).
• Observation noise is moderate and can be mitigated by local smoothing and statistical filtering.

2.2. Multi-Feature Analysis and Automatic Segmentation

2.2.1. Data Ingestion and Pre-modulation

Let a single demonstration be the discrete sequence $\mathcal{D}={\left\{\mathsf{\Delta }{\mathbf{x}}_t\right\}}_{t=1}^T,\mathsf{\Delta }{\mathbf{x}}_t\in {\mathbb{R}}^d$, where $d$ is the number of observed degrees of freedom (e.g., $d=3$ for UAV or mobile bases; $d=6$ for industrial manipulators). Cumulative sums yield the Cartesian trajectory of the tool center point

$$\mathbf{P}(t)={\left[p_x(t),p_y(t),p_z(t)\right]}^{\mathrm{\top }},t=1,\dots ,T,$$

(2.1)

We apply a one-dimensional Savitzky–Golay local polynomial smoother to $\mathbf{P}(t)$ before computing derivatives, suppressing high-frequency tele-operation jitter and stabilizing numerical differentiation. The filter is equivalent to an $(r+1)$-order Taylor truncation in time and ensures

$$\Vert \tilde{\mathbf{P}}-\mathbf{P}{\Vert }_{\mathrm{\infty }}\le c_W\epsilon ,c_W<1，$$

where $\epsilon$ bounds the measurement noise and the window $W\in \{\mathrm{5,7},9\}$. Unless stated otherwise, $\tilde{\mathbf{P}}$ denotes the smoothed trajectory [28].

Implementation note. Central differences are used for derivatives; all channels are time-synchronized at a fixed sampling period $\mathsf{\Delta }t$.

2.2.2. Feature Computation and Saliency Fusion

On the smoothed trajectory we compute four complementary, time-varying features that reveal when, where, and how kinematic changes occur ($\mathsf{\Delta }t\left(={\displaystyle \frac{1}{f_s}}\right)$ is fixed).

• Velocity. Let

$$\mathbf{v}(t)={\displaystyle \frac{\mathbf{P}(t+1)-\mathbf{P}(t)}{\mathsf{\Delta }t}},v(t)=\Vert \mathbf{v}(t){\Vert }_2.$$

(2.2)

2.

Acceleration.

$$\mathbf{a}(t)={\displaystyle \frac{\mathbf{v}(t+1)-\mathbf{v}(t)}{\mathsf{\Delta }t}},a(t)=\Vert \mathbf{a}(t){\Vert }_2.$$

(2.3)

Peaks indicate abrupt speed changes.

3.

Curvature. With ΔP(t)=P(t+1)-P(t),

$$\kappa (t)={\displaystyle \frac{2\Vert \mathsf{\Delta }\mathbf{P}(t-1)\times \mathsf{\Delta }\mathbf{P}(t)\Vert }{\Vert \mathsf{\Delta }\mathbf{P}(t-1)\Vert \Vert \mathsf{\Delta }\mathbf{P}(t)\Vert \Vert \mathsf{\Delta }\mathbf{P}(t)-\mathsf{\Delta }\mathbf{P}(t-1)\Vert }}$$

(2.4)

Curvature measures spatial bending and is naturally invariant to global time scaling due to the cubic velocity term [29,30].

4.

Direction-Change Rate (DCR). Define the unit-direction vector $\widehat{\mathbf{v}}(t)={\displaystyle \frac{\mathbf{v}(t)}{\Vert \mathbf{v}(t){\Vert }_2}}$. To avoid numerical issues at very low speed, introduce a threshold $v_{min}>0$ and set

$$\mathrm{D}\mathrm{C}\mathrm{R}(t)=\left\{\begin{array}{cc}\Vert \widehat{v}(t)-\widehat{v}(t-1){\Vert }_2,\hfill & \hfill \Vert v(t){\Vert }_2\ge v_{\min }\\ 0,\hfill & \hfill \Vert v(t){\Vert }_2<v_{\min }\end{array}\right.$$

(2.5)

5.

Dimensionless fusion. Apply min–max normalization to each feature to obtain $\tilde{v}(t),\tilde{a}(t),\tilde{\kappa }(t),\overline{\mathrm{D}\mathrm{C}\mathrm{R}}(t)\in \left[\mathrm{0,1}\right]$. For a weight vector $w={\left[w_v,w_a,w_{\kappa },w_d\right]}^{\mathrm{\top }}$ with $w_i\ge 0$ and ${\sum }_i w_i=1$, define the fused saliency

$$\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}(t;w)=w_v\tilde{v}(t)+w_a\tilde{a}(t)+w_{\kappa }\tilde{\kappa }(t)+w_d\widetilde{\mathrm{D}\mathrm{C}\mathrm{R}}(t)$$

(2.6)

The weights $w$ are learned without labels in Section 2.2.4.

2.2.3. Keyframe Extraction with Topological Simplification

The saliency $\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}(t)$ compresses multi-source information into a 1-D signal, but segmentation should rely on structural extrema (global landmarks), not every minor fluctuation. We adopt a bottom-up screening that contracts a dense set of local extrema into a sparse, stable keyframe set [1,31].

• Candidate extrema via quantile thresholds.

Let $\mathcal{Q}=\left\{q_1,\dots ,q_L\right\}\subset \left(\mathrm{0,1}\right)$be a grid of quantiles (e.g., uniformly in [0.60,0.95]). For each $q\in \mathcal{Q}$:

• set ${\tau }_q={\mathrm{quantile}}_q(\mathrm{Score})$;
• collect indices ${\tilde{\mathrm{E}}}_q=\left\{t\mid \mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}(t)>{\tau }_q\right\}$ and snap each $t$ to the nearest local extremum within a radius-3 neighborhood.

To pick a unique $q^{\star }$, minimize a sparsity–fidelity loss

$$L(q)=\underset{\mathrm{sparsity}}{\underbrace{\left|{\tilde{E}}_q\right|/T}}+\lambda \underset{\mathrm{reconstruction}}{\underbrace{\mathrm{M}\mathrm{S}\mathrm{E}\left(P,{\widehat{P}}_{{\tilde{E}}_q}\right)}},$$

where ${\widehat{P}}_{{\tilde{\mathrm{E}}}_q}$ is the spline reconstruction at ${\tilde{\mathrm{E}}}_q$ and $\lambda >0$ reflects the admissible reconstruction error. Set $q^{\star }=\mathrm{a}\mathrm{r}\mathrm{g}\underset{q\in \mathcal{Q}}{min} \mathcal{L}(q),\tilde{\mathrm{E}}={\tilde{\mathrm{E}}}_{q^{\star }}$.

2.

Persistence thresholding (scale-invariant importance).

For adjacent peak–valley pairs $\left(p_{max},p_{min}\right)$ of $\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}(t)$, define the persistence

$$\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{s}\left(p_{max},p_{min}\right)=\left|\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}\left(p_{max}\right)-\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}\left(p_{min}\right)\right|$$

(2.7)

Small persistence typically indicates noise or micro-tremor; large persistence corresponds to genuine kinematic transitions. Plot

$$g(\alpha )=\mid \left\{\left(p_{max},p_{min}\right)\mid \mathrm{pers}>\alpha \right\}\mid ,$$

which empirically exhibits a plateau–cliff–stable pattern; the elbow ${\alpha }^{\star }$ is detected by Kneedle. Keep

$${\mathrm{E}}^{†}=\left\{p_i\mid \mathrm{p}\mathrm{e}\mathrm{r}\mathrm{s}\left(p_i\right)>\right\}$$

Because persistence depends only on amplitude differences, it is invariant to vertical scaling and mild time stretching, facilitating cross-demo comparability [32,33].

3.

Non-maximum suppression (NMS).

To avoid peak clustering, scan ${\mathrm{E}}^{†}$ with a sliding window of ${\mathrm{N}\mathrm{M}\mathrm{S}}_w$ frames and retain an extremum only if it is the largest (same polarity) within the window. The final keyframe set is

$$\mathcal{K}={\mathrm{N}\mathrm{M}\mathrm{S}}_w\left({\mathrm{E}}^{†}\right)$$

The effect of persistence thresholding and NMS is illustrated in Figure 2, which contracts dense peak clusters into a sparse, stable keyframe set $\mathcal{K}$. The black solid curve is the original saliency signal $\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}(t)$; the dashed curve is the simplified signal after persistence-based filtering (elbow ${\alpha }^{\star }$ detected by Kneedle) and non-maximum suppression (window $\mathit{w}$). Green crosses mark the final retained keyframes $\mathcal{K}$. The procedure removes clustered minor extrema while preserving structural landmarks that drive segmentation and alignment.

With $\mathcal{K}$ in place, the reliability of subsequent segmentation and HSMM alignment still depends on the saliency weights $w$; Section 2.2.4 details a label-free, consistency-driven calibration.

2.2.4. Adaptive Feature-Weight Learning

In (2.6), the saliency $\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}(t;w)$ is a linear fusion of four heterogeneous features. Assigning heuristic, fixed weights $w={\left[w_v,w_a,w_k,w_d\right]}^{\mathrm{\top }}$ typically overfits a particular operator, platform, or task, degrading segmentation under distribution shift. We therefore treat $w$ as model parameters to be estimated without labels by enforcing cross-demonstration consistency.

• Consistency functiona.

For $M$ demonstrations ${\left\{{\mathbf{P}}^{(m)}\right\}}_{m=1}^M$ and a candidate $w$, define

$$\mathcal{F}(\cdot ;\mathit{w}):{\mathbf{P}}^{(m)}⟼{\widehat{\mathbf{P}}}^{(m)}(\mathit{w}),$$

(2.8)

denote the composition of saliency construction, keyframe extraction (Figure 2), HSMM-based alignment (Section 2.3), and resampling on shared semantic nodes. If $w$ is well chosen, the reconstructions $\left\{{\widehat{P}}^{(m)}(w)\right\}$ should be congruent in shape and timing. We quantify this by the mean point-to-point structural dispersion.

$$\mathrm{S}\mathrm{O}\mathrm{D}(\mathit{w})={\displaystyle \frac{2}{M(M-1)}}{\sum }_{1\le m<m^{\text{'}}\le M} {\displaystyle \frac{1}{T_{min}}}{\sum }_{t=1}^{T_{min}} {‖{\widehat{\mathbf{p}}}_t^{(m)}(w)-{\widehat{\mathbf{p}}}_t^{\left(m^{\text{'}}\right)}(w)‖}_2$$

(2.9)

here $T_{min}$ is the minimal resampled length. Smaller SOD means lower structural variance on the shared semantic time base.

2.

Objective and constraints.

To prevent dominance by any single channel and to improve identifiability, we regularize with a light ${\mathcal{l}}_2$ penalty and minimize over the probability simplex:

$$w^{\star }=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\mathit{w}\in {\mathsf{\Delta }}^3}{min} \left(\mathrm{S}\mathrm{O}\mathrm{D}(\mathit{w})+\lambda \Vert \mathit{w}{\Vert }_2^2\right),{\mathsf{\Delta }}^3=\left\{\mathit{w}\mid w_i\ge 0,\sum _i w_i=1\right\},$$

(2.10)

Because F contains non-smooth steps (extrema detection, discrete decoding), $\mathbf{S}\mathbf{O}\mathbf{D}(\mathit{w})$ is non-differentiable and gradient methods are unreliable.

3.

Solver and feasibility.

We employ the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) in an unconstrained space together with a softmax reparameterization

$$w_i(z)={\displaystyle \frac{\exp \left(z_i\right)}{{\sum }_{j=1}^4 \mathrm{e}\mathrm{x}\mathrm{p}\left(z_j\right)}},z\in {\mathbb{R}}^4$$

which enforces $w\in {\mathsf{\Delta }}^3$ at every iteration while preserving Gaussian updates in $z$-coordinates. CMA-ES is appropriate here because it (i) requires only function values of $J(w)=\mathrm{S}\mathrm{O}\mathrm{D}(w)+\lambda \Vert w{\Vert }_2^2$; (ii) adapts search scales via the evolving covariance; and (iii) is empirically robust on multi-modal, non-convex objectives [34,35]. We terminate when either the covariance trace shrinks below a preset ratio (indicating a local neighborhood) or the successive improvement $\left|J^{(k)}-J^{(k-1)}\right|$ falls below $\epsilon$. Under mild regularity on $J$ in the $z$-space, the iterates approach first-order stationarity (a near-KKT solution) with high probability [36].

4.

Computational profile.

Each evaluation of $J(w)$ entails one pass of keyframe extraction (linear in sequence length $T$) and one HSMM forward–backward/decoding pass per demonstration with complexity $O\left(NT{D}_{max}\right)$; overall cost scales linearly in $M$. For numerical stability we implement recursions in the log domain (log-sum-exp), apply a diagonal floor to GMM covariances, and cache feature streams across outer-loop calls.

2.3. Multi-Demo Alignment and Segmentation via a Duration-Explicit HSMM

Given the sparse, scale-invariant keyframes ${\mathcal{K}}^{(m)}={\left\{t_k^{(m)}\right\}}_{k=1}^{K_m}$ obtained from saliency (Section 2.2) and the learned fusion weights (Section 2.2.4), we seek a shared semantic time base across demonstrations. Directly matching wall-clock indices is unreliable due to operator-dependent pauses and backtracking; instead, we align demonstrations probabilistically via a Hidden semi-Markov model (HSMM) that explicitly models state durations (Figure 3) [37,38].

2.3.1. Model and Generative Mechanism

Let $\mathcal{S}=\left\{S_1,S_2,\dots ,S_N\right\},N\ll K_m$ be latent phases, each representing a macro-action (e.g., grasp, insert, lift). With parameters $\mathsf{\Theta }=\left\{\pi ,A,\left\{p_i\right\},\left\{b_i\right\}\right\}$, the HSMM generates $o_t\in {\mathbb{R}}^d$ at $t=1,\dots ,T$ as follows:

a)

Initial phase: $q_1\sim \mathrm{C}\mathrm{a}\mathrm{t}(\pi ),{\sum }_i {\pi }_i=1.$

b)

Duration: for the current phase $q_k$, sample dwell length $d_k\sim p_{q_k}(d)\in {\mathbb{N}}_{>0}$.

c)

Observations: for ${\tau }_{k-1}<t\le {\tau }_k$ with ${\tau }_k={\sum }_{j\le k} d_j$, $o_t\sim b_{q_k}(\cdot ),b_i(o)={\sum }_{m=1}^{M_i} {\pi }_{i,m}\mathcal{N}\left(o\mid {\mu }_{i,m},{\mathsf{\Sigma }}_{i,m}\right)$

d)

Transition: $\mathrm{P}\mathrm{r}\left(q_{k+1}=j\mid q_k\right)=a_{q_k}j$; terminal transitions end the sequence.

Writing $O=\left(o_1,\dots ,o_T\right),Q=\left(q_1,\dots ,q_n\right),D=\left(d_1,\dots ,d_n\right),$and $q^{*}(t)=q_k$ for ${\tau }_{k-1}<t\le {\tau }_k$, the joint density is

$$P(O,Q,D)={\pi }_{q_1}{p}_{q_1}\left(d_1\right)\prod _{k=1}^{n-1} \left[a_{q_k{q}_{k+1}}{p}_{q_{k+1}}\left(d_{k+1}\right)\right]\prod _{t=1}^T b_{q^{*}(t)}\left(o_t\right).$$

(2.11)

• Observation design. We concatenate kinematic descriptors (e.g., curvature, speed magnitude, DCR) and any available modalities (pose, force/tactile, depth) into $o_t$.
• Duration choices. We use either discrete $p_i(d)$ with support 1: $D_{max}$, or truncated Gaussian/Gamma families to accommodate unequal dwell times [39].

2.3.2. Parameter Estimation: Extended Baum–Welch

We maximize the total log-likelihood over demonstrations

$$\mathcal{L}(\mathsf{\Theta })=\sum _{m=1}^M \mathrm{l}\mathrm{o}\mathrm{g}P\left(O^{(m)}\mid \mathsf{\Theta }\right).$$

(2.12)

Explicit durations break first-order Markovity; forward–backward recursions therefore enumerate duration indices.

1. 

Forward variable (leaving $\mathit{i}$ at time $\mathit{t}$ ).

$${\alpha }_t(i)=\sum _{j=1}^N \sum _{d=1}^{min\left(D_{max},t\right)} {\alpha }_{t-d}(j)a_{ji}{p}_i(d)\prod _{r=0}^{d-1} b_i\left(o_{t-r}\right),{\alpha }_0(i)={\pi }_i.$$

(2.13)

2. 

Backward variable.

$${\beta }_t(i)=\sum _{j=1}^N \sum _{d=1}^{T-t} \left[a_{ij}{p}_j(d)\prod _{r=0}^{d-1} b_j\left(o_{t+1+r}\right){\beta }_{t+d}(j)\right],{\beta }_T(i)=1.$$

(2.14)

3. 

Posteriors (E-step).

$${\gamma }_t(i)={\displaystyle \frac{{\sum }_{d=1}^{min\left(D_{max},t\right)} {\sum }_{r=0}^{d-1} {\alpha }_{t-r}(i)p_i(d){\prod }_{u=0}^{d-1} b_i\left(o_{t-r+u}\right){\beta }_{t-r}(i)}{{\sum }_{k=1}^N {\alpha }_T(k)}}$$

(2.15)

$${\xi }_{t,d}(i,j)={\displaystyle \frac{{\alpha }_{t-d}(i)a_{ij}{p}_j(d){\prod }_{r=0}^{d-1} b_j\left(o_{t-r}\right){\beta }_t(j)}{{\sum }_{k=1}^N {\alpha }_T(k)}}$$

(2.16)

Sum $\gamma ,\xi$ over $m$ to obtain corpus-level sufficient statistics.

4. 

M-step (closed forms).

$$\begin{array}{cc}\hfill {\pi }_i^{\mathrm{new}}& ={\gamma }_{t=1}(i),\hfill \\ \hfill a_{ij}^{\mathrm{new}}& ={\displaystyle \frac{{\sum }_{t,d} {\xi }_{t,d}(i,j)}{{\sum }_t {\gamma }_t(i)}},\hfill \\ \hfill p_i^{\mathrm{new}}(d)& ={\displaystyle \frac{{\sum }_t {\xi }_{t,d}(i,\cdot )}{{\sum }_{t,d^{\text{'}}} {\xi }_{t,d^{\text{'}}}(i,\cdot )}},\hfill \\ \hfill {\mu }_{i,m}^{\mathrm{new}}& ={\displaystyle \frac{{\sum }_t {\gamma }_{t,m}(i)o_t}{{\sum }_t {\gamma }_{t,m}(i)}},\hfill \\ \hfill {\mathsf{\Sigma }}_{i,m}^{\mathrm{new}}& ={\displaystyle \frac{\sum _t {\gamma }_{t,m}(i)\left(o_t-{\mu }_{i,m}\right){\left(o_t-{\mu }_{i,m}\right)}^{\mathrm{\top }}}{\sum _t {\gamma }_{t,m}(i)}}\hfill \end{array}$$

(2.17)

with ${\gamma }_{t,m}(i)={\gamma }_t(i){\pi }_{i,m}\mathcal{N}\left(o_t\mid {\mu }_{i,m},{\mathsf{\Sigma }}_{i,m}\right)/b_i\left(o_t\right)$

5. 

Numerical stability.

All recursions are implemented in the log domain using log-sum-exp

$$\mathrm{L}\mathrm{S}\mathrm{E}\left(x_1,\dots ,x_K\right)=\mathrm{l}\mathrm{o}\mathrm{g}\sum _{k=1}^K e^{x_k},$$

(2.18)

and GMM covariances receive a diagonal floor $\delta I$. EM monotonicity guarantees convergence of${\mathcal{L}}^{(k)}$ to a stationary point [40,41].

2.3.3. Semantic Time Axis: Decoding and Outputs

After convergence, Viterbi decoding yields the MAP phase path ${\widehat{Q}}^{(m)}=\left(q_1^{(m)},\dots ,q_N^{(m)}\right)$ and durations ${\widehat{D}}^{(m)}=\left(d_1^{(m)},\dots ,d_N^{(m)}\right)$ for each demonstration, with ${\sum }_{k=1}^N d_k^{(m)}=T_m$. Define cumulative boundaries ${\tau }_k^{(m)}={\sum }_{j=1}^k d_j^{(m)},k=1,\dots ,N$, so segments $\left({\tau }_{k-1}^{(m)},{\tau }_k^{(m)}\right]$ correspond to the same semantic phase $S_k$ across demonstrations. Figure 4 and Figure 5 illustrate 3-D keyframe locations and HSMM reconstructions aligned to salient kinematic transitions.

2.3.4. Alignment Quality, Model Selection, and Robustness

Alignment metric. We quantify cross-demo temporal agreement by

$$\mathrm{A}\mathrm{E}(N)=\sum _{m=1}^M {‖\left({\tau }_1^{(m)},\dots ,{\tau }_N^{(m)}\right)-\left({\overline{\tau }}_1,\dots ,{\overline{\tau }}_N\right)‖}_2,{\overline{\tau }}_k={\displaystyle \frac{1}{M}}\sum _m {\tau }_k^{(m)}.$$

(2.19)

Model selection. To balance fit and parsimony, we use

$$\mathrm{B}\mathrm{I}\mathrm{C}(N)=-2\mathrm{l}\mathrm{o}\mathrm{g}\mathcal{L}\left({\mathsf{\Theta }}_N^{\star }\right)+{\kappa }_N\mathrm{l}\mathrm{o}\mathrm{g}\left(\sum _m T_m\right),$$

(2.20)

and select $N^{\star }$ by jointly considering $\mathrm{A}\mathrm{E}(N)$ and $\mathrm{B}\mathrm{I}\mathrm{C}(N)$.

Robustness to time warps. For any order-preserving reparameterization of time, $t\mapsto \tau (t)$, the decoded phase order $\left(S_1,\dots ,S_N\right)$ is invariant; only the tail behavior of duration distributions $p_i(d)$ is rescaled. This accounts for operator-dependent slowdowns, hesitations, or hovering.

Complexity. Keyframe processing is $O(T)$ per sequence; duration-explicit forward–backward is $O\left(NT{D}_{max}\right)$; training scales linearly in the number of demonstrations M. Outputs include $\left\{S_k\right\}$, and segment-wise statistics for direct use by planners (e.g., OMPL/MPC).

2.4. Statistical Motion Primitives and Probabilistic Generation

With cross-demonstration phases and boundaries ${\left\{{\mathcal{S}}_k\right\}}_{k=1}^N$, $\left\{{\tau }_k^{(m)}\right\}$obtained from the duration-explicit HSMM (Section 2.3), we model each phase within the shared semantic time base. For demonstration $m$, let

$${\mathcal{C}}_k^{(m)}=\left\{P^{(m)}(t);t\in \left({\tau }_{k-1}^{(m)},{\tau }_k^{(m)}\right]\right\}$$

denote the spatio-temporal segment for phase ${\mathcal{S}}_k$. Our objective is a segment-wise generative model that (i) captures cross-demo variability, (ii) supports conditioning and duration re-scaling, and (iii) yields calibrated uncertainty for downstream planning and safety assessment. We instantiate three complementary families: DMP, GMM/GMR, and ProMP. Unless otherwise noted, time within a segment is normalized to a phase variable $s\in \left[\mathrm{0,1}\right]$, ensuring a uniform interface across models.

2.4.1. Dynamic Movement Primitives (DMP)

• Single-segment dynamics.

For a one-degree-of-freedom trajectory $x(t)$, the classical DMP represents motion as a critically damped second-order system with a nonlinear forcing term:

$$\begin{array}{l}\tau \dot{z}={\alpha }_z\left({\beta }_z(g-x)-z\right)+f(s),\\ \tau \dot{x}=z,\tau \dot{s}=-{\alpha }_{s}s,\end{array}$$

(2.21)

where $g$ is the segment goal, $s\in (\mathrm{0,1}]$ is a phase variable, and

$$f(s)={\displaystyle \frac{{\sum }_{i=1}^B w_i{\psi }_i(s)}{{\sum }_{i=1}^B {\psi }_i(s)}}\left(g-x_0\right),{\psi }_i(s)=\mathrm{e}\mathrm{x}\mathrm{p}\left[-h_i{\left(s-c_i\right)}^2\right]$$

(2.22)

Given $\tau ,{\alpha }_z,{\beta }_z,{\alpha }_{s^{\text{'}}}$ the weights $\left\{w_i\right\}$are obtained by least squares (or locally weighted regression). Multi-DOF trajectories are modeled component-wise or via task-space coupling.

2.

Segment coupling and smoothness.

Let ${\overline{d}}_k$ be the mean duration of phase $S_k$ across demonstrations; set ${\tau }_k={\overline{d}}_k\mathsf{\Delta }t$. We fit ${\left\{w_i^{(k)}\right\}}_{i=1}^B$ per segment and compose ${\left\{\left({\mathrm{D}\mathrm{M}\mathrm{P}}_k,{\tau }_k\right)\right\}}_{k=1}^N$ along the decoded boundaries. Because the state$(x,z)$ is continuous across boundaries, the concatenation is $C^1$ -continuous without auxiliary velocity/acceleration matchers. DMPs thus offer low jerk and robust time-scaling at execution.

2.4.2. Gaussian Mixture Modeling and Regression (GMM/GMR)

Mixture modeling. For each segment, we pair normalized phase with position, $\left(s_t,y_t\right)\in {\mathbb{R}}^{1+d}$, and fit a K-component mixture

$$p(\mathbf{s},\mathbf{y})=\sum _{k=1}^K {\pi }_k\mathcal{N}\left({\mu }_k,{\mathsf{\Sigma }}_k\right)$$

(2.23)

typically with diagonal (or block-diagonal) ${\mathsf{\Sigma }}_k$ to denoise while retaining principal correlations.

Regression and uncertainty. At execution, Gaussian Mixture Regression (GMR) produces the conditional

$$p(\mathbf{y}\mid \mathbf{s})=\mathcal{N}(\widehat{\mu }(\mathbf{s}),\widehat{\mathsf{\Sigma }}(\mathbf{s})),$$

yielding a closed-form mean trajectory and a phase-indexed covariance band. Analytical derivatives of $\widehat{\mu }\left(\mathbf{s}\right)\mathrm{a}\mathrm{n}\mathrm{d}\widehat{\mathsf{\Sigma }}(\mathbf{s})$ facilitate integration with MPC for online terminal corrections and constraint handling.

2.4.3. Probabilistic Movement Primitives (ProMP)

Bayesian representation. Using basis functions ${\mathsf{\Phi }}_t={\left[{\varphi }_1(t),\dots ,{\varphi }_B(t)\right]}^{\mathrm{\top }}$, a segment trajectory is modeled as

$${\mathit{y}}_t={\Phi }_t^{\top }\mathit{w}+\epsilon ,\epsilon \sim \mathcal{N}\left(0,{\Sigma }_y\right),\mathit{w}\sim \mathcal{N}\left({\mu }_w,{\Sigma }_w\right)$$

(2.24)

where $\left({\mu }_w,{\Sigma }_w\right)$ encode the distribution over shapes learned from multiple demonstrations (EM or closed-form under Gaussian assumptions).

Conditioning and coordination. Linear constraints—endpoints, waypoints, or partial observations—are imposed by conditioning the weight distribution:

$$w\mid \left(C{Y}_T=g\right)\sim \mathcal{N}\left({\mu }_w+K_T\left(g-{\mathsf{\Phi }}_T{\mu }_w\right),{\mathsf{\Sigma }}_w-K_T{\mathsf{\Phi }}_T{\mathsf{\Sigma }}_w\right),$$

(2.25)

with $K_T={\mathsf{\Sigma }}_w{\mathsf{\Phi }}_T{\left({\mathsf{\Phi }}_T^{\mathrm{\top }}{\mathsf{\Sigma }}_w{\mathsf{\Phi }}_T+{\mathsf{\Sigma }}_{\epsilon }\right)}^{-1}$. Sampling from this posterior and stitching segments at decoded boundaries yields constraint-consistent trajectories with explicit predictive uncertainty.

2.4.4. Model Choice and Complementarity

No single primitive dominates: DMP excels in real-time low-jerk execution; GMM/GMR gives closed-form means/covariances and gradients; ProMP supports exact linear-Gaussian conditioning for multi-goal tasks. Table 1 summarizes properties. Figure 6 shows HSMM-aligned ProMP generation with calibrated uncertainty bands:

• DMP excels at real-time execution and low jerk with simple time scaling;
• GMM/GMR offers closed-form means and covariances over phase and is convenient for planners needing analytic gradients;
• ProMP provides a distribution over shapes with exact linear-Gaussian conditioning, ideal for multi-goal tasks and collaboration.

In practice, we first fit GMM/GMR to obtain mean/covariance; use the GMR mean to initialize DMP weights for smooth execution; and, when hard/soft constraints or multi-goal adaptation are required, overlay ProMP conditioning on top of the DMP nominal to reconcile smoothness with constraints.

Interface to planning and safety. Segment-wise covariances (from GMR or ProMP) propagate to risk metrics and constraint tightening in MPC/OMPL. Because segments share a semantic time base, the uncertainty is comparable across operators and platforms, enabling principled safety margins (e.g., $2\sigma$ envelopes) and priority-aware blending when composing multi-segment tasks.

Computational profile. Per segment, DMP fitting is linear in the number of basis functions; GMM/GMR training is $O(KT)$ per EM iteration; ProMP estimation is closed-form/EM with state dimension B. Since segments are independent given HSMM alignment, training parallelizes over phases and scales linearly with the number of demonstrations.

3. Experiments and Results

3.1. Objectives and Evaluation Protocol

This section subjects the proposed end-to-end pipeline to a cross-domain, reproducible evaluation covering the three core components introduced in Section 2: (i) unsupervised segmentation (multi-feature saliency with topological persistence), (ii) duration-explicit alignment (HSMM), and (iii) probabilistic in-segment generation (GMM/GMR, ProMP, optionally DMP). We target three complementary questions:

• Segmentation robustness. Do multi-feature saliency and topological persistence yield sparse yet structurally stable keyframes under heterogeneous noise and tempo variations?
• Semantic alignment quality. Does duration-explicit HSMM reduce cross-demonstration time dispersion when non-geometric dwelling is present (e.g., hover, wait)?
• Generator calibration. On the shared semantic time base, do segment-wise probabilistic models achieve low reconstruction error, nominal uncertainty coverage, and dynamically schedulable executions?

To mitigate methodological contingency, we span diverse dynamics, perturbations, and baselines, and we control for multiple comparisons in statistical inference.

3.2. Tasks, Datasets, and Testable Hypotheses

We consider three representative domains:

• Domain A—UAV-Sim (multi-scene flight). 100 Hz sampling. Subtasks include take-off–lift–cruise–drop and gate-pass–loiter–gate-pass. Six subjects, 20–30 segments per task. Observations: tool-center position (optional yaw). Figure 7 shows the environment and demonstrations.
• Domain B—AV-Sim (CARLA/MetaDrive urban). 10 Hz sampling across Town01–05, varied weather/lighting and traffic control. Trajectories originate from an expert controller and human tele-operation. Observations: $(x,y,\theta )$. See Figure 8.
• Domain C—Manip-Sim (robomimic/RLBench assembly). 50–100 Hz sampling. Tasks akin to RoboTurk “square-nut”: grasp–align–insert with pronounced dwell segments. Observations: end-effector position. See Figure 9.

3.3. Metrics and Statistical Inference

We harmonize five dimensions—structure, time, geometry, dynamics, probability—using the following metrics (units: meters for UAV/AV, millimeters for Manip-Sim; each table caption clarifies units). Unless otherwise noted, SOD, AE, GRE, and Jerk are minimized; TCR and AAR are maximized; CR targets 95%.

• SOD (Eq. (2.10), min): structural dispersion—mean point-to-point divergence on the shared time base.
• AE (Eq. (2.17), min): Euclidean dispersion of phase end times across demonstrations.
• AAR (max): action acquisition rate. Given reference key actions $\left\{{\overline{\tau }}_k\right\}$ (expert consensus / common boundaries) and detected $\left\{{\tau }_k^{\text{'}}\right\}$, we count a hit if $\left|{\tau }_k^{\text{'}}-{\overline{\tau }}_k\right|\le \delta$, with $\delta =max(5\mathrm{frames},0.2\mathrm{s})$ (i.e., 50 ms at 100 Hz; 0.5 s at 10 Hz).
• GRE (min): geometric reconstruction error (RMSE).
• TCR (max): time compression rate.
• Jerk (min): $J=\sum _t \Vert \stackrel{⃛}{P}(t){\Vert }^2$ (normalized).
• CR (target ≈ 95%): nominal $2\sigma$ coverage. For each semantic segment we sample uniformly in time; if $y_t\in \left[{\mu }_t\pm 2\sqrt{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\mathsf{\Sigma }}_t\right)}\right]$, it is counted as covered; segment-level coverage is averaged and then length-weighted globally.

Statistical inference. We apply Shapiro–Wilk normality and Levene homoscedasticity tests. When met, paired t-tests are used; otherwise Wilcoxon signed-rank tests. Multiple comparisons are controlled via Holm–Bonferroni. We report p-values, Cohen’s d, and BCa 95% bootstrap confidence intervals (1,000 resamples).

3.4. Baselines and Fairness Controls

We benchmark four families to isolate contributions at the signal, boundary, and latent levels:

(i) Signal-level: single-feature curvature + quantile threshold; and multi-feature equal weights without TDA simplification / NMS (no persistence, no cluster suppression).

(ii) Boundary-level: BOCPD (Bayesian Online Changepoint Detection).

(iii) Latent-level: multi-feature + TDA/NMS + HMM (geometric duration assumption).

(iv) Full method: $w^{\star }$ (consistency-learned weights) + TDA/NMS + HSMM (duration explicit) + segment-wise generator (default ProMP; we also compare DMP/GMR on the same segmentation when isolating generation quality, § 3.9).

Fairness controls. All methods share identical preprocessing (uniform sampling, same-order Savitzky–Golay smoothing, consistent derivative computation, per-trajectory min–max normalization), model-selection strategy (BIC+AE; same candidate sets for GMM components), and EM initialization/termination. Segmentation quality is evaluated using each model’s MAP/Viterbi boundaries; for pure generator comparison (§ 3.9), we fix HSMM boundaries across methods to remove boundary confounds.

3.5. Overall Results

Cross-domain evidence appears in Figure 10 and Table 2. Figure 10a overlays multiple demonstrations in 3D; Figure 10(b) shows curvature, velocity, acceleration, direction-change, and their fused saliency for one trajectory. Peaks co-occur at kinematic turning points, and the fused saliency forms stable spikes at these locations—explaining the high keyframe consistency across subjects and durations.

• Domain A—UAV-Sim (Table 2)

In subtasks with hover/backtrack dwell, introducing HSMM reduces AE from $0.41\pm 0.09$ s to $0.28\pm 0.07$ s (−31%, p < 0.01, Cohen’s $d=0.86$). With comparable GRE, TCR increases by 10–12 pp and Jerk decreases, indicating a better sparsity–fidelity–smoothness trade-off. The probabilistic output attains CR ≈ 95% at the nominal $2\sigma$ level.

• Domain B—AV-Sim (Table 3)

In segments with non-geometric dwell (e.g., slow-down–wait–turn), AE drops from $0.37\pm 0.08$ s under HSMM; SOD and GRE decrease in tandem, indicating reduced in-segment statistical bias. ProMP coverage at $2\sigma$ is $95.1\pm 2.5\mathrm{\%}$.

Table 3. AV-Sim (CARLA/MetaDrive) (units: SOD/GRE in meters).

Method	SOD (m)	AE (s)	GRE (m)	TCR (%)	AAR (%)	Jerk	CR (%)
Curvature + quantile	0.172 ± 0.030	0.70 ± 0.14	0.247 ± 0.041	32.7 ± 5.2	66.0 ± 6.7	1.18 ± 0.09	–
Multi-feat (equal), no TDA/NMS	0.160 ± 0.029	0.63 ± 0.12	0.231 ± 0.038	39.5 ± 5.0	71.6 ± 6.1	1.14 ± 0.08	–
Multi-feat + TDA/NMS + HMM	0.148 ± 0.027	0.55 ± 0.11	0.214 ± 0.035	47.4 ± 4.7	78.8 ± 5.4	1.08 ± 0.07	–
Ours: $w^{\star }$+TDA/NMS+HSMM+ProMP	0.112 ± 0.022	0.37 ± 0.08	0.191 ± 0.033	47.5 ± 4.6	86.3 ± 5.0	1.00 ± 0.06	95.1 ± 2.5

Table 3. AV-Sim (CARLA/MetaDrive) (units: SOD/GRE in meters).

Method	SOD (m)	AE (s)	GRE (m)	TCR (%)	AAR (%)	Jerk	CR (%)
Curvature + quantile	0.172 ± 0.030	0.70 ± 0.14	0.247 ± 0.041	32.7 ± 5.2	66.0 ± 6.7	1.18 ± 0.09	–
Multi-feat (equal), no TDA/NMS	0.160 ± 0.029	0.63 ± 0.12	0.231 ± 0.038	39.5 ± 5.0	71.6 ± 6.1	1.14 ± 0.08	–
Multi-feat + TDA/NMS + HMM	0.148 ± 0.027	0.55 ± 0.11	0.214 ± 0.035	47.4 ± 4.7	78.8 ± 5.4	1.08 ± 0.07	–
Ours: $w^{\star }$+TDA/NMS+HSMM+ProMP	0.112 ± 0.022	0.37 ± 0.08	0.191 ± 0.033	47.5 ± 4.6	86.3 ± 5.0	1.00 ± 0.06	95.1 ± 2.5

• Domain C—Manip-Sim (Table 4)

Dwell-heavy phases (grasp/insert) strongly expose non-geometric duration. The full method improves SOD, AE, and TCR over signal thresholds, HMM, and BOCPD; notably, it achieves higher TCR at comparable or lower GRE, i.e., fewer keyframes suffice to reconstruct high-fidelity shapes.

Table 4. Manip--Sim (robomimic/RLBench) (units: SOD/GRE in millimeters).

Method	SOD (mm)	AE (s)	GRE (mm)	TCR (%)	AAR (%)
Curvature + quantile	1.12 ± 0.27	0.33 ± 0.09	1.02 ± 0.19	30.0 ± 4.8	65.3 ± 6.3
Multi-feat (equal), no TDA/NMS	0.94 ± 0.24	0.28 ± 0.08	0.86 ± 0.16	33.0 ± 5.0	69.1 ± 5.8
Multi-feat + TDA/NMS + HMM	0.87 ± 0.22	0.29 ± 0.08	0.91 ± 0.17	35.1 ± 5.1	70.0 ± 5.7
BOCPD	0.98 ± 0.25	0.31 ± 0.09	1.05 ± 0.20	45.0 ± 5.6	71.9 ± 5.4
Ours: $w^{\star }$+TDA/NMS+HSMM +ProMP	0.72 ± 0.18	0.24 ± 0.07	0.79 ± 0.15	49.2 ± 4.8	83.4 ± 5.1

Table 4. Manip--Sim (robomimic/RLBench) (units: SOD/GRE in millimeters).

Method	SOD (mm)	AE (s)	GRE (mm)	TCR (%)	AAR (%)
Curvature + quantile	1.12 ± 0.27	0.33 ± 0.09	1.02 ± 0.19	30.0 ± 4.8	65.3 ± 6.3
Multi-feat (equal), no TDA/NMS	0.94 ± 0.24	0.28 ± 0.08	0.86 ± 0.16	33.0 ± 5.0	69.1 ± 5.8
Multi-feat + TDA/NMS + HMM	0.87 ± 0.22	0.29 ± 0.08	0.91 ± 0.17	35.1 ± 5.1	70.0 ± 5.7
BOCPD	0.98 ± 0.25	0.31 ± 0.09	1.05 ± 0.20	45.0 ± 5.6	71.9 ± 5.4
Ours: $w^{\star }$+TDA/NMS+HSMM +ProMP	0.72 ± 0.18	0.24 ± 0.07	0.79 ± 0.15	49.2 ± 4.8	83.4 ± 5.1

Semantic time alignment. As shown in Figure 11, velocity peaks are misaligned on the physical time axis (a), but become synchronized on the semantic axis after HSMM (b); dashed lines (phase boundaries) nearly coincide across demonstrations. This mirrors the systematic AE reduction in Table 2, Table 3 and Table 4 and evidences the benefit over geometric-duration HMM.

3.6. Contribution Attribution: Ablation Study (UAV-Sim)

We perform stepwise ablations under LOSO with $N=6$ (other settings as in § 3.3–3.4). Table 5 reports relative changes w.r.t. the full method ($w^{\star }$+TDA/NMS+HSMM+ProMP); positive values are worse (e.g., AE↑), negative better.

• Remove TDA (keep NMS): SOD +12.3%, AE +9.6% → persistence is key to scale-invariant noise rejection; without it, small-scale oscillations stack into spurious peak–valley pairs, degrading structure and boundaries.
• Remove NMS (keep TDA): AE +21.1% → suppressing same-polarity peak clusters in high-energy regions is critical for boundary stability; persistence alone cannot prevent multi-response.
• Fix equal weights (no ${\mathit{w}}^{\star }$): SOD +18.6%, AAR −7.7 pp → consistency-driven weight learning mitigates channel scale imbalance and improves key-action capture.
• Replace HSMM with HMM: AE +31.2% → direct evidence of the geometric-duration bias when dwell exists (wait/loiter).

Differences on AE and SOD remain statistically significant after Holm–Bonferroni (p < 0.05); effect sizes are medium-to-large. In sum: TDA+NMS stabilize the input structure, $w^{\star }$ provides cross-demo self-calibration, and HSMM addresses duration bias mechanistically.

3.7. Robustness: Time-Warping, Noise, and Missing Data (AV-Sim)

We evaluate robustness on AV-Sim (Table 6. We apply monotone time warps $t\mapsto t^{0.8}$ (speed-up) and $t\mapsto t^{1.2}$ (slow-down), add white noise $\sigma =0.05\mathrm{m}$ to positions, and randomly drop 20% of samples. Compared to HMM, HSMM reduces AE by 36%, 34%, 30%, and 35%, respectively, while preserving phase order. AE grows roughly linearly with noise amplitude but, under order-preserving time warps, manifests mainly as duration redistribution—consistent with the “weak time-distortion invariance” discussed in § 2.4.

ProMP $2\sigma$ coverage stays within 94–96% under these perturbations, indicating stable uncertainty calibration.

3.8. Model Selection and the Sparsity–Fidelity Trade-off

Figure 12 shows how the persistence threshold $\alpha$ controls keyframe sparsity and reconstruction GRE. As $\alpha$ increases, keyframe count decreases monotonically; GRE stabilizes near the Kneedle elbow ${\alpha }^{\star }$ (vertical dashed line), reflecting an effective sparsity–fidelity compromise with BCa 95% CIs (shaded).

To determine HSMM state number $N$, we train models on $N\in \{2,\dots ,8\}$ with identical hyper-parameters and initialization. On the validation set we compute AE and BIC (Eq. (2.17) and Eq. (2.20)), normalize both to $\left[\mathrm{0,1}\right]$, and average across subjects. Figure 13 shows that AE decreases rapidly then plateaus at $N\ge 6$, whereas BIC increases monotonically. The joint criterion $J(N)=\widetilde{\mathrm{A}\mathrm{E}}(N)+\widetilde{\mathrm{B}\mathrm{I}\mathrm{C}}(N)$ yields the optimal number of states $\widehat{N}$ or UAV/AV and 5 for Manip-Sim (star markers). Deviating by ±1 does not change the overall trend; we fix $N^{\star }$ in all experiments to balance alignment accuracy and parsimony.

3.9. In-Segment Generation: Accuracy, Smoothness, and Calibration

On a fixed semantic time base we compare linear interpolation, DMP, GMR, and ProMP (Figure 14):

• Geometry: ProMP ≈ GMR < DMP < linear in GRE.
• Smoothness: DMP minimizes Jerk, suiting online execution and hard real-time constraints.
• Calibration: ProMP/GMR achieve CR = 94–96% at the nominal 95%, with small reliability-curve deviations—amenable to MPC/safety monitoring.

Segments are concatenated at HSMM boundaries with $C^1$ continuity (see Eqs. (2.19)–(2.20)); the critically damped second-order dynamics preserve velocity/acceleration continuity across segments.

3.10. Summary of Findings

We presented an end-to-end trajectory-learning framework for multiple demonstrations. Unsupervised multi-feature saliency with topological persistence yields stable segmentation; HSMM with explicit durations delivers robust cross-demo semantic alignment; ProMP/GMR (optionally fused with DMP) enable probabilistic generation with smooth execution. Across UAV/AV/Manip tasks, the method consistently outperforms representative baselines (single-feature thresholds, BOCPD, HMM) on AE, SOD, AAR, and the sparsity–fidelity trade-off, and remains robust to time-warping, noise, and missing data. Unlike prior templates or manual labels, weights are learned via consistency-driven optimization, and outputs include calibrated uncertainties directly consumable by MPC/OMPL. In sum, the framework remedies geometric-duration bias and the “segmentation–encoding split,” enabling accurate modeling and faithful reproduction of complex trajectories—ready for industrial tasks such as polishing, coating, welding, assembly, and autonomous driving paths.

4. Discussion

• Evidence vs. hypotheses.

• (H1) Segmentation. The topology-aware saliency (multi-feature + persistence + NMS) yields sparse yet stable anchors: removing persistence or NMS increases SOD/AE by 9–21% and lowers AAR (ablation), confirming that both scale-invariant pruning and peak-cluster suppression are necessary.
• (H2) Alignment. Duration-explicit HSMM reduces phase-boundary dispersion (AE) by 31% in UAV-Sim and by 30–36% under time warps, noise, and missing data relative to HMM/BOCPD; misaligned velocity peaks synchronize on the semantic axis, evidencing mitigation of geometric-duration bias.
• (H3) Generation. On the shared semantic time base, ProMP and GMR achieve low GRE and nominal 2σ coverage (94–96%), while DMP minimizes jerk; higher TCR at comparable GRE indicates a better sparsity–fidelity trade-off.

Positioning to prior work. Rule- or template-based segmentation and BOCPD lack geometric/semantic guarantees; DTW aligns sequences but offers no generative dwell model; HMMs assume geometric durations. Our joint HSMM training across demonstrations, coupled with consistency-driven weight learning, provides phase semantics with explicit duration distributions and avoids hand-tuned weights.

Practical implications. Segment-wise covariances (GMR/ProMP) propagate directly to MPC/OMPL for risk-aware execution; the entire pipeline is label-free and computationally linear in the number of demonstrations; Figure 12 and Figure 13 give operational choices for the persistence threshold and state number.

Limitations. (i) Weakly featured tasks may yield shallow saliency peaks; (ii) the fixed feature family can miss critical modalities (e.g., tactile); (iii) unimodal duration models may underfit multi-modal dwell; (iv) Kneedle/NMS windows are still heuristics; (v) sim-to-real effects may perturb calibration.

Future work. Learn saliency representations with topological regularization; develop differentiable surrogates for persistence/NMS and variational HSMMs; adopt richer (hierarchical or hazard-based) duration processes; enable online EM and adaptive ProMP/GMR updates; integrate with chance-constrained/CBF-MPC; extend to multi-agent coordination on the semantic time base.

5. Conclusions

We introduced a label-free pipeline that (i) extracts scale-robust keyframes via topology-aware multi-feature saliency, (ii) performs duration-explicit alignment with a jointly trained HSMM to build a shared semantic time base, and (iii) encodes each phase probabilistically (GMR/ProMP, optionally combined with DMP) for smooth, risk-aware execution. Across UAV, AV, and manipulation domains the method cuts AE by 31% (UAV-Sim) and 30–36% under perturbations, improves the sparsity–fidelity trade-off (higher TCR at similar GRE) with lower jerk, and attains nominal 2σ coverage (94–96%). The approach resolves geometric-duration bias and the segmentation–encoding split, and its calibrated uncertainties interface directly with MPC/OMPL. Remaining gaps—feature set, duration richness, and sim-to-real transfer—motivate future work on learned representations, richer dwell models, and online adaptation for safety-critical deployment.