In the real-time systems community, it is well-known [1] that due to hardware-level stochasticity, multiple executions of the same software on the same hardware with identical initial conditions and parameters, result in different execution times. A common way to account for this is to analyze some lumped variable such as (worst case or probabilistic) execution time [2], [3], [4] for a given software-hardware combination. In contrast, the ability to directly learn in the joint space of LLC, processor and memory availability, could enable the design of more fine-grained dynamic resource schedulers. This is of particular interest in safety-critical controlled cyber-physical systems that operate in resource-constrained environments because high-performance controllers [e.g., model predictive controller (MPC)] are more compute-intensive and have more pronounced stochasticity than computationally benign controllers [e.g., proportional-integral-differential (PID)].

Learning in the joint space, however, is technically challenging because variables such as processor availability, LLC, and memory bandwidth are not only statistically correlated, but their correlations also change over time, i.e., they cannot be assumed as (neither strict nor wide sense) stationary stochastic processes for learning purposes. It is also impractical to fit parametric statistical models such as a Markov chain by gridding the multidimensional space of computational resource state since the resource state space is not countable. In this work, we leverage the recent progress in multimarginal Schrödinger bridges to show the feasibility of nonparametric joint stochastic learning of computational resource usage from hardware-software profile data.

Related Works: Multimarginal Schrödinger bridge problems (MSBPs) are entropy regularized variants of the multimarginal optimal mass transport (MOMT) problems [5], [6]. The latter has been applied to learning and inference problems in chemical physics [7], [8], team matching [9], fluid dynamics [10], and risk management [11]. For recent applications of MSBPs to learning and tracking problems, see [12], [13], [14]. The MSBPs can also be seen as generalized variants of the classical (also known as bimarginal) Schrödinger bridge problems (SBPs) which we will explain in Section IV-A. In the control literature, there are growing works on the stochastic optimal control interpretations [15], [16], generalizations [17], [18], [19], and applications [20], [21], [22] of the SBPs.

Prior work has applied machine learning to performance prediction and workload modeling, especially in data center and cloud environments; e.g., the surveys in [23] and [24] and references therein. Unlike our work, these solutions focus on predicting coarse-grained characteristics-such as request rates, CPU utilization, memory and disk I/O usages, bandwidth, or energy consumption--and they do not consider the interdependent relationship among shared resources. In realtime and embedded systems, learning techniques have been used for estimating timing behaviors [25], [26], [27]; however, existing work focuses on the worst case or probabilistic execution time of the entire software, rather than the dynamic resource usage patterns during the execution of the software. To the best of the authors' knowledge, our work is the first to use learning for predicting the stochastic dynamical runtime behavior of computational resource usage that jointly considers interdependent resources on multicore platforms (such as CPU, cache, and memory bandwidth) as well as multithreaded applications.

Contributions: This work builds on our preliminary work [28] but significantly extends the same by considering multicore resources and more general graph structures than path trees. The specific contributions are threefold.

1) Mathematical formulations (Sections III and IV) for the stochastic learning of computational resource usage as graph-structured MSBPs where the graph structures arise from single-core or multicore computational platforms. These formulations are motivated by maximum likelihood interpretations in the path space. 2) Numerical algorithms (Section V) for solving the proposed graph-structured MSBPs. We show that unlike generic graph structures, the MSBPs induced by singlecore or multicore computational resource usage lead to structured tensor optimization problems that can be leveraged to reduce the learning algorithms' computational complexity from exponential to linear in number of observations. 3) Numerical experiments (Section VI) to illustrate the proposed learning framework. This includes a singlecore MPC benchmark and a synthetic multicore software benchmark.

In this section, we fix notations and background ideas that will be used subsequently.

For any natural number ν, we use the finite set notation [[ν]] := {1, 2, . . . , ν}. We denote the cardinality of set S as |S|. We will use the following.

Lemma 1 [29,Th. 19.2]: Given finite sets S 1 ⊆ S 2 with

We use unboldfaced (resp. boldfaced) small letters to denote scalars (resp. vectors). Unboldfaced capital letters denote matrices and bold capital letters denote tensors of order three or more. We occasionally make exceptions for standard notations such as (2) and (3).

We use square braces to denote the components. For example, [X i 1 ,...,i r ] denotes the (i 1 , . . . , i r )th component of an order r tensor X, where (i 1 , . . . , i r ) ∈ N r . We use the r fold tensor product space notation (R d ) ⊗r := R d ⊗ . . . ⊗ R d " ƒ‚ … r times . For two given tensors X and Y of order r, their Hilbert-Schmidt inner product is

The operators exp(•), log(•), , and are all understood elementwise, i.e., denote the elementwise exponential, logarithm, multiplication, and division, respectively. Vector or matrix transposition is denoted by the superscript . We use diag(v) to denote a diagonal matrix with entries of vector v along its main diagonal, and 1 to denote all one's column vector of suitable length.

A probability measure µ over some polish space X satisfies R X dµ = 1. For a pair of probability measures µ and ν defined over two polish spaces X and Y, respectively, their product measure is µ ⊗ ν, and

where (dµ/dν) denotes the Radon-Nikodym derivative, and µ ν is a shorthand for "µ is absolutely continuous with respect to ν."

The Wasserstein distance W between two probability measures µ and ν, supported, respectively, on X , Y ⊆ R d , is

xy 2 2 dπ(x, y)

where the infimum in ( 3) is over all joint couplings of µ and ν, i.e., Π(µ, ν)

Unlike the Kullback-Leibler divergence (2), the Wasserstein distance (3) is a metric on the space of probability measures.

1) Hilbert Projective Metric: The Hilbert projective metric [30], [31], [32]

The convergence plots with respect to Hilbert metric for our numerical experiments will be given in Section VI. Formula (4) is a special case of the general definition of Hilbert metric d H between any two elements of a pointedfoot_0 convex cone K in a real vector space

Note that (5) reduces to (4) for K ≡ R n >0 .

We collect distributional data for the computational resource availability state thought of as a continuous-time stochastic process over time τ ∈ [0, t]. In what follows, we will make these ideas precise. For now, let us begin by assuming that the data are collected at s ∈ N, s ≥ 2 snapshots

Define the snapshot index set [[s]] := {1, 2, . . . , s}.

We consider multicore computing hardware with J ∈ N cores and define the core index set [[J]] := {1, 2, . . . , J}. For j ∈ [[J]] and for a given time horizon of interest [0, t], we think of the computational resource state vector ξ j (τ), 0 ≤ τ ≤ t, as a continuous time R d -valued stochastic process.

For each j ∈ [[J]], a (random) vectorial sample path ξ j (τ) is referred to as a profile. For the single-core (J = 1) case, we drop the core index superscript j and simply denote the sample path or profile as ξ(τ).

As a concrete example, consider the case

with the components of ξ j as follows:

In this example, the three elements of ξ j denote the number of CPU instructions, the number of LLC requests, and the number of LLC misses in the last time unit (e.g., in the last 10 ms if the profiling frequency is 100 Hz), respectively, in the core j ∈ [[J]]. We will use this specific ξ j ∈ R 3 in our numerical experiments (both single core and multicore). However, we will develop the proposed stochastic learning framework for generic ξ j ∈ R d . This will allow generalizability in the sense that the proposed method can be adapted to an application at hand with a custom definition of the d ∈ N components of the stochastic state ξ j . The distributional (i.e., measure valued) observations collected at instances (6) comprise a sequence of joint state probability distributions or measures

In this work, we are interested in learning the hardware-level stochasticity, i.e., the stochasticity in the state ξ j arises from dynamic resource variability. In other words, repeated execution of the same (e.g., control) software with the same initial condition and same parameters result in different profiles ξ j (τ), and 9 as exemplar single-core and multicore profiles, respectively. Intuitively, when the hardware-level stochasticity is negligible, then the distributions µ j σ will be approximately Dirac deltas supported on the graph of a single path. When the stochastic variability is significant, µ j σ will have significant dispersion.

In practice, the probability measures {µ j σ } ( j,σ)∈[[J]]× [[s]] are only available empirically from a fixed, say n ∈ N profiles. Let the sample or profile index i ∈ [[n]]. We accordingly set

where δ(ξ jξ i, j (τ σ )) denotes the Dirac delta at the ith sample location ξ i, j (τ σ ) ∈ R d for a fixed index pair ( j, σ)

For any fixed pair ( j, σ)

With the basic notations in place, the informal statement for our stochastic learning problem is as follows.

Most Likely Distributional Learning Problem (Informal): Given distributional snapshots (8), predict the most likely distribution of the computational resource state

In Section IV, we will formalize this problem statement and discuss the related information graph structures.

To motivate the mathematical formulation, we start by outlining the classical (bimarginal) SBP in Section IV-A. Then, in Section IV-B, we make the informal statement mentioned at the end of the previous section rigorous using the framework of large deviation principle [33]. The resulting MSBP formulation takes the form of a maximum likelihood problem in the space of probability measure-valued curves, generalizing a similar formulation for the classical SBP.

We then consider specific cases of this MSBP formulation which result from the information graph structures induced by the stochastic profiles in single-core and multicore computing hardware. In particular, Section IV-C details MSBP over a path tree which is the information graph structure that arises in the single-core profiling. For the multicore case, we discuss two different information graph structures in Sections IV-D and IV-E.

In 1931 and 1932, Schrödinger [34], [35] formulated which is now called the classical SBP, in two works: one written in German and another in French. For a relatively recent survey, see [36]. The classical SBP corresponds to the bimarginal (i.e., s = 2) maximum likelihood problem: its solution finds the most likely measure-valued curve µ τ , where τ ∈ [0, t] connecting the given endpoint measures µ 1 and µ 2 at τ 1 = 0 and τ 2 = t with respective supports over subsets of R d , see Fig. 1(a).

The classical SBP formulation proceeds as follows. Let X 1 , X 2 ⊆ R d be the supports of the given endpoint measures µ 1 and µ 2 , respectively. Letting X :=

) is a distance between the random vectors ξ(τ 1 ), ξ(τ 2 ). Let M(X ) be the manifold of joint probability measures on the product space X . For fixed but not necessarily small

In words, (11) seeks to compute a joint probability measure that minimizes the entropy-regularized transportation cost (11a) subject to prescribed marginal constraints (11b) and (11c). Note that M(X ) is a convex set. The objective (11a) is strictly convex in M, thanks to the ε-regularized negative entropy term R X ε log M dM. Constraints (11b) and (11c) are linear.

To clarify the maximum likelihood interpretation of (11

and consider the bimarginal Gibbs kernel

Recall that the minimizer of problem ( 11) is an optimal coupling M opt ∈ M(X ). Proposition 1 next formalizes why the solution of ( 11) corresponds to the most likely measure-valued path [Fig. 1(a)] consistent with the observed measure-valued snapshots µ 1 and µ 2 . Its proof uses Sanov's theorem [37]; for details, we refer the readers to the references below.

The optimal coupling M opt ∈ M(X ) in ( 11) corresponds to the optimal measure-valued path π opt ∈ Π(µ 1 , µ 2 ) solving the (scaled) relative entropy minimization problem min

Remark 1 (Existence Uniqueness of the Solution for ( 14)): Under the stated assumptions on the ground cost C, the existence of minimizer for ( 14) is guaranteed [40], [41]. The uniqueness of the minimizer follows from strict convexity of the map π → D KL (π ν) for fixed ν.

Remark 2: Intuitively, Proposition 1 links the solution of a static optimization problem (11) with that of a dynamic optimization problem (14) in the sense that the minimization in ( 14) is performed over all continuous measure-valued curves connecting the endpoints µ 1 and µ 2 .

Our setting in (10) requires generalizing these ideas for the multimarginal (s ≥ 2) case discussed next.

In the single-core (J = 1) case, Fig. 1(b) illustrates how the MSBP of our interest generalizes the classical, also known as bimarginal SBP, as shown in Fig. 1(a).

In general, for any J ∈ N, we start by defining

the Cartesian product

and a ground cost C : X → R ≥0 . Let M(X ) be the manifold of probability measures on X , and for a fixed pair ( j, σ)

For a fixed ε > 0, the MSBP, generalizes the bimarginal problem (11) as follows:

The maximum likelihood interpretation for the bimarginal problem given in Proposition 1 generalizes for MSBP (18) as stated next in Proposition 2. Similar to the case in Proposition 1, its proof is by direct computation and omitted. Proposition 2: Given the sets ( 15) and ( 16), ground cost C : X → R ≥0 and fixed ε > 0, define

and the multimarginal Gibbs kernel

Then, the optimal coupling M opt ∈ M(X ) in ( 18) corresponds to the optimal measure-valued path π opt ∈ Π(µ 1 1 , . . . , µ J s ) solving the (scaled) relative entropy minimization problem min

The existence uniqueness of the minimizer for (21), and thus for (18), follows from the same strict convexity argument as in Remark 1.

Motivated by the maximum likelihood interpretation ( 21) for (18), we propose to solve (18) for learning the stochastic computational resource state (10). The minimizer of ( 18), M opt (ξ 1 (τ 1 ), . . . , ξ J (τ s )) can be used to compute the optimal bimarginal probability mass transport plans between any ( j 1 , σ 1 ), (

where

The coupling (22) will find use in implementing our multimarginal Sinkhorn recursions in Section V.

We next discuss the discrete formulations of ( 18) that arise from the information graph structures in single-core and multicore computing hardware for finite scattered data

] as in (9). We will see that the corresponding formulations lead to strictly convex problems over finite dimensional nonnegative tensors. To reduce the notational overload, we use the same boldfaced symbols for the continuous and discrete versions of the tensors.

For single-core (J = 1) computing hardware, the stochastic process ξ(τ) is indexed only by time τ ∈ [0, t]. Therefore, the information graph structure is also induced by time, i.e., by sequential in time measure-valued observations {µ σ } σ∈ [[s]] . In other words, the information graph in this case is a path tree shown in Fig. 2. We refer to an MSBP specialized to such path tree as the path-structured MSBP.

Thanks to the path tree structure, the ground cost C in this case can be written as follows:

where we choose the squared Euclidean distance sequential cost between two consecutive snapshot indices, i.e., c σ (•,

To formulate the discrete version of the corresponding MSBP, note that the ground cost in (18a), in general, is an order

that encodes the cost of transporting a unit amount of mass for an s tuple (i 1 , . . . , i s ). However, the path-structured cost (23) implies that the s tuple in ( 24) equal to P s-1 σ=1 c σ (ξ i σ , ξ i σ+1 ), which is a sum of suitable elements of s -1 different (in our case, Euclidean) distance matrices.

The discrete mass tensor M ∈ (R n ) ⊗s ≥0 has components M i 1 ,...,i s = M ξ i 1 , . . . , ξ i s (25) where [M i 1 ,...,i s ] denotes the amount of transported mass for an s tuple (i 1 , . . . , i s ). Furthermore, for path-structured MSBP, we can drop the core index superscript j to simplify constraints (18b) as follows: Z

In the discrete version, the empirical marginals

Then, the LHS of constraints (26) are projections of the mass tensor (25) on the σth marginal µ σ ∀σ ∈ [[s]]. We denote this projection as proj σ (M), which is a mapping proj σ : (R n ) ⊗s ≥0 → R n ≥0 , and is given componentwise as follows:

Similarly, the discrete version of ( 22) in the path-structured case, is the projection of M ∈ (R n ) ⊗s ≥0 on the (σ 1 , σ 2 )th marginal denoted as proj σ 1 ,σ 2 (M), i.e., proj

≥0 and is given componentwise as follows: Then, the discrete path-structured version of ( 18) becomes min

Note that ( 29) is a strictly convex problem in n s decision variables and is computationally intractable in general. In Section V, we will design algorithms to solve them efficiently.

For multicore (J > 1) computing hardware, the stochastic process ξ j (τ) is indexed by both time τ ∈ [0, t] and CPU core j ∈ [[J]]. Thus, the information graph structure is induced jointly by time and cores, i.e., by measure-valued observations

. Unlike Section IV-C, we now have J measurevalued snapshots at each time index. To account for this, we propose two different MSBP formulations, the first of which termed BC MSBP is discussed next.

The main idea behind our proposed BC MSBP formulation is to imagine a phantom CPU core whose statistics is the average (i.e., barycenter) of the measure-valued snapshots for all cores at any fixed time index σ ∈ [[s]]. Accordingly, in this formulation, we consider the core index j ∈ {0} ∪ [[J]] where j = 0 refers to the phantom BC CPU. The corresponding graph structure for J = 2 is shown in Fig. 3. While a BC graph as shown in Fig. 3 is more general than the path tree as shown in Fig. 2, we note that the treewidth [42, pp. 354 and 355] for BC graphs equals unity as in the path tree case. It is known that the computational complexity for graph-structured MSBP problems grows with the treewidth [43], [44]. We will discuss the specific complexity of the proposed algorithms in Section V.

The BC graph structure implies that ( 24) is no longer equal to P s-1 σ=1 c σ (ξ i σ , ξ i σ+1 ) as shown in Section IV-C, but instead becomes a sum of two types of ground transport costs, viz., the cost of transport between consecutive-in-time barycenters and the cost of transport between BC and actual CPU cores. Denoting the BC random vectors as {ξ 0 (τ σ )} σ∈[[s]] , and letting c j,σ (•, •) be the corresponding ground transport costs for all j ∈ {0} × [[J]], the cost (24) for BC MSBP equals

Similar ideas appeared in [12, Sec. 3.3] in a different context. Defining the BC index set

we note that |Λ BC | = (J + 1)s, and thus for the BC MSBP, the tensors

. In summary, the proposed BC MSBP is min

In Section V-B, we will discuss the computation of projections (18b) and ( 22) for the minimizer of ( 31).

Different from the BC MSBP in Section IV-D, we now propose another MSBP formulation for the multicore (J > 1) case, based on the SP information graph structure in Fig. 4.

This SP graph-structured MSBP is motivated by the observation that many software on multicore computing hardware have notions of input and output terminals. In such applications, the input and output of the computational task are aggregated to a single core. Thus, the information graph structure is induced by measure-valued observations {µ j σ } ( j,σ) , where the tuple ( j, σ) belongs to the index set

Although the SP graph in Fig. 4 is not a tree, it has treewidth at most two. Therefore, such MSBPs remain computationally tractable as before (see details in Section V). Unlike the formulation in Section IV-D, the SP graphstructured MSBP formulation does not involve any phantom CPU core. From Fig. 4, the cost (24) now equals

By Lemma 1, we find |Λ SP | = J(s-2)+2, and hence for the SP graph-structured MSBP, the tensors C, M ∈ (R n ) ⊗J(s-2)+2 ≥0 . In summary, the proposed SP graph-structured MSBP is min

In Section V-C, we will discuss the computation of projections (18b) and ( 22) for the minimizer of (33).

This section provides algorithmic details and computational complexities to solve the discrete MSBPs ( 29), (31), and (33).

These MSBPs are strictly convex tensor optimization problems in n s , n (J+1)s , and n J(s-2)+2 decision variables, respectively, and computationally intractable in general. By leveraging an interplay between duality and graph structures, we will see that it is possible to reduce the computational complexity from exponential to linear in s.

Recognizing that ( 29), (31), and ( 33) are instances of the generic-structured optimization problem min

for suitable index set Λ (see Table I), the following is a consequence of strong Lagrange duality. Proposition 3 (See [45], [46], [47]): Given problem (34), let λ j σ ∈ R n be the Lagrange multipliers for the equality constraints (34b) for all ( j, σ) ∈ Λ. Let

Then, the multimarginal Sinkhorn recursions

have guaranteed linear rate of convergence, 2 and the minimizer M opt for (34) is given in terms of the converged U as follows:

Once M opt is computed, its bimarginal projections of the form (28) yield the probability mass transport matrices between any two marginals. With this, we are able to interpolate between any two marginals to obtain a predicted distribution. 2 The proof sketch for discrete state space is as follows. The MSBP (34) can be recast [45,Sec. 4.1] as a Kullback-Leibler projection to a convex set that is an intersection of |Λ| hyperplanes given by (34b). That this iterative Kullback-Leibler projection has guaranteed convergence follows from the seminal result on iterative Bregman projection [48] and from the fact that the Kullback-Leibler divergence is an instance of Bregman divergence. That the rate is linear follows from the coordinate descent analysis by Luo Remark 3: In our problem, the useful application of this interpolation is as follows: given a time τ ∈ [0, t) and CPU index j ∈ [[J]], find σ ∈ [[s]] such that τ σ ≤ τ < τ σ+1 and use the expressions for the bimarginal projections above to compute the bimarginal bridge M j,σ := proj ( j,σ),( j,σ+1) (M opt ) :

Using this projection, we can interpolate between µ j σ and µ j σ+1 to obtain our estimate b µ j τ for the computational resource usage distribution for CPU j at time τ, as follows:

where b ξ j (τ, ξ r, j (τ σ ), ξ ℓ, j (τ σ+1 )) := (1λ)ξ r, j (τ σ ) + λξ ℓ, j (τ σ+1 ), and

For numerically solving bimarginal SBPs (i.e., the case s = 2 and J = 1), the Sinknorn recursions (36) have become the standard [50], [51], [52], [53]. However, applying the same is challenging for MSBPs because computing proj ( j,σ) (K U) requires O(n |Λ| ) operations. The same issue arises for computing proj ( j 1 ,σ 1 ),( j 2 ,σ 2 ) (K U). For both projections, this complexity can be reduced by exploiting the structure of the Hilbert-Schmidt inner product K, U . Specifically, we make use of the following results from [12].

Lemma 2: Let the tensor U be as in (35c), and consider a tensor K ∈ (R n ) ⊗|Λ| .

1) [12, Lemma 1]: For a fixed ( j, σ) ∈ Λ, if K, U = w 1 diag(u j σ )w 2 for some vectors w 1 , w 2 ∈ R n that do not depend on u j σ , then

2) [12, Lemma 2]: For fixed ( j 1 , σ 1 ), ( j 2 , σ 2 ) ∈ Λ, if K, U = w 1 diag(u j 1 σ 1 )Φdiag(u j 2 σ 2 )w 3 for some vectors w 1 , w 3 ∈ R n and matrix Φ ∈ R n×n , where w 1 , Φ, and w 3 do not depend on u j 1 σ 1 and u j 2 σ 2 , then

The special inner product structures in Lemma 2 arise from structured tensors K which are in turn induced by the graph structures discussed earlier. Thus, expressing K, U in the appropriate forms helps compute the desired projections. Below, we show how the imposition of structure on K allows this to be done efficiently.

For software running on a single-CPU core, we have J = 1 as in Section IV-C. The corresponding MSBP is (29). Here, the cost tensor C in (24) has a path structure

where the matrix C σ ∈ R n×n ≥0 encodes the cost of transporting unit mass from {ξ i (τ σ )} n i=1 to {ξ i (τ σ+1 )} n i=1 . This allows us to write K in (35a) as follows:

which leads to the following expressions for the marginal projections. Proposition 4 (See [12, Proposition 2] and [28, Proposition 1]): If C has the form (39), K σ := exp(-C σ /ε) ∈ R n×n ≥0 , K as in (40), and U as in (35c); then, ( 27) and ( 28) can be expressed as follows:

and proj

Observe that even the Naïve computation of ( 41) is dominated by 2s -4 matrix-vector multiplications (by cancellation, (36) can be computed by s -1 such multiplications, see [28,Remark 3]). Since such multiplications have O(n 2 ) complexity, each Sinkhorn iteration has O((s-1)n 2 ) complexity-linear in s, and a great improvement from the general O(n s ) complexity of the method. Our method's O((s -1)n 2 ) complexity is sharp for exact computation. The recent work [54] reduces the complexity in n from quadratic to linear at the expense of approximate computation.

Remark 4: While it is clear from their expressions that the bimarginal projection (42) has similar order-of-magnitude complexity to the unimarginal projection (41), hereafter we focus only on the complexity of the latter, as these are performed every Sinkhorn iteration until the method converges. Following this, bimarginal projections of the solution M opt onto each pair of marginals of interest are to be performed only once a posteriori, and so exact floating-point operational count for these projections is not critical.

For software running on multiple CPU cores, we have J > 1. The corresponding BC MSBP (31) formulated in Section IV-D involves a tree that is more general than a path.

In this section, we let n 0 denote the number of samples in the BC CPU. For the sake of generality, here we derive the complexity for the projections in terms of n and n 0 . Then,

For n 0 n, formulation ( 34) and ( 35) applies mutatis mutandis.

Recall that here the index set Λ = Λ BC as in (30). The cost tensor C for BC MSBP takes the form

where the matrices C 0,σ ∈ R n 0 ×n 0 ≥0 are the ground cost matrices between barycenters ξ 0 (τ σ ) and ξ 0 (τ σ+1 ) for σ ∈ [[s -1]], whereas C j,σ ∈ R n 0 ×n ≥0 are those between the CPU marginals ξ j (τ σ ) and their barycenters ξ 0 (τ σ ), for ( j, σ)

. So now, we can write K in (35a) as follows:

With this, the projections can be computed as follows (proof in Appendix A).

Proposition 5: If C has the form ( 43), K j,σ := exp(-C j,σ /ε), K as in (44), and U as in (35c), then the projections ( 27) and ( 28) can be expressed as follows:

Authorized licensed use limited to: Iowa State University. Downloaded on December 03,2025 at 09:28:11 UTC from IEEE Xplore. Restrictions apply.

where

Similar to the single-core case in Section V-A, the computational complexity for all projections in Proposition 5 are linear in J and s. Specifically, the computation of the s p σ vectors in Proposition 5 requires a total of Js matrixvector multiplications, and the total number of floating-point operations when projecting onto a barycenter [as shown in (45a)] is

The number of floating-point operations when projecting onto a non-BC marginal [as shown in (45b)] is

The SP graph-structured MSBP (33) formulated in Section IV-E involves an SP graph which is not a tree. Recall that the corresponding index set Λ = Λ SP as in (32). For any fixed j ∈ [[J]], let

Then, the cost tensor C for SP graph-structured MSBP takes the form

where each C j ∈ (R n ) ⊗s >0 is the path-structured cost tensor along the jth CPU's path, i.e., h

wherein the matrices C j,1 ∈ R n×n ≥0 are the ground cost matrices between marginals ξ 1 (τ 1 ) and ξ j (τ 2 ) for j ∈ [[J]], and similarly C j,s-1 maps between ξ j (τ s-1 ) and ξ 1 (τ s ). When σ ∈ [[s-2]]\{1}, the matrices C j,σ map between ξ j (τ σ ) and ξ j (τ σ+1 ).

Consequently, we write K in (35a) as follows:

and the associated projections can be expressed as follows (proof in Appendix B). Proposition 6: If C has the form (49), K j,σ := exp(-C j,σ /ε), K as in (51), and U as in (35c), then letting

and

, projection (27) takes the form

where ( j, σ) ∈ S k∈[[J]] Λ k SP in (53c). Furthermore, projection (28) can be expressed as follows:

Directly computing A k in (52) requires 1 + 2(s -2) matrixmatrix multiplications, each in this case involving n 3 floatingpoint operations. Thus, (53a) and (53b) each involve

operations, and (53c) involves

operations.

Remark 5: For n 0 = n, s ≥ 2, J ≥ 2, we note that the expression (48) = (Js + 2s)n 2 + (Js + 4)n ≥ 8n 2 + 8n.

In contrast, expression (56)= (2Js -3J + 3)nfoot_2 + (J -1)n 2 + n ≥ n 3 + n 2 + n. Thereby for n large, Propositions 5 and 6 enable computation of unimarginal projections with O(Jsn 2 ) complexity for the BC case, and O(Jsn 3 ) for the SP case. The n 3 complexity of the latter is a consequence of the treewidth (=2) of the SP graph structure. Both cases have complexity linear in Js.

Having shown that the computation of the projections is tractable in the number of marginals, we are able to efficiently solve (29) for all of our proposed graph structures: path, BC, and SP. In applying MSBP for computational resource prediction, the choice of problem graph structure can affect the quality of the prediction. For single-core software, all graph structures necessarily degenerate to the path structure. But for multicore software, there are two structurally distinct options: BC and SP. In Section VI, we evaluate the merits of these structures.

Regardless of the choice of graph structure, our overall methodology is as follows.

Step 1: Execute the software of interest n times over [0, t], generating hardware resource state snapshots

Our marginal distributions µ j σ are then as per (9).

Step 2: Pick an appropriate graph structure (path, BC, or SP). From the marginals in Step 1, construct the bimarginal Euclidean distance matrices C j,σ . If the BC structure is used, compute the barycenters µ 0 σ . Run the Sinkhorn recursion (36) until convergence, i.e., until the Hilbert projective metric (4) between the new and old u j σ is less than some numerical tolerance for all marginals. The resulting M opt given by (37) is the solution for our MSBP.

Step 3: Given any time τ ∈ [0, t) and CPU index j, use (38) to predict our estimate b µ j τ for the hardware resource usage distribution of CPU j at time τ.

In Sections IV and V, we explained how to formulate and efficiently solve graph-structured MSBPs to compute the optimal mass transport plan M opt between marginal distributions (9) which in our application correspond to the statistics of computational resource usage at times (6).

In Remark 3, we outlined how the computed M opt can be used to make predictions on a software's resource usage at any (possibly out of sample) given time. To evaluate the quality of these predictions, in this section, we perform two experiments: one with a single-core software and one with a multicore software. We quantify the quality of our predictions against known resource usage data.

For all experiments herein, benchmark programs were profiled on an Intel Xeon E5-2683 v4 processor with 16 cores and a 40-MB LLC running Ubuntu 16.04.7. We leveraged Intel's Cache Allocation Technology (CAT) [55] and Memguard [56] to control the amount of LLC and memory bandwidth available to the benchmark programs.

The implementation of the algorithms described in Section V, the solution of the MSBP (34) itself, and all other processing of profiled data, were done in MATLAB R2019b on a Debian 12 Linux machine with an AMD Ryzen 7 5800× CPU.

For this experiment, we used a custom single-core software 3 which we wrote in C language, implementing a path-following nonlinear model predictive controller (NMPC) for a kinematic bicycle model (KBM) [21], [57]. At each control step, the NMPC used the IPOPT nonlinear program solver [58] to determine the control decision which minimizes the sum of various measurements of the deviation from the desired path while promoting smoothness of the control inputs. For details on this control software, we refer the readers to [28, Sec. V].

1) Profiling: To gather the execution profiles for our NMPC control software, we ran our application on an isolated CPU and used the Linux perf tool [59], version 4.9.3, to sample the computational resource state vector ξ as in (7) every 10 ms.

For each profile, we ran the NMPC software for n c := 5 control steps with a fixed reference trajectory and a fixed combination of LLC and memory bandwidth partitions (30 MB and 30 MB/s, respectively, both allocated in blocks of 2 MB) to ensure a fixed execution context.

For our chosen context, we profiled the controller n := 500 times. Fig. 5 overlays all of the profiles, split into its three components as per (7).

2) Numerical Results: To account for the asynchrony among the profiles, we analyzed the statistics of the control cycle end times, i.e., the (random) times needed to complete one pass of the NMPC feedback loop. For details, we refer the readers to [28, Fig. 4 and Table I]. We then consider snapshots of the profiling data at times corresponding to the average control cycle end times as well as at s int equispaced intervals between control cycles. This results in s = 1 + n c (s int + 1) total marginals forming our path-structured MSBP (29) and completes Step 1 in Section V-D.

Setting the regularizer ε = 0.1 and a numerical tolerance of 10 -14 , we solve the discrete path-structured MSBP (29) with squared Euclidean cost C as in (39). Fig. 6 depicts the linear convergence of the Sinkhorn iterations (as in Proposition 3) shown with respect to the Hilbert metric (4). The path structure of the problem yields efficient computation of the Sinkhorn iterations leading to rapid convergence. In particular, setting s int := 4 (and so s = 26), we solve the MSBP (29) with n s = 500 26 decision variables in approx. 10 s and ≈120 Sinkhorn iterations with no optimizations to the code. Note that problems of this size are impractical to load in memory, let alone to solve efficiently using off-the-shelf solvers. This completes Step 2 in Section V-D.

For TABLE II FOR THE SINGLE-CORE EXPERIMENT IN SECTION VI-A, THE NUMBER OF INTRACYCLE MARGINALS s int VERSUS WASSERSTEIN DISTANCES W b σ AS IN (57) duration of the third control cycle, i.e., between τ 2(s int +1)+1 and τ 3(s int +1)+1 , with b τ b σ = τ 2(s int +1)+1 + Â τ 3(s int +1)+1τ 2(s int +1)+1 s int + 2 Ã b σ where the index b σ ∈ [[s int + 1]]. Since b τ b σ ∈ [τ 2(s int +1)+b σ , τ 2(s int +1)+b σ+1 ], we used (38) with σ = 2(s int + 1) + b σ to arrive at the predictions.

In Fig. 7, we compare our predictions against the observed empirical distributions. Fig. 7 shows that our predictions capture the mode(s) of the measured distributions quite well. Further improvements can be made by placing marginals closer together in time. In this experiment, we do this by increasing the value of s int , i.e., by placing more marginals equispaced between control cycle boundaries. In Table II, we report the Wasserstein distances W(•, •) as in (3) between the corresponding predicted and measured distributions We computed (57) as the square roots of the optimal values of the corresponding Kantorovich linear programs [60, Ch. 3.1] that results from specializing (29) with s = 2 and ε = 0. As expected, placing the marginals closer together in time results in higher accuracy in predictions.

For this experiment, we used the Canneal benchmark from the PARSEC suite [61]. Canneal is a resource heavy, multithreaded application that simulates an anneal workload to minimize routing costs for chip design. Specifically, Canneal pseudorandomly picks pairs of input elements to swap in a tight loop, leading to a heavy reliance on both cache and memory bandwidth. This can be seen from Fig. 8, where the average number of LLC misses decreases as the cache size increases. The dashed vertical lines represent the cache allocation used in our profiling.

1) Profiling: To enable multicore profiling, we made a small modification to Canneal's source code to pin each created thread to its own core. We then use perf, as with the single-core experiment, to sample the computational resource state ξ j , j ∈ [[J]], as in (7) for every 10 ms.

Desiring to examine how Canneal behaves when its resource allocation varies across its cores, we profile using a single context defined by J = 4 cores, and resource allocations of (24, 10, 4, 2) MB/s of LLC and (125, 25, 5, 1) MB/s of memory bandwidth, ordered left-to-right by increasing CPU number. Therefore, CPU 1 has the highest resource allocation while CPU 4 has the lowest. For this choice of context, we collected a total of n = 400 profiles for the Canneal software; these profiles are shown in Fig. 9. For the BC formulation, we fix n 0 = 600.

2) Numerical Results: We placed snapshots at times τ σ ∈ {0.0, 0.5, 1.5, 2.5, 5.0, 9.5, 10.5} (so s = 7) and interpolated at times b τ b σ ∈ {0.8, 2.2, 7.0, 9.0, 10.0} (so b σ ∈ [[5]]). Both BC and SP methods were then used to solve the MSBPs (31) and (33), respectively, therefrom we used (38) to estimate the resource usage distributions at the desired times.

Once again, the code for solution of ( 31) and ( 33), including the computation of the projections in Propositions 5 and 6, were implemented in MATLAB with some minor optimizations due to the large size of the problem (the BC and SP formulations have 400 35 and 400 22 decision variables, respectively; see Section V). TABLE III WASSERSTEIN DISTANCES (58) BETWEEN THE MEASURED DISTRIBU-TIONS AND THOSE PREDICTED IN THE BC CASE (TOP) AND SP CASE (BOTTOM) FOR THE CANNEAL BENCHMARK IN SECTION VI-B AT EACH τb σ∈[[5]]

With the regularizer ε = 0.05 and a numerical tolerance of 10 -13 , the BC and SP algorithms converged in 0.4810 and 0.5055 s, with 110 and 127 iterations, respectively (see Fig. 10). The predicted distributions, marginalized to each component of ξ j , j ∈ [[J]], are shown in Fig. 11. Finally, Table III lists the comparison of the Wasserstein distances of the predicted distributions for each CPU from those measured from our profiles, i.e.,

These figures show the same behavior as in the singlecore experiment-primarily, that our predictions accurately capture the mode(s) of the resource usage distribution at all interpolation times. Note that at b σ = 1 and b τ b σ = 0.8, all CPUs except the CPU 1 are idle. Our predictions match the measured distribution exactly for the CPUs 2-4 (i.e., the Dirac delta at 0 for all components), but not so for CPU 1, which has a nontrivial distribution.

This work explores a new vision for learning and predicting stochastic time-varying computational resource usage from hardware-software profile data in both single-core and multicore platforms. We propose to formulate the problem as a distributional learning problem directly in the joint space of correlated computational resources such as processor availability and memory bandwidth. This leads to graph-structured MSBPs, where the specific graph structures are induced by underlying singl-core or multicore nature of computation. At first glance, such formulations for scattered profile data appear computationally intractable because they involve tensorial convex optimization problems where the number of primal variables is exponential in the number of observational snapshots. By leveraging strong duality and graph structures, we show that the computational complexities can be reduced from exponential to linear without approximation, and these problems can in fact be solved efficiently via nonparametric computation with maximum likelihood guarantees in the space of distributions. This enables us to predict the most likely joint (and hence all marginal) computational resource usage distribution at any user-specified query time. We emphasize that our proposed algorithms and the benchmark results reported here, directly work with the scattered profile data without gridding the joint space of computational resource (here, instructions retired, LLC requests, and LLC load misses). For the singlecore case, we illustrate the computational details for an NMPC benchmark. For the multicore case, we provide the results for a benchmark software from real-time systems literature.

The learning-scheduling co-design that builds on the proposed method, will comprise our future work. Another direction of interest is to account for the asynchronous nature of the profiles by formulating the corresponding MSBPs with joint space-time stochasticity, i.e., as a distributional generalization of optimal stopping problems initiated in the recent work [62]. Our multimarginal formulation also opens up the possibility to incorporate additional flow rate constraints as in [63].

With K as in (44) and U as in (35b), consider the Hilbert-Schmidt inner product K, U = ! K 0,s-1 p s (60) where (59) follows by the same argument as for Proposition 4 (see also [28,Proposition 1]). Then, the unimarginal projections (45a) and (45b) follow by applying Lemma 2 to (60).

For the bimarginal projections, let ( j, σ) ∈ [[J]] × [[s]] and σ 1 , σ 2 ∈ [[s]] such that σ 1 < σ 2 . Starting with (59), we have

= p 1 K 0,1

Authorized licensed use limited to: Iowa State University. Downloaded on December 03,2025 at 09:28:11 UTC from IEEE Xplore. Restrictions apply.

•

Then, (46a) follows from (61), and (46b) follows from expanding diag( p σ 1 ) and diag( p σ 2 ) in ( 62), followed by the application of Lemma 2.

With K as in (51) and U as in (35b), we start with the Hilbert-Schmidt inner product K, U = = trace diag(u 1 1 ) X j σ diag(u j σ )Y j σ diag(u 1 s )B j = 1 diag diag(u j σ )Y j σ diag(u 1 s )B j diag(u 1 1 )X j σ = 1 "ƒ‚…

where from projections (53a) and (53b) follow by applying Lemma 2 to (63). Similarly, (53c) follows by applying Lemma 2 to (64).

Next, for the bimarginal projections, we write

where w 1 = w 3 = 1. Now, (54a) follows from (65); (54b) from (66); and (54c) from (67).