UATING DISCRETE LATENT TEMPORAL STRUCTURE
Year: 2019

Abstract:
Clustering methods and latent variable models are often used as tools for pattern mining and discovery of latent structure in time-series data. In this work, we consider the problem of learning procedural abstractions from possibly highdimensional observational sequences, such as video demonstrations. Given a dataset of time-series, the goal is to identify the latent sequence of steps common to them and label each time-series with the temporal extent of these procedural steps. We introduce a hierarchical Bayesian model called PRISM that models the realization of a common procedure across multiple time-series, and can recover procedural abstractions without supervision. We also bring to light two characteristics ignored by traditional evaluation criteria when evaluating latent temporal labelings (temporal clusterings) – segment structure, and repeated structure – and develop new metrics tailored to their evaluation. We demonstrate that our metrics improve interpretability and ease of analysis for evaluation on benchmark time-series datasets. Results on benchmark and video datasets indicate that PRISM outperforms standard sequence models as well as state-of-the-art techniques in identifying procedural abstractions.

heading: 1 INTRODUCTION
text: A fundamental problem in machine learning is the discovery of structure in unsupervised data. In this work we are particularly interested in uncovering the latent structure in procedural data – potentially high-dimensional observational time-series resulting from some latent sequence of events. For example, a video showing how to change a tire might involve jacking up the car, removing one nut from the wheel, removing a second nut, etc. There exists an enormous wealth of videos of procedures available and inferring the latent structure of such videos could be useful in a huge range of applications, from supporting search queries (”find all segments where someone jacks up the car”) to adding to nascent work on learning from observation (Borsa et al., 2017; Stadie et al., 2017; Torabi et al., 2018) in which robots may be taught merely by observing a procedure performed in a video. Procedural learning is related to activity recognition and the broader field of latent temporal dynamical learning, but focuses on the simpler but still ubiquitous setting where the latent activity sequence is a fixed procedure.
As we started to develop methods for performing latent procedure inference from temporal data, we considered how to evaluate our resulting methods. One important evaluation protocol for unsupervised learning methods is external evaluation, where predicted labelings of the discrete variables are compared to ground-truth labels (Rosenberg & Hirschberg, 2007). But the precise criteria for this evaluation is critical. Though there are an enormous number of methods for modeling and inferring latent structure in dynamical systems, including Hidden Markov models, dynamic Bayesian networks, and temporal clustering for time-series data e.g. Fox et al. (2008b; 2014); Krishnan et al. (2015); Linderman et al. (2017); Zhou et al. (2008; 2013); Vidal & Favaro (2014), the external evaluation for such approaches is typically done using the same clustering metrics that are used for non-temporal data: metrics like the normalized mutual information (NMI) (Strehl & Ghosh, 2002) or by computing cluster correspondences using the Munkres algorithm (Munkres, 1957). Unfortunately such metrics disregard temporal information and can therefore make it hard to assess important considerations that arise in latent temporal structure extraction, including in procedural inference.
To address the limitations of prior evaluation criteria for the external evaluation of discrete latent temporal structure, our first contribution is the creation of new criteria. We identify at least two key aspects of temporal latent structure discovery that are not captured well with existing non-temporal metrics: segment structure, the temporal durations between transitions from one discrete latent state to another, and repeated structure, where a particular latent variable may appear in multiple segments. Figure1 shows an example latent sequence associated with a person performing cardiopulmonary resuscitation (cpr) – the colors denote the activity type (such as blowing air into the lungs) and the temporal duration of a fixed activity type is shown as a continuous block of color. Our new metrics can evaluate both segment structure and repeated structure, and we later demonstrate how these metrics may be of broader interest in latent temporal dynamical learning, as they can illustrate and illuminate different strengths and weaknesses of existing work on common benchmark datasets.
Equipped with new measures for evaluating success, we return to our key goal, inferring the latent procedure observed in a set of high dimensional demonstrations. Past expressive models like Hidden Markov models can struggle in this setting, since the resulting procedure may not be Markov: while the setting can be forced to appear Markov by employing the history as part of the state, this can reduce sample efficiency (by artificially increasing the latent space size). Additionally assuming the system is Markov when it is not can mean inferring stochastic latent dynamical structure when it is in fact deterministic. Recent work (Sener & Yao, 2018) tackles a closely related setting, but does not allow repeated primitives, which is a key aspect of many procedures. Instead we introduce a new generative Bayesian model (PRISM) for identifying procedural abstractions without supervision. Results on datasets for procedural tasks – including surgical demonstrations and how-to videos – show that our method can learn procedural abstractions effectively compared with state-of-the-art methods and popular baselines, and outperform prior work on both our new metrics and classic metrics when repeated structure is present.

heading: 2 DEFINITIONS AND SETTING
text: Temporal Clustering. A time-series Xi = (xi1,xi2, . . . ,ximi), i ∈ [n] is a sequence of feature vectors (or items) of length mi with xit ∈ Rd. A temporal clustering C = (c1, . . . , cm) is a sequence of cluster labels with ct ∈ ΓC where ΓC is a set of cluster labels. Temporal clusterings map each item in the feature trajectory X to a cluster label in ΓC . We use Ci as shorthand to denote the temporal clustering generated by some method for Xi. The ground truth temporal clustering for Xi will be denoted by Gi, with corresponding ground truth label set ΓG . For simplicity, we define Ca:b = (ca, ca+1, . . . , cb) where 1 ≤ a ≤ b ≤ m are time indices. For α ∈ ΓC , let C[α] = {t|ct = α, 1 ≤ t ≤ m} be the set of time indices in C whose cluster labels equal α. C[α] is the cluster or partition (in the standard sense) corresponding to the label α. Temporal clusterings are distinct from time-series segmentation (Chung et al., 2004) and changepoint detection (Killick et al., 2012), which consider only identifying boundaries of different dynamic regimes in the time-series.
Segment. A segment is a contiguous sequence of identical cluster labels in C. Each cluster in C may be split across multiple segments. We represent each segment as a pair of time indices (a, b), where a represents the start time-index of the segment while b the end time-index (both included). Formally, we let S(C, α) = {(a, b)|Ca:b = (ca = α, ca+1 = α, . . . , cb = α), 1 ≤ a ≤ b ≤ m} be a function
that maps a temporal clustering C and cluster label α ∈ Γ to the set of segments associated with α in C. We let SC = ⋃ α∈ΓC S(C, α) be the set of all segments in C.
Procedure. Let P be a function that removes running duplicates from a sequence of labels, yielding a sequence of tokens. For instance, P(A,A,A,B,B,B,C,C,C,A,A,B) = (A,B,C,A,B) yielding 5 tokens. For a temporal clustering C, we define its corresponding procedure to be P(C).
The procedure captures the primitives (as tokens) present in the temporal clustering and the sequence in which they occur. We say that C exhibits repeated structure if atleast one primitive is reused in C. This corresponds to having multiple tokens of the same cluster label in P(C). The procedure (A,B,C,A,B) contains repeated structure since both A and B are reused. We also define a weight function W that counts the number of running duplicates in a sequence e.g. W(A,A,A,B,B,B,C,C,C,A,A,B) = (3, 3, 3, 2, 1). Notice that W computes the length of each segment in the sequence. Weights in W(C) thus have a one-to-one correspondence with tokens in P(C),
equaling the length of each token in the procedure. Taken together W(C),P(C) can be used to recover C. Lastly, we letH(P1,P2,W1,W2) be a similarity measure between 2 temporal clusterings represented by P1,W1 and P2,W2 respectively. Fig 2 illustrates these ideas for a temporal clustering defined over the label set Γ = {A,B,C,D}. The procedure exhibits repeated structure since B has (atleast) 2 tokens/segments corresponding to it.

heading: 3 EVALUATING LATENT TEMPORAL STRUCTURE
text: External clustering evaluation is the process of comparing a clustering to a ground-truth labeling in order to evaluate it. Unfortunately, as we discuss below, using traditional clustering evaluation criteria is problematic for evaluating temporal clusterings.
Existing evaluation criteria. The use of standard clustering metrics for evaluation is prevalent in prior work that attempts to evaluate latent temporal structure (Zhou et al., 2008; 2013; Fox et al., 2008b; 2009; 2014; Krishnan et al., 2015; Sener & Yao, 2018; Hoai & De la Torre, 2012). For a temporal clustering C the widely used purity (Rosenberg & Hirschberg, 2007) metric penalizes the presence of items from different ground truth labels in the same cluster. This is a desirable property that a clustering evaluation criterion should have. Purity is defined as purity = 1m ∑ α∈ΓC maxβ∈ΓG |G[β] ∩ C[α]|, A related metric is homogeneity (Rosenberg & Hirschberg, 2007), which captures a similar idea in a different way: homogeneity = 1− H(ΓG |ΓC)H(ΓG) where H(ΓG) is the entropy of the clustering G when partitioned by the labels ΓG , H(ΓG |ΓC) is the conditional entropy of G given C. Intuitively the conditional entropy term looks inside every cluster in C, and checks the entropy of the items inside in terms of the ground truth labels. The conditional entropy will be low (and homogeneity high) if every cluster in C contains items of only a single ground truth label. Another important metric is completeness (Rosenberg & Hirschberg, 2007), which prefers clusterings where all items from a ground-truth label lie in the same cluster i.e. the ground truth label is not split across clusters. Completeness is defined identically to homogeneity, but with the clusterings swapped G ←→ C.
Both completeness and homogeneity are considered important criteria for clustering. Maximizing either at the expense of the other is bad – high homogeneity and low completeness imply a clustering that is too fine-grained (e.g. every item in its own cluster), while low homogeneity and high completeness imply a clustering that is too coarse (e.g. a single cluster containing all items). Therefore, metrics that balance these 2 criterion often perform best (Rosenberg & Hirschberg, 2007) – normalized mutual information (NMI) (Strehl & Ghosh, 2002) and V-measure (Rosenberg & Hirschberg, 2007) are examples of such metrics.
Another popular suite of metrics uses the Munkres algorithm (Zhou et al., 2008) – a one-to-one correspondence between clusters and ground-truth labels is computed so that any classification criterion (commonly accuracy) under this correspondence is maximized.
Characterizing latent temporal structure. Two major characteristics distinguish the evaluation of temporal clusterings from the standard clustering case.
• Repeated structure. The presence of repeated structure in ground-truth – ground-truth labels that are split across several segments in the trajectory. Ideally, we want that this repeated structure is identified correctly e.g. in Fig 3 identifying that the segments A0:2,A4:6 are repetitions of the same primitive.
• Segment structure. The locations and sizes of the segments in ground-truth e.g. in Fig 3 we would like temporal clusterings that switch between primitives at exactly the time points t = 2, 4 with the same 3 segment structure.
Neither repeated or segment structure are seen in standard clusterings, since they treat data as a bag of items.
Limitations of existing criteria. Existing metrics cannot score repeated or segment structure because they convert temporal clusterings to standard clusterings – a contingency matrix (Rosenberg & Hirschberg, 2007) – before evaluation. This erases the temporal information important to score these ideas. In fact, qualitatively distinct temporal clusterings can be mapped to the same contingency matrix, and receive the same score, which is clearly undesirable.
We use Fig 3 as an illustrative example to highlight these shortcomings. In Fig 3, C1 perfectly matches G, while both C2 and C3 split the cluster G[A] – C2 into C2[X] and C2[Z] and C3 into C3[X] and C3[Z]. In addition, C2 also misses the repeated structure in A while C3 does not. C3 captures this repeated structure as XZ, with X0:1Z1:2 corresponding to A0:2 and X4:5Z5:6 corresponding to A4:6. However, it does so at the cost of creating more segments than C2 – ignoring labels, C2 has a segmentation exactly like G. Intuitively, we prefer C1 > C3 > C2 if we care more about repeated structure, and C1 > C2 > C3, if we care more about the segmentation quality. However, all 3 temporal clusterings C1, C2, C3 in Fig 3 are assigned a perfect purity or homogeneity, since as seen in Fig 3b, no cluster in any of the 3 contains items from more than one ground-truth label. Moreover, no traditional clustering criteria can distinguish C2 from C3 since their cluster composition is identical, as seen in Fig 3b. We now discuss new evaluation criteria that can systematically evaluate and score both repeated and segment structure.

heading: 3.1 EVALUATING REPEATED STRUCTURE
text: We describe a new evaluation criteria, called the repeated structure score (RSS) below. Algorithmically, the calculation of the RSS requires the following steps,
1. Pick a ground-truth label β ∈ ΓG and find the set of segments S(G, β) marked with β in G. These segments all exhibit repeated structure in G.
2. Find Ca:b for every (a, b) ∈ S(G, β). This is the predicted temporal clustering restricted to a particular segment (a, b) in G.
3. For every such pair Ca1:b1 , Ca2:b2 , use a scoring functionH to computeH(Ca1:b1 , Ca2:b2). 4. Aggregate these scores across all such pairs for every ground-truth label β ∈ ΓG ,
RSS =
∑ β∈ΓG ∑ (a1,b1),(a2,b2)∈S(G,β)H(Ca1:b1 , Ca2:b2)
2 ∑ β∈ΓG |S(G, β)| ∑ (a,b)∈S(G,β)(b− a+ 1)
where the denominator normalizes RSS to lie in [0, 1].
(a) Example temporal clustering.
(b) Scoring C0:2, C4:7 with P(C0:2) = (X), W(C0:2) = (2), P(C4:7) = (X,Z), W(C4:7) = (1, 2) and running H yields heaviest sub-sequence (X) with score 3. C5:7 (dotted) does not contribute to the score.
Figure 4: Example for repeated structure scoring.
As a concrete example, in Fig 4a, we could pick A in G for step 1 with S(G, β = A) = {(0, 2), (4, 7)}. For step 2, we find C0:2 = (X,X) and C4:7 = (X,Z,Z). For step 3, we must choose an appropriate scoring function H. While many choices for H are possible, we let H(Ca1:b1 , Ca2:b2) be the total weight of the heaviest common sub-sequence or substring (Jacobson & Vo, 1992) in Ca1:b1 and Ca2:b2 . Fig 4b uses this choice ofH to score the pair of segments we picked – we run H ((X,X), (X,Z,Z)) and it returns a score of 3. Intuitively,H tries to match procedures P(Ca1:b1) to P(Ca2:b2) – finding the (heaviest) sequence of tokens common to both P(Ca1:b1) and P(Ca2:b2), while respecting the temporal ordering of both procedures. At one extreme, if H evaluates to 0 then no overlapping tokens exist in P(Ca1:b1) and P(Ca2:b2). This implies that C failed to identify repeated structure in the two segments and a score of 0 is appropriate. At the other extreme, if both segments follow identical procedures P(Ca1:b1) = P(Ca2:b2), they receive the maximum possible score according toH; reflecting that C identified repeated structure perfectly. H can be computed efficiently using a dynamic program in O(|P1||P2|).

heading: 3.2 EVALUATING SEGMENT STRUCTURE
text: We first introduce a new information-theoretic criteria – the label-agnostic segmentation score (LASS). The main purpose of LASS is to score how well C predicts the location of transition points between segments in G. We want to identify and score cases where C either oversegments (introduces extra segments) or undersegments (omits segments) the time-series compared to G. In any temporal clustering, the transition boundaries do not depend on the actual labeling, but only on where the labels switch. Therefore, LASS is label-agnostic and only requires the transition points between segments.
Prior work in changepoint detection and time-series segmentation (Killick et al., 2012) has attempted to quantify this using detection criteria – e.g. if a changepoint close to one in ground-truth is detected, a true positive is recorded. However, these criteria (i) require specifying the tolerance interval for detection, and (ii) don’t degrade gracefully with the segmentation, which limits their applicability to temporal clusterings.
Recall that SG ,SC are the set of all segments in G and C. To compute LASS, the main quantity we rely on is H (SC |SG) – a conditional entropy over segments in C. This can be written as,
H (SC |SG) = − ∑
(a,b)∈SG
(b− a+ 1) n ∑ (ã,b̃)∈SCa:b (b̃− ã+ 1) (b− a+ 1) log (b̃− ã+ 1) (b− a+ 1)
Intuitively, for each ground-truth segment (a, b) ∈ SG , we compute the (weighted) entropy of C restricted to (a, b). This entropy will be > 0 only if C creates new segments in (a, b). On the other hand, to detect under-segmentation, we can simply invert the roles of G and C by calculating H (SG |SC). Finally, noting that H (SC |SG) ≤ H (SC), we define LASS := 1 − H(SC|SG)+H(SG |SC)H(SC)+H(SG) . Oversegmentation (and under-segmentation) can also be measured separately by defining LASS−O := 1− H(SC|SG)H(SC) and LASS−U := 1− H(SG |SC) H(SG)
. The algebraic form chosen for the criteria relates it back to existing information-theoretic criteria such as the NMI (Vinh et al., 2010).
Next, we extend the completeness metric defined earlier. While completeness is desirable for any temporal clustering, using completeness directly has an unintended drawback – it penalizes repeated structure. An example of this can be seen in Fig 3a, where C3 is assigned the same completeness score as C2 despite capturing repeated structure. Since the RSS already provides a systematic way of evaluating presence/absence of repeated structure, we modify the completeness score slightly. We define segmental completeness to be 1− H(ΓC|SG)H(ΓC) , where the conditioning is now on ground-truth
segments rather than on ground-truth labels. Equivalently, we can view this as relabeling each segment in ground-truth to be distinct, and then computing completeness. In this way, segmental completeness serves a similar function to completeness without interfering with repeated structure scoring. Homogeneity can also be extended to define segmental homogeneity similarly.
Lastly, we define an aggregate segment structure score (SSS),
SSS := 1− H(SC |SG) +H(SG |SC) +H(ΓC |SG) +H(ΓG |SC) H(SC) +H(SG) +H(ΓC) +H(ΓG)

heading: 3.3 A COMBINED EVALUATION CRITERIA
text: We have introduced several new metrics: (i) the RSS systematically scores repeated structure; (ii) the LASS scores how well the transition structure is captured; (iii) the SSS provides a unified metric for assessing segment structure. Note that all of these metrics obey the useful property of n-invariance (Rosenberg & Hirschberg, 2007; Meilă, 2007) – scaling the number of items (stretching the temporal clustering) does not change the metric values. Finally, we provide a single evaluation measure for temporal clusterings, the temporal structure score (TSS) that balances the RSS with the SSS.
TSS := (1 + β) · RSS · SSS β · RSS + SSS
While it is preferable to examine the constituent scores individually to determine the strengths and weaknesses of methods, a unified metric allows us to choose methods that balance these criteria. β can be tuned to change the influence of the constituent metrics according to the desired application – we set β = 1.0 in a problem-agnostic fashion in this paper. This choice is in line with previous work that defines compound criteria, such as Rosenberg & Hirschberg (2007).

heading: 4 IDENTIFYING PROCEDURAL ABSTRACTIONS
text: So far, we have argued the importance of evaluating repeated and segment structure correctly in temporal data. Our new metrics enable a systematic evaluation of these ideas, whereas existing metrics did not specifically evaluate or emphasize them. Given this, we now return to our original goal of identifying procedural abstractions in time-series data, which we term procedure identification.
Concretely, we assume a dataset of time-series such that each demonstrates the same task or phenomenon. Crucially, we assume that all time-series share this single, latent procedure – different time-series are distinct realizations of this underlying procedure. The goal is to label the temporal structure present in each time-series while identifying the common, latent procedure. For example, this may involve identifying the common sequence of steps in a set of demonstrations that follow the same cooking recipe as well as their temporal extent in each demonstration.
Procedure identification is related to recent efforts in computer vision that relate to unsupervised activity recognition in videos (Wu et al., 2015; Yang et al., 2013; Krishnan et al., 2015; Sener & Yao, 2018) as well as prior work in latent variable modeling (Yang et al., 2013; Fox et al., 2008b; 2014; Johnson & Willsky, 2013; Linderman et al., 2017). These papers typically assume a general setting where activities or latent tokens may be reused within or across multiple time-series but
there is no onus on methods to discover this structure. Unlike prior work, procedure identification requires that methods also return a common procedure in addition to the temporal clustering of each time-series. This additional requirement places a burden on the method to explicitly, and not implicitly reason about repeated structure.
Formally, finding a common procedure P corresponds to finding temporal clusterings Ci for each time-series, such that P = P(C1) = P(C2) = · · · = P(Cn). Since a single procedure consistent with all temporal clusterings is desired, shared primitives must be identified across all videos. This requirement also makes evaluation in this setting particularly suited for our metrics. Methods that
cannot relate primitives across demonstrations (identify repeated structure) or reason about their temporal extents (identify segment structure) will be unable to recover the underlying procedure and should therefore be penalized. There may also be additional repeated structure within a time-series, for example if someone mixes cake batter, adds ingredients, and then mixes the cake batter again.
A closely related setting was tackled by the method outlined in Sener & Yao (2018), who assume that each time-series is a permutation of a shared bag of activities. However, as they acknowledge, a major limitation of their method is that they allow each activity to occur atmost once within a time-series, which limits repeated structure. Our problem setting allows for this possibility, which is often seen in real-world datasets e.g. the JIGSAWS surgical dataset (Gao et al., 2014) contains videos of surgeons doing multiple repetitions of primitives such as needle-passing or suturing. However, we do not allow the procedure to be permuted in our problem setting, which is a useful extension that we hope to explore in future work.

heading: 4.1 PRISM: PROCEDURE IDENTIFICATION WITH A SEGMENTAL MIXTURE MODEL
text: We now describe a statistical model called PRISM for the procedure identification problem. PRISM is a generative hierarchical Bayesian model that explicitly separates reasoning about the procedure from its realization in each time-series in the dataset. It utilizes this to create a mixture distribution over primitives for every observation in the dataset. Fig 5 illustrates the model.
Unlike popular sequence models based on the Hidden Markov model, PRISM does not make a restrictive Markov assumption and can handle any procedure. And in contrast to a Gaussian mixture model, PRISM shares statistical strength across time-series by using a structured prior over the (discrete) local variables that generate each observation – amortizing the procedure across many timeseries. As we show later in experiments, this structure can help PRISM outperform these methods, especially when the observations are noisy.
Modeling the procedure. PRISM assumes that each time-series Xi follows the same underlying procedure p = (p1, . . . , ps), where s is the number of steps in the procedure. p is a vector of discrete, global latent variables in the model. This implies that each time-series will split into s segments, each segment corresponding to a single step in the procedure.
Next, let |ΓC | = K be the number of cluster labels (primitives). Each step in the procedure will correspond to one of these K labels i.e. pr ∈ [K], 1 ≤ r ≤ s. Thus, p is an ordered sequence of high-level primitives e.g. to make tea, we would boil water → get cup→ pour water and so on. Formally, we assume a generation process,
λr ∼ Dirichlet(α1, . . . , αK), 1 ≤ r ≤ s pr ∼ Categorical(λr), 1 ≤ r ≤ s
Each token in the procedure is independently generated from a separate Dirichlet prior λr, which allows fine-grained control over the distribution of tokens at each step. Note that we do not disallow self-transitions, which can be utilized by the model to coalesce adjacent segments.
Modeling realizations in each time-series. Different time-series can be vastly different in how they implement the procedure e.g. they can be videos of different individuals making tea. These individuals, while following the same procedure, may spend varying amounts of time on each step (or even skip steps).
In PRISM, the realization of the procedure in the ith time-series (having mi time-steps) is given by wi = (wi,1, . . . , wi,mi), a sorted sequence of discrete random variables. Each wi,j ∈ [s] is an index into the procedure p. Its value dictates the step being followed in the procedure at the jth time-step in the ith time-series. We assume wi is generated as follows,
γ ∼ Dirichlet(β1, . . . , βs) ui,j ∼ Categorical(γ), 1 ≤ i ≤ n, 1 ≤ j ≤ mi wi = Sort(ui,1, . . . , ui,mi), 1 ≤ i ≤ n
Mathematically, to generate wi, we first generate an unordered bag of mi independent categorical variables ui,j and then sort them (in ascending order). Each ui,j is generated from the same Dirichlet prior γ. Intuitively, if we use the stick-breaking metaphor for the Dirichlet distribution, γ is a prior over how long each token in the procedure is – γr is the expected relative length of pr. Thus, γ allows us to impose a prior on the relative length of each step of the procedure, in every time-series. An equivalent view into our model is that each wi is directly generated by a Multinomial distribution
with mi draws on probabilities dictated by γ. However, we find that the view that we adopt yields a far more efficient inference scheme (see Appendix).
Generating observations. For every time-series, we can combine the common procedure p and the realization of the procedure wi to recover a temporal clustering zi for that time-series. We write this temporal clustering as zi = (zi,1, . . . , zi,mi) where zi,j ∈ [K] is a discrete primitive. For each primitive k ∈ [K], we assume a Gaussian observation model (µk,Σk) (with a Normal-Inverse Wishart prior) shared by all time-series. PRISM then generates each xij independently from its local assignment zi,j and the appropriate observation model. We use Gibbs sampling (see Appendix) to perform posterior inference of the latent variables ({pr}r∈[s], {ui,j}i∈[n],j∈[mi], {µk}k∈[K], {Σk}k∈[K]) in the model.

heading: 5 EXPERIMENTS WITH EVALUATION CRITERIA
text: We first compare our metrics to existing criteria by evaluating several competing methods on real datasets. We note that the methodology for validating our evaluation criteria is similar to that followed by Rosenberg & Hirschberg (2007), with the exception that we compare criteria on real-world data with a large number of methods. We compare to homogeneity (Hom), completeness (Com), NMI since they give the best combination of performance for clustering evaluation (Rosenberg & Hirschberg, 2007). We also compare to accuracy computed using the Munkres method since it is widely used e.g. (Sener & Yao, 2018; Zhou et al., 2008) and the adjusted Rand index (ARI) (Vinh et al., 2010).
Methods. We use several methods to generate temporal clusterings: Hierarchical Clustering (AGG); Hidden Markov Models (HMM); Hierarchical Dirichlet Process HMMs (HDP-HMM) (Teh et al., 2005); Sticky HDP-HMM (SHDP-HMM) (Fox et al., 2008a); HDP Hidden Semi-Markov Models (HDP-HSMM) (Johnson & Willsky, 2013); Switching Linear Dynamical Systems (SLDS) (Fox et al., 2008b); HDP-SLDS (Fox et al., 2008b); Sticky HDP-SLDS (SHDP-SLDS) (Fox et al., 2008b); Ordered Subspace Clustering (OSC) (Tierney et al., 2014); Temporal SC (TSC) (Li et al., 2015); Low-Rank SC (LRSC) (Vidal & Favaro, 2014); Sparse SC (SSC) (Elhamifar & Vidal, 2013) and Spatial SC (SPATSC) (Guo et al., 2013).
Datasets. We use 2 common benchmark datasets (Fox et al., 2008b; 2009; 2014; Zhou et al., 2008; 2013). BEES consists of 6 time-series, each recording the dancing behavior of a bee. Each time-series is built up of 3 primitive movements – waggle, left-turn and right-turn – repeated several times. For BEES, we treat each sequence as a separate dataset (BEES1,. . . ,BEES6) and run each method on them individually, as done in Fox et al. (2008b). Here, we focus on results on BEES1 but results on MOCAP6, a motion-capture dataset and all BEES datasets can be found in the Appendix.
Results. For each of our metrics, we first illustrate that they qualitatively capture the characteristics that they were designed for. Fig 6 visualizes the best and worst performing methods on BEES1 for each of our metrics. In Fig 6a, SSC is worst on LASS-O since it heavily over-segments the time-series by introducing many transitions, while TSC introduces far fewer transitions. In Fig 6b, under-segmentation is penalized and we can see that AGG should be worse since it contains longer segments and far fewer transitions than SSC. Combining these ideas in Fig 6c yields that HDPHSMM has the best LASS – observe that it contains fewer closely spaced transition points and creates a more even segmentation at the right level of granularity, explaining its score.
Fig 6e shows the result of evaluating with just the RSS – we can clearly see that SpatSC is able to capture far more repeated structure than SHDP-HMM, which is unable to label segments from the same ground-truth label consistently. Lastly, we see that SpatSC is judged to be the best method on BEES1 by TSS in Fig 6f since it has the highest RSS and a SSS of 0.81 which is close to the highest.
Improving scores on evaluation metrics is common practice to demonstrate performance gains due to a method. We now highlight cases where TSS disagrees with NMI, ARI and/or Munkres in the assessment of competing methods, since this discrepancy has significant impact on which method is picked. We find that such disagreements happen very often across all datasets. An illustrative example for BEES1 is shown in Fig 6h: TSS and Munkres agree that SSC is better than HDP-HSMM, but NMI and ARI strongly disagree. Before diving into why this may be happening, we note that examining the RSS and SSS allows us to understand why TSS considers SSC to be better than HDP-HSMM – it captures more repeated structure. On the other hand Munkres, while it agrees with
TSS, is unable to give us this insight. We find that NMI is highly sensitive to the number of clusters – SSC only uses 2 clusters, which means it is heavily penalized by NMI, in spite of capturing some repeated structure. ARI is less sensitive than NMI, but is far less interpretable.
Another example is shown in Fig 6g, where the TSS for both SLDS and HDP-HMM is comparable. Interestingly, NMI/ARI and Munkres are all highly skewed and disagree strongly on the better method. Since HDP-HMM has more clusters, the clusters are more homogenous and Hom is much higher, pushing up NMI. On the other hand, since Munkres finds a one-to-one correspondence between clusters and ground-truth labels, having more clusters means some cannot be matched. This causes the score to strongly downweigh suchh clusterings. Overall, our metrics provide both interpretability and the ability to analyze where the method is falling short.
Lastly, we conduct a sensitivity analysis for β as follows: assume that the most appropriate value for evaluation is β′, but we instead use the default value of β = 1.0. Our main idea is to compute how well the ranking of methods under β = 1.0 captures that under β′. To do this, we use the Normalized Discounted Cumulative Gain (NDCG), a common criteria from information retrieval. Varying β′, we see in Fig 6i that β = 1.0 outputs similar rankings to a relatively wide range of other β values, other than at the extremes. This confirms that most settings of β will be captured well by our choice of β, unless there is a strong reason to prefer either the RSS or the SSS only.

heading: 6 EXPERIMENTS WITH PRISM
text: We now evaluate PRISM on complex datasets with noisy features, including video data. For PRISM, we set α = 1.0,β = 0.1 (except BEES) and we set K (number of mixture components) to equal the number of ground-truth labels for all methods. Results are averaged over 10 random seeds. Experimental details including hyperparameters can be found in the Appendix.
Before our main results, we briefly conduct a simple check on performance for the BEES dataset introduced earlier. Though this only has a single example per procedure, we find that our approach still does well (Table 1) compared to two simple baselines – all approaches use the same observation model with Gaussian emissions.
Leveraging a fixed procedure. We now evaluate the benefit of our approach on the JIGSAWS dataset (Gao et al., 2014), which consists of surgeons with varying skill demonstrating common procedures, such as suturing, needle-passing and knot-tying. We use kinematic features that correspond to the grippers controlled by the surgeon for all methods. For each surgical task, we select the 10 demonstrations that correspond to the 2 expert surgeons since they adhere most closely to the prescribed procedures. Even for a single expert, there is considerable variability across the 5 demonstrations they generate, so this is far from an idealized setting. Once again, all methods leverage the same observational model with the same features. Our baselines are two standard approaches – a Gaussian mixture model (GMM), and variants of a Bayesian Hidden Markov model.
As expected, PRISM outperforms these baselines by a considerable margin. On the suturing task, performance is improve by almost double-digits on both metrics, with performance gains also visible on the other tasks. Leveraging the common procedure allows PRISM to share statistical strength
across demonstrations and make consistent predictions where possible, whereas the other methods must rely more strongly on the observation model to infer the latent sequence. Without very cleanly separated features for different primitives, they can often be led astray.
We also experiment with a smaller subset of 5 demonstrations provided by a single expert on the suturing task. This subset of data has the least amount of variability in terms of the expert’s procedure. Because of this, we see that PRISM gives an improvement of ∼ 18 points on NMI and ∼ 10 points on the TSS. To illustrate PRISM’s sensitivity to hyperparameters, we vary s and K in Fig 7 on this single-expert suturing task. We find that PRISM’s performance is quite stable with respect to both – this holds as long as s is larger than the number of segments in ground-truth, and K is close to the true number of ground-truth clusters.
Noisy procedures with video data. Next, we consider labeling demonstrations that are very noisy in their adherence to the procedure, with a great deal of variability in both the steps being followed, their order, as well as their relative lengths. Recent work by Sener & Yao (2018) achieved state-of-the-art performance without any form of supervision on two large video datasets – Breakfast actions (Kuehne et al., 2014) and INRIA instructional videos (Alayrac et al., 2016). We expect our method to perform worse on the Breakfast actions dataset, as it has no repeated structure, in which case our model has no advantages over this prior model, and indeed uses a simpler observation model. Perhaps encouragingly, in the Breakfast actions dataset, we find that our method only performs slightly worse than but close to Sener & Yao (2018), achieving a Munkres score of 33.5 compared to their score of 34.6. It is likely that this can be improved by tuning the number of segments for each procedure (we set s = 20 for all procedures).
However a key limitation of their work is that they do not model repeated structure. Therefore in domains with repeated structure we expect our approach to do well. Such structure is present in some of the activities in the INRIA dataset. We compare to Sener & Yao (2018)’s results on the INRIA dataset in Table 2. PRISM outperforms their best method on 3 out of the 5 activities despite not explicitly modeling background information, with a significant difference on cpr and jump car. On their comparable method which does not model background, PRISM is significantly better across all procedures. This is despite PRISM using a much weaker observation model (they simultaneously learn a feature representation for the videos) and no permutation model. An example for the cpr procedure is shown in Fig 1 and illustrates that unlike Sener & Yao (2018)’s method, PRISM can successfully recover repeated structure both within and across videos. We expect that incorporating a background model in PRISM could enable our approach to yield even stronger performance.

heading: 7 CONCLUSION
text: We developed PRISM, a hierarchical Bayesian model to identify latent procedures in time-series data. Results on several datasets show that PRISM can be used to learn procedural abstractions more effectively than competing methods. We also introduced new evaluation criteria for the problem of externally evaluating temporal clusterings. Our metrics address gaps in temporal clustering evaluation, and provide both interpretability and the potential for easier analysis. In future work, we hope to extend PRISM by incorporating richer observation models, and relaxing the assumption that a single strictly-ordered procedure is followed by all time-series. Links to code can be found in the Appendix.

heading: 8 ACKNOWLEDGMENTS
text: The authors would like to thank the Schmidt Foundation and a NSF Career Grant for their support.

heading: A ALGORITHM FOR REPEATED STRUCTURE SCORE
text: Detailed pseudocode for reproducing the computation of the RSS is given in Algorithm 1.
Algorithm 1 Calculation of RSS Require: Ground truth G, temporal clustering C, similarity functionH
1: \\ Optional Pruning Step 2: for all α ∈ ΓC do β∗α ← arg maxβ∈ΓG |G[β] ∩ C[α]| . find β∗α which most overlaps α 3: for all (a, b) ∈ S(G, β), β ∈ ΓG do . for every segment of every β in G 4: P(Ca:b),W(Ca:b) := (p1, p2, . . . , pr), (w1, w2, . . . , wr) . locate segment in C 5: w̃i ← wi × I[β = β∗pi ] . trim weights based on purity 6: W̃(Ca:b)← (w̃1, w̃2, . . . , w̃r) . construct new weights 7: \\ Repeated Structure Scoring 8: score← ∑ β∈ΓG ∑ (a1,b1),(a2,b2)∈S(G,β)H(P(Ca1:b1),P(Ca2:b2), W̃(Ca1:b1), W̃(Ca2:b2))
9: . score every repeated structure pair for every ground-truth label β 10: max score← 2 ∑ β∈ΓG |S(G, β)| ∑ (a,b)∈S(G,β)(b− a+ 1) 11: return RSS← scoremax score

heading: B EXPERIMENTS WITH EVALUATION CRITERIA
text: The results presented for the evaluation criteria developed in this paper can be reproduced using code available at https://github.com/StanfordAI4HI/ICLR2019_evaluating_ discrete_temporal_structure. The temporal clusterings predicted by each method, as well as ground-truth data is available at the link.
In addition, an implementation of only the newly proposed temporal clustering evaluation criteria can be found at https://github.com/StanfordAI4HI/tclust-eval.
For completeness, we include the result of evaluation across all methods, datasets and criteria.

heading: C DETAILS FOR PRISM
text: This section contains details related to both the modeling and inference performed by PRISM.
C.1 GENERATION IN PRISM
The generative process used in the hierarchical Bayesian model is given below,
λr ∼ Dir(α1, . . . , αk), 1 ≤ r ≤ s γ ∼ Dir(β1, . . . , βs) pr ∼ Cat(λr), 1 ≤ r ≤ s ui,j ∼ Cat(γ), 1 ≤ i ≤ n, 1 ≤ j ≤ mi wi = Sort(ui,1, . . . , ui,mi), 1 ≤ i ≤ n zi = TemporalClustering(p,wi), 1 ≤ i ≤ n (µk,Σk) ∼ NIW(µ0,Σ0, κ0, ν0), 1 ≤ k ≤ K xi,j ∼ Normal(µzi,j ,Σzi,j ), 1 ≤ i ≤ n, 1 ≤ j ≤ mi
C.2 INFERENCE IN PRISM
We use Gibbs sampling to perform posterior inference of the latent variables ({pr}r∈[s], {ui,j}i∈[n],j∈[mi], {µk}k∈[K], {Σk}k∈[K]) in the model.
Sampling ui,j . Let Φ = ( {µk}k∈[K], {Σk}k∈[K] ) , and write the complete conditional for ui,j ,
Pr(ui,j = r|u−i,u−ji , . . . ) ∝ Pr(ui,j = r|β,u−i,u −j i ) mi∏ j=1 Pr(xi,j |ui,j = r,u−ji ,p,Φ)
where in the first term we have collapsed the Dirichlet prior γ. The main difficulty is in the datalikelihood terms Pr(xi,j |ui,j = r, . . . ), where we must reason about the effect of changing ui,j . To reason about this effect, we must look at w(r)i = Sort(ui,j = r,u −j i ), combine it with p to yield z (r) i and then compute the observation-likelihood terms. Due to the sorting step that yields w(r)i , changing a single ui,j affects the entire sequence z (r) i . Naively, we could compute Pr(xi,j |ui,j = r, . . . ) for every j ∈ [mi], r ∈ [s] but this would be extremely-inefficient (compute s ·mi terms).
Instead, we note that the assignments z(r)i will be almost identical for every r ∈ [s] since u −j i is held fixed and only ui,j is varied. The only points at which the z (r) i are not identical are those right before and right after segment boundaries – O(s) such locations. We only need to compute the data-likelihood at these locations since no other locations are changed by sampling ui,j . Using this insight, we only need to compute s2 terms, a big reduction since we expect s << mi. We sample all the ui,j’s sequentially using this update.
Sampling pr. To sample pr, we write the complete conditional,
Pr(pr = k|p−r,α, . . . ) ∝ Pr(pr = k|α) n∏ i=1 mi∏ j=1 Pr(xi,j |pr = k,p−r,wi,Φ)
Let z(k)i be the temporal clustering generated by combining (pr = k,p−r) with wi. Note that in each time-series, setting pr = k affects only those observations that are assigned to pr by wi. We
thus need to compute the observation-likelihood of only these terms,
Pr(pr = k|p−r,α, . . . ) ∝ Pr(pr = k|α) n∏ i=1 mi∏ j=1 zi,j=k Pr(xi,j |µk,Σk)
We sample the pr’s sequentially using this update.
Sampling µk,Σk. Sampling for the observation model can be done directly from the posterior of the Normal-Wishart Inverse distribution. We set uninformed priors for the observation model, µ0 = 0d,Σ0 = Id×d, κ0 = 1, ν0 = d+ 2.

heading: D EXPERIMENTS WITH PRISM
text: The results presented for PRISM can be reproduced using code available at https://github. com/StanfordAI4HI/ICLR2019_prism, which contains a Python implementation of PRISM.
D.1 DATASETS
We now provide details for preprocessing each dataset, with data loaders available at the linked code.
JIGSAWS. For the JIGSAWS dataset, we use the first 38 kinematic features. We subsample every 3rd timestep. To preprocess the data, we perform Principal Components Analysis (PCA) for dimensionality reduction to 12 features. We then stack features from 10 time-steps and perform PCA again to reduce dimensionality to 10 features. Lastly, we standardize features to have zero mean and unit variance.
Breakfast actions. We use the Fisher vector features taken from the ‘stereo01’ camera without any other preprocessing. We use the coarse segmentation for evaluation.
INRIA instructional videos. We use PCA to reduce the dimensionality of the features to 64.
Bees. We use 4 features – x, y, cos θ and sin θ.
D.2 EXPERIMENT DETAILS
We use the hyperparameter settings given in Table 10 for PRISM. The Bayesian HMM has a single hyperparameter α which represents the hyperparameter for the Dirichlet prior over the transition matrix.
D.3 RESULTS ON BREAKFAST ACTIONS
Detailed results for the Breakfast actions dataset are summarized in Table 11.
D.4 NON-MARKOV SIMULATION STUDY
We generate 10 time-series using the following simple non-Markov sequence,
AABBAACCAADDAAEEAAFFAAGGAAHHAAAAAAAA
There are 8 distinct latent tokens in the sequence, and we generate noisy Gaussian 2-dimensional observations for each such token, using a mean of (i, i) for the ith token, and a standard-deviation of 0.35 in either dimension.
We compare the Hidden Markov model variants to PRISM and the Gaussian mixture model. All methods use the same observational model with Gaussian emissions.
We expect that the main advantage of using an HMM vis-a-vis a GMM for recovering the latent temporal clustering, is in cases where the observations are noisy. In such cases, the learned transition dynamics can override the (noisy) observational likelihood to yield the correct sequence of tokens, as long as the procedure is Markov. However, this benefit should disappear in cases where the procedure is not Markov, like the sequence constructed by us above.
In contrast to these methods, PRISM makes no Markov assumption about the underlying procedure. This allows it to learn any procedure, Markov or otherwise. In addition, by sharing this procedure across multiple time-series, PRISM can avoid being led astray by noisy observations (unlike the GMM, which relies purely on the observation to infer the local latent token).
This argument is confirmed by the results in Table 12. PRISM outperforms all methods on both evaluation metrics, while the HMM variants suffer due to their inability to abandon the Markov assumption that is baked into the model.

Review summaries:
In "Learning procedural abstractions and evaluating discrete latent temporal structure" the authors develop a hierarchical Bayesian model for patterns across time in video data. They also introduce new metrics for understanding structure in time series (completeness and homogeneity). This work is appropriate for ICLR. They provide some applications to robotics, suggesting that this could be used to teach robots to act in environments by learning from videos.