Path Structured Multimarginal Schrödinger Bridge for Probabilistic
Learning of Hardware Resource Usage by Control Software

We emphasize that our proposed method is not limited by what specific components comprise $\bm{\xi}$. To highlight this flexibility, we describe the proposed approach for $\bm{\xi}\in\mathbb{R}^{d}$ with suitable interpretations for the specific application.

For a time interval $[0,t]$ of interest, we think of time-varying $\bm{\xi}$ as a continuous time vector-valued stochastic process over subsets of $\mathbb{R}^{d}$. Suppose that $s\in\mathbb{N},s\geq 2$ snapshots or observations are made for the stochastic state $\bm{\xi}(\tau)$, $0\leq\tau\leq t$, at (possibly non-equispaced) instances

Consider the snapshot index set $\llbracket s\rrbracket:=\{1,\ldots,s\}$. For a fixed context $\bm{c}$, the snapshot observations comprise a sequence of joint probability measures $\{\mu_{\sigma}\}_{\sigma\in\llbracket s\rrbracket}$ satisfying $\int{\rm{d}}\mu_{\sigma}(\bm{\xi}(\tau_{\sigma}))=1$. In other words,

In our application, the data $\{\mu_{\sigma}\}_{\sigma\in\llbracket s\rrbracket}$ comes from control software execution profiles, i.e., by executing the same control software for the same $\bm{c}$ with all parameters and initial conditions fixed. So the stochasticity in $\bm{\xi}(\tau_{\sigma})$ stems from the dynamic variability in hardware resource availability.

In particular, for finitely many (say $n$) execution profiles, we consider empirical distributions

where $\delta(\bm{\xi}-\bm{\xi}^{i}(\tau_{\sigma}))$ denotes the Dirac delta at sample location $\bm{\xi}^{i}(\tau_{\sigma})$ where $i\in\llbracket n\rrbracket$, $\sigma\in\llbracket s\rrbracket$. At any snapshot index $\sigma\in\llbracket s\rrbracket$, the set $\{\bm{\xi}^{i}(\tau_{\sigma})\}_{i=1}^{n}$ is scattered data.

Given the data (8)-(9), we would like to predict the most likely hardware resource state statistics

Without the qualifier “most likely”, the problem is overdetermined since there are uncountably many measure-valued continuous curves over $[0,t]$ that are consistent with the observed data (8)-(9).

Following [6, Sec. 3], let

For $\varepsilon\geq 0$ (not necessarily small), the multimarginal Schrödinger bridge problem (MSBP) is the following infinite dimensional convex program:

In particular, $\mathcal{M}(\bm{\mathcal{X}})$ is a convex set. The objective (12a) is strictly convex in $\bm{M}$, thanks to the $\varepsilon$-regularized negative entropy term $\int_{\bm{\mathcal{X}}}\varepsilon\bm{M}\log\bm{M}$. The constraints (12b) are linear.

In this work, the measures $\{\mu_{\sigma}\}_{\sigma\in\llbracket s\rrbracket}$ correspond to sequential observation, and we therefore fix the path structured (Fig. 1) ground cost

In particular, we choose the squared Euclidean distance sequential cost between two consecutive snapshot indices, i.e., $c_{\sigma}(\cdot,\cdot):=\|\cdot-\cdot\|_{2}^{2}$ $\forall\sigma\in\llbracket s\rrbracket$. MSBPs with more general tree structured ground costs have appeared in [7].

When the cardinality of the index set $\llbracket s\rrbracket$ equals $2$, then (12) reduces to the (bi-marginal) Schrödinger bridge problem (SBP) [8, 9]. In this case, the solution of (12) gives the most likely evolution between two marginal snapshots $\mu_{1},\mu_{2}$. This can be established via the large deviations [10] interpretation [11, Sec. II] of SBP using Sanov’s theorem [12]; see also [13, Sec. 2.1].

Specifically, let $\mathcal{C}\left([\tau_{1},\tau_{2}],\mathbb{R}^{d}\right)$ denote the collection of continuous functions on the time interval $[\tau_{1},\tau_{2}]$ taking values in $\mathbb{R}^{d}$. Let $\Pi(\mu_{1},\mu_{2})$ be the collection of all path measures on $\mathcal{C}\left([\tau_{1},\tau_{2}],\mathbb{R}^{d}\right)$ with time $\tau_{1}$ marginal $\mu_{1}$, and time $\tau_{2}$ marginal $\mu_{2}$. Given a symmetric ground cost (e.g., Euclidean distance) $C:\mathcal{X}_{1}\times\mathcal{X}_{2}\mapsto\mathbb{R}_{\geq 0}$, let

and consider the bimarginal Gibbs kernel

Then, the bimarginal SBP solves

i.e., the most likely evolution of the path measure consistent with the observed measure-valued snapshots $\mu_{1},\mu_{2}$.

Under the stated assumptions on the ground cost $c$, the existence of minimizer for (16) is guaranteed [14, 15]. The uniqueness of minimizer follows from strict convexity of the map $\pi\mapsto D_{\rm{KL}}(\pi\|\nu)$ for fixed $\nu$.

This relative entropy reformulation, and thereby “the most likely evolution consistent with observed measures” interpretation, also holds for the MSBP (12) with $s\geq 2$ snapshots. Specifically, for $\bm{C}:\bm{\mathcal{X}}\mapsto\mathbb{R}_{\geq 0}$ as in (12)-(13), we generalize (14) as

and define the multimarginal Gibbs kernel

Problem (16) then generalizes to

where $\Pi(\mu_{1},\ldots,\mu_{s})$ denotes the collection of all path measures on $\mathcal{C}\left([\tau_{1},\tau_{s}],\mathbb{R}^{d}\right)$ with time $\tau_{\sigma}$ marginal $\mu_{\sigma}$ $\forall\sigma\in\llbracket s\rrbracket$. The equivalence between (12) and (19) can be verified by direct computation. Thus solving (19), or equivalently (12), yields the most likely evolution of the path measure consistent with the observed measure-valued snapshots $\mu_{\sigma}$ $\forall\sigma\in\llbracket s\rrbracket$.

We propose to solve the MSBP (12) for learning the time-varying statistics of the hardware resource state $\bm{\xi}$ as in (10). We next detail a discrete formulation to numerically solve the same for scattered data $\{\bm{\xi}^{i}(\tau_{\sigma})\}_{i=1}^{n}$ where $n$ is the number of control software execution profiles.

The minimizer of (12), $\bm{M}_{\rm{opt}}\left(\bm{\xi}(\tau_{1}),\ldots,\bm{\xi}(\tau_{s})\right)$ can be used to compute the optimal coupling between snapshot index pairs $(\sigma_{1},\sigma_{2})\in\{\llbracket s\rrbracket^{\otimes 2}\mid\sigma_{1}<\sigma_{2}\}$ as

where

This will be useful for predicting the statistics of $\bm{\xi}(\tau)\sim\mu_{\tau}$ at any (out-of-sample) query time $\tau\in[0,t]$.

With slight abuse of notations, we use the same symbol for the continuum and discrete version of a tensor. The ground cost in discrete formulation is represented by an order $s$ tensor $\bm{C}\in\left(\mathbb{R}^{n}\right)^{\otimes s}_{\geq 0}$, with components $\left[\bm{C}_{i_{1},\ldots,i_{s}}\right]=\bm{C}\left(\bm{\xi}_{i_{1}},\ldots,\bm{\xi}_{i_{s}}\right)$. The component $\left[\bm{C}_{i_{1},\ldots,i_{s}}\right]$ encodes the cost of transporting unit mass for a tuple $(i_{1},\ldots,i_{s})$.

Likewise, consider the discrete mass tensor $\bm{M}\in\left(\mathbb{R}^{n}\right)^{\otimes s}_{\geq 0}$ with components $\left[\bm{M}_{i_{1},\ldots,i_{s}}\right]=\bm{M}\left(\bm{\xi}_{i_{1}},\ldots,\bm{\xi}_{i_{s}}\right)$. The component $\left[\bm{M}_{i_{1},\ldots,i_{s}}\right]$ denotes the amount of transported mass for a tuple $(i_{1},\ldots,i_{s})$.

For any $\sigma\in\llbracket s\rrbracket$, the empirical marginals $\bm{\mu}_{\sigma}\in\mathbb{R}^{n}_{\geq 0}$ are supported on the finite sets $\{\bm{\xi}^{i}(\tau_{\sigma})\}_{i=1}^{n}$. We denote the projection of $\bm{M}\in\left(\mathbb{R}^{n}\right)^{\otimes s}_{\geq 0}$ on the $\sigma$th marginal as $\operatorname*{proj}_{\sigma}(\bm{M})$. Thus $\operatorname*{proj}_{\sigma}:\left(\mathbb{R}^{n}\right)^{\otimes s}_{\geq 0}\mapsto\mathbb{R}^{n}_{\geq 0}$, and is given componentwise as

Likewise, denote the projection of $\bm{M}\in\left(\mathbb{R}^{n}\right)^{\otimes s}_{\geq 0}$ on the $(\sigma_{1},\sigma_{2})$th marginal as $\operatorname*{proj}_{\sigma_{1},\sigma_{2}}(\bm{M})$, i.e., $\operatorname*{proj}_{\sigma_{1},\sigma_{2}}:\left(\mathbb{R}^{n}\right)^{\otimes s}_{\geq 0}\mapsto\mathbb{R}^{n\times n}_{\geq 0}$, and is given componentwise as

We note that (22) and (23) are the discrete versions of the integrals in (12b) and (20), respectively.

With the above notations in place, the discrete version of (12) becomes

The primal formulation (24) has $n^{s}$ decision variables, and is computationally intractable. Recall that even for the bimarginal ($s=2$) case, a standard approach [25] is to use Lagrange duality to notice that the optimal mass matrix $M_{\rm{opt}}$ is a diagonal scaling of $K:=\exp(-C/\varepsilon)\in\mathbb{R}^{n\times n}_{>0}$, i.e., $M_{\rm{opt}}=\mathrm{diag}(\bm{u}_{1})K\mathrm{diag}(\bm{u}_{2})$ where $\bm{u}_{1}:=\exp(\bm{\lambda}_{1}/\varepsilon)$, $\bm{u}_{2}:=\exp(\bm{\lambda}_{2}/\varepsilon)$, and $\bm{\lambda}_{1},\bm{\lambda}_{2}\in\mathbb{R}^{n}$ are the Lagrange multipliers associated with respective bimarginal constraints $\operatorname*{proj}_{1}(M)=\bm{\mu}_{1}$, $\operatorname*{proj}_{2}(M)=\bm{\mu}_{2}$. The unknowns $\bm{u}_{1},\bm{u}_{2}$ can be obtained by performing the Sinkhorn iterations

with guaranteed linear convergence [26] wherein the computational cost is governed by two matrix-vector multiplications.

The duality result holds for the multimarginal ($s\geq 2$) case. Specifically, the optimal mass tensor in (24) admits a structure $\bm{M}_{\rm{opt}}=\bm{K}\odot\bm{U}$ where $\bm{K}:=\exp(-\bm{C}/\varepsilon)\in\left(\mathbb{R}^{n}\right)^{\otimes s}_{>0}$, $\bm{U}:=\otimes_{\sigma=1}^{s}\bm{u}_{\sigma}\in\left(\mathbb{R}^{n}\right)^{\otimes s}_{>0}$, $\bm{u}_{\sigma}:=\exp(\bm{\lambda}_{\sigma}/\varepsilon)$, and $\bm{\lambda}_{\sigma}\in\mathbb{R}^{n}$ are the Lagrange multipliers associated with the respective multimarginal constraints (24b). The unknowns $\bm{u}_{\sigma}$ can, in principle, be obtained from the multimarginal Sinkhorn iterations [27]

which generalize (25). However, computing ${\rm{proj}}_{\sigma}\left(\bm{K}\odot\bm{U}\right)$ requires $\mathcal{O}\left(n^{s}\right)$ operations. Before describing how to avoid this exponential complexity (Sec. III-F), we point out the convergence guarantees for (26).

More recent works have also established [29] guaranteed convergence for the continuous formulation (12) with linear rate of convergence [30].

We next outline how the solution tensor $\bm{M}_{\rm{opt}}=\bm{K}\odot\bm{U}$ obtained from the converged Sinkhorn iterations can be used together with (28), to make stochastic predictions of the most likely hardware resource state in the form (10).