extfieldsofstudy:
Computer Science
Engineering
Mathematics

text:
End-to-End Learning of Semantic Grasping

We consider the task of semantic robotic grasping, in which a robot picks up an object of a user-specified class using only monocular images. Inspired by the two-stream hypothesis of visual reasoning, we present a semantic grasping framework that learns object detection, classification, and grasp planning in an end-to-end fashion. A"ventral stream"recognizes object class while a"dorsal stream"simultaneously interprets the geometric relationships necessary to execute successful grasps. We leverage the autonomous data collection capabilities of robots to obtain a large self-supervised dataset for training the dorsal stream, and use semi-supervised label propagation to train the ventral stream with only a modest amount of human supervision. We experimentally show that our approach improves upon grasping systems whose components are not learned end-to-end, including a baseline method that uses bounding box detection. Furthermore, we show that jointly training our model with auxiliary data consisting of non-semantic grasping data, as well as semantically labeled images without grasp actions, has the potential to substantially improve semantic grasping performance.

Introduction
Robotic manipulation presents a number of challenges as well as exciting opportunities for computer vision. Robots that interact with the world must recognize and localize objects with a very high degree of precision, which places a high demand on the performance of computer vision systems. The substantial domain shift between objects and scenes seen by a robot and those in labeled object recognition datasets makes it nontrivial to generalize conventional computer vision systems to robotic tasks. On the other hand, robots have the ability to autonomously collect their own data as they perform manipulation tasks in the real world. By developing methods for self-supervision, we can utilize robots to autonomously collect large self-labeled datasets that can help us exceed the performance of standard vision systems that are trained on fixed-size, human-labeled datasets.
In this paper, we explore this in the context of robotic grasping with semantic labels. The purpose of our method is to allow a robot to successfully grasp objects in its environment based on highlevel user commands, which specify the type of object to pick up. This task aims to combine both spatial and semantic reasoning: the robot must determine which of the objects in the scene belong to the requested class, and how those objects might be grasped successfully. We aim to answer two questions: first, how can we design a model that effectively combines both spatial and semantic reasoning, and second, how can we exploit the capacity of the robots to collect self-supervised data while also making use of modest numbers of human-provided labels to learn about semantics. Answering these questions successfully could lead not only to more effective and useful robot systems, but to advances in spatial and semantic visual reasoning and automated data gathering.
Previous work [17,21] has demonstrated that a non-semantic grasp outcome predictor can be trained in an entirely self-supervised manner via autonomous data collection. However, those works did not Figure 1: (a) We consider the task of learning to pick up objects from 16 object classes, conditioned on a user-specified class. (b) Our robotic setup is an arm with a 2-finger gripper that servos using monocular camera images. (c) Inspired by the two-stream hypothesis of vision [26], we propose a model that is divided into a dorsal and ventral stream: the dorsal stream selects which actions result in good grasps, while the ventral stream selects the class of the object being grasped. Best viewed in color.
consider object class in grasp planning. In order to solve the semantic grasping task, we require semantically labeled data and a model that is capable of reasoning about object identity.
In designing our method, we take inspiration from the "two-stream hypothesis" of human vision, where visual reasoning that is both spatial and semantic can be divided into two streams: a ventral stream that reasons about object identity, and a dorsal stream that reasons about spatial relationships without regard for semantics [26]. We train a deep neural network for performing the semantic grasping task that is decomposed into these two streams, as illustrated in Figure 1. The dorsal stream predicts whether a particular grasping action will result in a successful grasp. The ventral stream is trained to predict what type of object will be picked up. We tackle the problem of assigning semantic labels to large amounts of robot-collected data via a label propagation scheme: if the robots successfully grasp an object, they can then bring it up to their cameras for closer inspection. A modest number of such human-labeled images can be used to train a classifier for these images, and then propagate labels to the images used for learning grasping in a semi-supervised fashion.
The main contribution of our work is a two-stream neural network model for learning semantic grasping from off-policy data. We use a combination of self-supervised robotic data collection and label propagation to construct large object classification dataset for training our models. Finally, we show that both dorsal and ventral streams can benefit independently from domain transfer, by integrating auxiliary spatial information from non-semantic grasping datasets and auxiliary semantic information from image classification datasets.
Our experimental results demonstrate that: (1) Our proposed two-stream model substantially outperforms a more standard detection pipeline for grasping, as well as single-stream approaches. (2) Integrating auxiliary non-semantic data substantially improves grasping performance. (3) Auxiliary semantically labeled data without grasping actions has the potential to improve generalization, but this improvement is highly dependent on the degree of domain shift between the auxiliary data and the grasping images.
While most grasping methods focus on the spatial aspects of the problem, a number of works have sought to incorporate semantic information into grasping, including affordances and task-specific knowledge [2,14,9], instance-based object models [19], known shape information about objects [6], and prior knowledge that objects of the same type should be grasped in similar ways [18]. However, while these prior methods incorporate semantic knowledge to aid spatial reasoning [15], our approach instead aims to combine both spatial and semantic reasoning to pick up objects of a particular category in cluttered scenes. Incorporating additional feedback from semantic knowledge into the spatial problem of grasping, as in prior work, is also a very valuable research direction, but largely orthogonal to the main focus of this work. In this sense, our method is complementary.
Self-supervision and large-scale data collection have recently been proposed for image-based robotic grasping [21,17]. Our work also follows this paradigm for training the dorsal stream of our semantic grasping network, with large numbers of grasp attempts collected with real robotic manipulators and labeled as successful or failed grasps automatically. However, we also incorporate a ventral stream into our system that is trained on semantic categories, with knowledge transfer from object classification datasets and automatic label propagation from a small number of human-provided labels. In this sense, our aims are orthogonal to this prior work; our goal is not to simply build a better grasping system, but to explore how semantics and label propagation interact with real-world robotic manipulation tasks.

Deep Semantic Grasping
In our setup, a robotic arm with a two-finger gripper learns to grasp particular objects from a bin based on a user-provided object class. We train our grasping system end-to-end, using off-policy data collected by a non-semantic grasping policy and a modest number of semantic labels provided by humans. We propose several deep neural net architectures for learning semantic grasping, discuss our methodology for collecting large datasets for training our models, and describe how information from external datasets can be leveraged to improve grasping performance.

Two-Stream Value Function
Our policy critic, visualized in Figure 1c, takes as input images I t , I 0 , which correspond to the current image seen by the robot's camera, and the initial image seen during the episode. The initial image is recorded before the grasping episode begins and is included to allow the system to handle self-occlusions caused by the robot arm. The network also receives as input a candidate task-space motion command a t , which is parameterized by a Cartesian displacement of the gripper along with a vertical rotation of the gripper. The network outputs p(g, c|I t , a t ), the joint probability that executing the command a t will result in g, the event that an object is grasped and c, the event that the grasped object is the correct semantic label. Our servoing function samples grasps randomly and uses the neural critic to choose the grasp with the highest probability of picking up the desired object. See Appendix C for additional details on the model architecture and the servoing policy.
A straightforward way to implement a semantic neural critic is to extend the final layer of the model proposed in prior work on self-supervised grasping [17] to output a softmax distribution over object classes, plus an additional "failure" class for when no objects are grasped (Figure 2b). This model, which we call "Single Stream" in our experiments, lacks a dorsal stream that explicitly learns what it means to "pick up any object". Instead, it merges spatial and semantic inference into a single output, predicting grasp success for each object class independently.
Additional supervision of spatial reasoning capabilities can be obtained by distinguishing between failures of classification and failures of detection. We decompose the prediction problem into a dorsal stream p(g|I t , a t ) that is trained on non-semantic grasping data (whose labels are success or failure), while the ventral stream p(c|I t , a t , g) is trained to classify successful grasps whose labels are one of 16 object classes. The model learns to infer where actions will send the gripper, and implicitly attends to that location in order to extract its spatial relevance to grasping (Figure 2a    Hand-eye coordination binds spatial awareness with visual semantics, enabling the policy to learn object classification and detection end-to-end. (b) Single Stream model that predicts grasp success for each class independently, without explicitly modeling object semantics. (c) Two-Stream model that shares model parameters between a grasp branch and class branch, which comprise the dorsal (cyan box) and ventral streams (magenta box). (d) Our best two-stream model uses the class branch from (c) for the ventral stream, but uses a separate dorsal stream trained only on grasp prediction. Best viewed in color.
We consider Two-Stream models as not a single neural net architecture, but rather a family of architectures that functionally decouple the prediction of grasp success (dorsal stream) and object class (ventral stream). An example of a Two-Stream model is a "Two-Branch" (2B) architecture that shares inputs and action processing before branching out to grasp-specific and class-specific parameters ( Figure 2c). Another example is a "Separated" architecture which takes advantage of the fact that Two-Stream models can mix-and-match dorsal and ventral streams at inference time. We train both branches of the 2B architecture as before, but only use the class branch as the ventral stream. The dorsal stream of the Separated architecture is an entirely separate convolutional network trained only to predict grasp success ( Figure 8b). Because this separate dorsal stream outperforms nonsemantic grasping capabilities of our 2B models, we use this separated architecture for most of our Two-Stream models. 2 . Sharing the dorsal stream between Two-Stream architectures also allows for a controlled comparison between different ventral stream architectures.
Instead of using a deep CNN that implicitly attends to an object using the action vector, we wondered whether an explicit attention mechanism would improve the model's ability to extract a representation of visual features that was invariant to the object's location in the bin. In addition to our 16-layer CNN, we also try a shallower 9-layer network that bottlenecks predictions through soft keypoint locations [4]. Both architectures perform similarly, though the soft keypoint architecture has substantially fewer parameters and therefore trains more efficiently. Further details on both architectures are provided in Appendix B.1.

Parallel Data Collection & Model Training
We chose 16 object categories found in typical kitchen and office settings (Table 1). We obtained 30-40 object instances per category and divided these 500 objects into a set of training and test objects   (which are used to evaluate generalization to unseen class instances). To collect the data used to train our semantic grasping models, we first learn a non-semantic grasping policy that is capable of picking up objects in a self-supervised manner, similar to the approach from [17]. We utilized 8 KUKA iiwa arms to collect data in parallel and introduce diverse lighting and variations in camera calibration into the dataset. Labels for grasp success/failure can be generated autonomously via a simple threshold on an image delta, as described in [17].
Next, we use two approaches to annotate the collected dataset with object labels, in order to supervise learning of the ventral branch. In the first approach, each robot's bin is filled with objects sharing the same semantic label. This makes it trivial to assign labels to large amounts of data, for which we gather 10,000 such examples. Object collections are frequently shuffled between the different robots so that our models do not overfit to a specific robot's idiosyncratic visual cues. In our second approach, we program the robot to lift grasped objects close to the camera at the end of a grasping sequence for futher inspection, at which point a close-up image of the object being held in the gripper is recorded (e.g. Figure 3). Presenting objects in this manner isolates the grasped objects for the purpose of easier classification, somewhat analogously to how children pick up and examine objects in order to isolate it from surrounding visual context. We label 20,000 present images using Amazon Mechanical Turk and use this information to train a CNN classifier. By classifying the remaining 300,000 present images and propagating those labels to their corresponding grasping images, we can generate a large autonomously collected dataset for semantic grasping.

Transfer from External Datasets
Although the robots can gather large amounts of data autonomously and then label them via automatic label propagation, the diversity of object instances in such datasets is limited. Fortunately, our two-stream decomposition of spatial and semantic perception affords us a straightforward mechanism for incorporating images from either (1) auxiliary grasping data with no semantic labels, including objects that are difficult to classify or do not belong to the semantic classes, or (2) auxiliary semantic data of different object instances without grasping actions.
For auxiliary grasping data, we jointly train our two-stream models on all the grasping data we collected, including objects that do not belong to our semantic categories, supervising only the grasp branch of the model. In all, we had 11 million grasping images that we could use to train the dorsal stream (corresponding to about 2 million grasps), of which 330,000 had semantic labels and could be included for training the ventral stream. This auxiliary data serves to not only improve our model's ability to generalize to the test objects, but also benefits classification by improving the early visual representations learned by the shared CNN. To incorporate auxiliary semantic data, we include images from object classification datasets directly into the grasping model with a random action.
In our experiments, we present results with two auxiliary semantic datasets, which we call "S1" and "S2." S1 contains images from ImageNet and JFT [3,11] and with present images of training objects. S2 contains present images of objects in the test set (but are not paired with grasping actions). The rationale behind including S2 in our evaluations is to determine how much improvement can be had from a classification dataset that perfectly represents the distribution of objects seen at test time, from a different viewpoint and with a different camera. This represents a "best case" scenario for Figure 4: In our evaluation, the robots are tasked to pick up specific categories from bins containing 5 randomly selected objects from different classes. Top row: Training objects. Bottom row: Test objects. Evaluation objects were selected at random to cover all object categories uniformly.

Experimental Results
We evaluate semantic grasping performance on real-world robots by assigning 5 randomly selected objects from different categories to each robot ( Figure 4) and commanding the robot to pick up the objects in a random order. For a sequence of 5 grasp attempts, the robot removes grasped objects one at a time from the bin and drops them outside the workspace. This is repeated 10 times per robot with 6 robots in total, for a set of 30 training objects and 30 unseen testing objects. As shown in Appendix D, we evaluate policies based on the fraction of grasp attempts that were successful at picking up some object (Objects/Attempts) and the fraction of all grasp attempts that picked up the correct object (Class/Attempts). We demonstrate through baseline comparisons and ablation studies the effect of various architectural decisions in our semantic grasping models.
Summarily: (1) End-to-end learning of semantic grasping outperforms a traditional detectionclassification-grasping method. (2) Two-stream decompositions of grasp and class prediction outperform a single-stream model. (3) A separated architecture is better at grasping than two-branch architectures. (4) A 9-layer CNN with attention performs comparably to a 16-layer CNN without attention. (5) Auxiliary grasping data improves classification of our two-stream architectures. (6) Auxiliary semantic data that reflects the distribution of objects we are evaluating on improves class grasping accuracy.

Baseline Comparisons
The primary challenge in moving from non-semantic grasping to semantic grasping is that the policy must now extract a specific object from clutter, instead of just reaching for the object that is easiest to grasp. Furthermore, the actions that result in good grasps on a target object may not be the best actions that suggest the object's semantic identity (e.g., a handle on a shovel is easy to grasp but could also be the handle for a toilet brush). Our "Dorsal Only" baseline lacks semantic reasoning capabilities; while it can pick up an object 73.1% and 73.5% of the time on training and test objects, the fraction of successful grasps being of the right class is random (20.8% on training objects and 17.2% on test objects, Table 2, "Dorsal Only") 3 .
We compare our end-to-end semantic grasping policies to a semantic grasping approach in which object detection, classification, and grasp planning are separated into separate pipelines, which we refer to as the "Bounding Box Model". Since our grasping dataset lacks annotated segmentation masks for training region proposal models, we cannot train an object detector or classifier on grasping images; instead, we train an Inception-Resnet-V2 classifier [25] on the S1 dataset (Section 3.3) to obtain a classifier with 95% accuracy on held-out "present images" of training objects and 35% accuracy on "present images" test objects collected during the data-collection phase. For object detection, we adopt a simplified version of RCNN [5,20] that performs a grid search within the initial image I 0 to find the fixed-size crop which most likely corresponds to the target class, as predicted by the trained classifier. The crop bounds in screen space are then projected to the workspace plane using precise camera extrinsic measurements, forming a quad that bounds the detected object. The non-semantic dorsal stream described in Section 3.1 is then instructed to grasp an object constrained to these quad bounds ( Figure 6). The "Bounding Box" model has 41.6% and 32.4% class accuracy among successfully grasped objects, which is a significant improvement over the random behavior of a non-semantic dorsal grasping model. However, this baseline still had lower classification performance than all of our end-to-end learned models (> 51.7%). Failure modes of this model include (1) failing to center the object in the detected crop, causing a poor grasp location on the object and (2) including more than one object in the crop, which may result in the arm picking up another object close to the target. These errors are primarily caused by a large domain shift from canonical image classification datasets to grasping images, despite the classifier being trained on present images of the same objects. Although learning a region-proposal model would result in better canonical crops for grasping, acquiring these labels is a labor-intensive procedure. In contrast, our models exploit closed loop robotic interaction with visual feedback to circumvent the need to acquire manually-annotated region-proposal annotations, and uses learned hand-eye coordination rather than camera calibration.

Two-Stream Models
We trained and evaluated 8 two-stream models, with 2 using a Two-Branch architecture and 6 using Separated architecture. Separated architectures performed better than Two-Branch models that jointly predicted class probability and grasp success (Appendix D, rows 4,5 vs. rows 6,10). Figure 1a visualizes the value function learned by our two stream models. In the left image, the image is superimposed with the actions that correspond to the highest probabilities of picking up any object. On the right, the image is superimposed with actions which are color-coded by the most likely object class at that location. The colors roughly correspond to the locations of objects, showing our model's ability to jointly learn visual grounding of semantic concepts.
A shallower model with soft keypoint attention (Section B.1) yields about the same semantic grasping performance as a deeper convolutional model while using half as many parameters (Appendix D, rows 6,7 vs. rows 8,9). The effect of auxiliary grasping data has similar effects on both attention and non-attention models, suggesting that attention allows us to obtain comparable accuracy with fewer parameters. We use "Separated + Attention" as the base model for further experiments with transfer learning from auxiliary semantic datasets.

Effect of External Dataset Transfer
We first examine the effect of auxiliary grasping data that is used to provide additional supervision to the grasp branch of our Two Branch architecture. In our "Separated" architectures, we train both grasp and class branches from the Two Branch model but only use the latter as the ventral stream of the two-stream model. Even when the grasp branch is replaced by a separately-trained model and not used for grasp prediction at evaluation, the auxiliary data still helps to learn low-level sensorimotor representations that are useful for classifying semantic objects. The Two-Stream Separated + Aux Grasp data yields a 3.0% improvement on Class/Attempts over the Two-Stream Separated model, and the corresponding attention-based model yields a 2.2% improvement on training objects and 3.3% improvement on unseen objects.
Next, we consider the effect of transfer from the S1 and S2 datasets described in Section 3.3. While our models do not seem to benefit from external domain transfer of ImageNet/JFT, our experiments suggest that training on the "optimal" auxiliary dataset does improve classification success: S1 improves classification performance on training objects only, while training jointly on present images of test objects (S2) improves grasping accuracy on both training/test objects. These mixed results on domain transfer suggest that either JFT/ImageNet do not provide much information about our random selection of evaluated training and test objects, or that a more powerful domain adaptation architecture is required to generalize recognition of our object categories from these classification datasets. From a practical standpoint, this transfer is still useful as it enables robots to learn about object semantics without having to grasp them. However, it remains an open question as to how we can perform more effective knowledge transfer from datasets like ImageNet.
Overall, the results in Table 2 point to two-stream architectures being better suited to semantic grasping than single-stream models or a pipelined approach using an object detector. The overall success rate of the methods still leaves room for improvement, which indicates the difficulty of semantic grasping in clutter. This task combines the physical challenge of picking up diverse and varied objects with the finegrained perception challenge of predicting which object will be picked up. We hope that future work will continue to improve on these results, and that end-to-end semantic grasping will emerge as a useful benchmark task for both robotic control and situated object detection.

Discussion
In this work we propose a learning-based system for visual robotic grasping that combines spatial and semantic reasoning into a single neural network policy. Our network is decomposed into two streams, with the "ventral" stream being selective for the commanded object class while the "dorsal" stream selects for good grasps, independent of object identity.
Semantic grasping is not only a compelling research area for studying the intersection of spatial and semantic reasoning, but is also an as-yet unsolved task of considerable practical value. Semantic grasping capabilities could enable a household robot to clean a room, perform industrial sorting tasks on food products, or pack an ever-changing inventory of products into boxes. In each of these situations, the robot is likely to encounter objects that it has never seen or grasped before, and must appropriately generalize both spatial and semantic understanding to handle object diversity.
More broadly, robotics present an unique opportunity for cheaply computing labels for large amounts of complex data via active perception. We address the difficult problem of assigning semantic labels to large amounts of cluttered data via label propagation: during data collection, the robotic arm isolates the grasped object so that an uncluttered view of the object may be captured. This can be done by either lifting up the object to the camera for closer inspection, or using an auxiliary wrist-mounted camera to view the contents of the gripper directly. Such "present images" are centered, uncluttered, and easy to classify, and the temporal coherence of a grasping sequence allows us to retroactively assign the resulting semantic label to an earlier cluttered grasping scene. Our setup enables us to collect additional present poses per grasped object (e.g. multiple orientations), though we leave this for future work. The Amazon Mechanical Turk service was used to label 20,000 present images for label propagation. The raters' choices included the 16 categories in Table 1 as well as No Objects, Multiple, and None of these. We used a simple majority consensus to obtain the final label of each grasp sequence. Each task took under 20 seconds and 2 cents to complete on average.  Figure 8: Two-stream model for semantic grasping. ReLU activations and batch normalization [12] are used at all but the final layer of each model.

B.1 Soft Keypoint Attention
In soft keypoint attention [4], the output a cij of a convolutional layer is normalized across the spatial extent via a softmax s cij = e acij /τ / i j e a ci j /τ . The spatial softmax for each channel is then used as weights for spatially averaging a mesh grid f c = ( i i · s cij , j j · s cij ), which yields expected screen coordinates for salient visual features in that channel. These features are normalized to mean 0 and standard deviation 1 prior to classification, since class identity should be invariant to where the spatial features are located and their scale. Figure 9: Visualization of soft keypoint attention. Each point corresponds to the mean screen coordinates of a spatial softmax over channels in a convolution layer.
Since soft keypoint attention returns a low-dimensional representation of attended features in the form of coordinates, we can directly model the dorsal and ventral branches of our model using a few fully-connected layers 8. Although our models experimentally suggest that soft keypoint attention does not contribute to substantially better performance than simply using a deeper CNN, it is a convenient way to reduce the number of parameters in the model and accelerate training.

C Evaluation Details
Each grasp attempt consists of up to 10 time steps and at each time step, the robot records the current image I t and samples an action a t along which to move the gripper. At each step, action is optimized for several iterations via CEM [24], which optimizes the neural critic described in Section 3.1 with respect to the proposed action. At the final grasp step T , the robot closes the gripper and evaluates the success of the grasp by recording feedback from the gripper. Data collection initially is bootstrapped on top of a policy that issues random motor commands at each time step. Off-policy data was collected T = 10, with 2 iterations of CEM per step. At inference time, we used T = 4 with 3 iterations of CEM per step to speed up the grasping system. The total numbers of grasps vary due to hardware failures that required restarting the evaluation protocol, resulting in some robots/objects having more attempts than others. However, grasps for a model were only scored after all evaluations were completed.
Occasionally, the robot picks up more than one object. If one of the objects grasped was the correct commanded one, we still count the grasp as a success (though the removal of two objects may result in another grasp failing).
Our semantic grasping models are trained using off-policy data from a non-semantic grasping policy, which tends to servo to the nearest object during data collection. Consequently, our dorsal and ventral streams are not trained on actions that are distal from the gripper, and thus cannot accurately reason about semantics and spatial reasoning for distant actions. Instead of the random approach move in [17], we modify the approach trajectory to hover the end-effector over 5 equally-spaced regions of the bin, so that the action search is able to consider moves from a variety of starting points (one of which will be close to the best candidate object). For each "hover pose", we optimize for the action that maximizes the semantic probability p(c|I t , a t , g), and select the maximum-likelihood action among all such poses as the actual approach pose, in order to bring the gripper closer to the object.
Small probabilities from the ventral stream can destabilize predictions from the dorsal stream, even if the value function correctly identifies the right object to grasp. At inference time, we smooth the ventral probability via min(p(c|I t , a t , g) + 0.5, 1), prior to multiplying with the dorsal probability, which allows the robot to pick up an object even if it is unsure about the semantic identity of an object.


extfieldsofstudy: Engineering Mathematics

text:
Interpolation and Extrapolation of Toeplitz Matrices via Optimal Mass Transport

In this work, we propose a novel method for quantifying distances between Toeplitz structured covariance matrices. By exploiting the spectral representation of Toeplitz matrices, the proposed distance measure is defined based on an optimal mass transport problem in the spectral domain. This may then be interpreted in the covariance domain, suggesting a natural way of interpolating and extrapolating Toeplitz matrices, such that the positive semi-definiteness and the Toeplitz structure of these matrices are preserved. The proposed distance measure is also shown to be contractive with respect to both additive and multiplicative noise, and thereby allows for a quantification of the decreased distance between signals when these are corrupted by noise. Finally, we illustrate how this approach can be used for several applications in signal processing. In particular, we consider interpolation and extrapolation of Toeplitz matrices, as well as clustering problems and tracking of slowly varying stochastic processes.

Although many methods rely on stationarity for modeling signals, such assumptions are typically not valid over longer time horizons. Therefore, tools for interpolation and morphing of covariance matrices are important for modeling and fusing information. Several such tools for interpolating covariances have recently been proposed in the literature, for example methods based on g-convexity [9], optimal mass transport [18], and information geometry [19]. An alternative approach for such interpolation is to relax, or "lift", the covariances and instead consider interpolation between the lifted objects. For example, in [20] (see also [21]), interpolation between covariance matrices is induced from the optimal mass transport geodesics between the Gaussian density functions with the corresponding covariances.
The topic of optimal mass transport (see, e.g., [22], [23]) was originally introduced in order to address the problem of, in a cost efficient way, supplying construction sites with building material and has been used in many contexts, such as, e.g., economics and resources allocation. Lately, it has also gained interest in application fields such as image processing [24], [25] signal processing [26]- [28], computer vision and machine learning [29]- [33]. In this work, we will utilize optimal mass transport to model changes in the covariance structure of stochastic processes, or signals. To this end, we propose a new lifting approach, where the lifting is made from the covariance domain to the frequency domain, using the fact that any positive semi-definite Toeplitz matrix has a spectral representation. We combine this approach with the frequency domain metric based on optimal mass transport, proposed in [27], in order to define pairwise distances between Toeplitz matrices. This is done by considering the minimum distance, in the optimal mass transport sense, between the sets of power spectra consistent with each of the Toeplitz matrices. The proposed distance measure is shown to be contractive with respect to additive and multiplicative noise, i.e., it reflects the increased difficulty of discriminating between two stochastic processes if these are corrupted by two realizations of a noise process. Also, we show that the proposed distance measure gives rise to a natural way of interpolating and extrapolating Toeplitz matrices. The interpolation method preserves the Toeplitz structure, the positive semi-definiteness, as well as the zeroth moment of the interpolating/extrapolating matrices.
The proposed optimal mass transport problem is in its original form an infinite-dimensional problem. As an alternative to finding solutions using approximations based on discretizations of the underlying space, we show that certain formulations of the problem allows for approximations by a semi-definite program using a sum-of-squares representation. Also, we illustrate how the method can be used for interpolation, extrapolation, tracking, and clustering. This paper is organized as follows. In Section II, we provide a brief background on the moment problem, i.e., determining the power spectrum from a partial covariance sequence or finite covariance matrix, as well as introduce the problem of optimal mass transport. Section III introduces the proposed distance notion for positive semi-definite Toeplitz covariance matrices. Here, the dual problem is derived, and properties of the proposed distance notion are described. In Section IV, we describe applications of the proposed distance notion, such as induced interpolation, extrapolation, tracking, and clustering. Section V formulates a sum-of-squares relaxation of the dual problem. Section VI provides numerical illustrations of the proposed distance notion, as well as the described applications. Finally, Section VII concludes upon the work.

Notation
Let M n denote the set of Hermitian n × n matrices and let (·) T denote the transpose, (·) H the Hermitian transpose, and (·) the complex conjugate. Let T = (−π, π] and let C perio (T) denote the set of continuous and 2π-periodic functions on T. The set of linear bounded functionals on C perio (T), or equivalently, the dual space of C perio (T), is the set of generalized integrable functions on the set T, here denoted by M(T). Thus, M(T) includes, e.g., functions containing singular parts such as Dirac delta functions [34]. 1 Further, we let M + (T) denote the subset of such functions that are non-negative. We use Φ, f to denote the application of the functional Φ on f , e.g., if f ∈ C perio (T) and Φ ∈ M(T). For Hilbert spaces, ·, · is the standard inner product, e.g., when X and Y are vectors or matrices, then X, Y = trace(XY H ). We denote matrices by boldface upper-case letters, such as X, whereas vectors are denoted by boldface lower-case letters, such as x.

A. Stochastic processes and spectral representations
We will in this work consider complex-valued discrete time stochastic processes, or signals, y(t) for t ∈ Z. These will be assumed to be zero mean and wide sense stationary (WSS), i.e., E(y(t)) = 0 for all t ∈ Z, and the covariance being independent of t. Here, E(·) denotes the expectation operator. The frequency content of the process y(t) may then be represented by the power spectrum, Φ, i.e., the non-negative function on T whose Fourier coefficients coincide with the covariances: for k ∈ Z (see, e.g., [36,Chapter 2]). Typically in spectral estimation, one considers the inverse problem of recovering the power spectrum Φ from a given set of covariances r k , for k ∈ Z, with |k| ≤ n − 1. The condition for any such reconstruction to be valid is that Φ should be consistent with the covariance sequence {r k } |k|≤n−1 , i.e., (1) should hold for |k| ≤ n − 1. The corresponding n × n covariance matrix is which is a Hermitian Toeplitz matrix, since y(t) is WSS. Thus, expressed in the form of matrices, a spectrum is consistent with an observed partial covariance sequence, or, equivalently, a fi- and is the Fourier vector. Note that Γ(Φ) is a Toeplitz matrix, since a(θ)a(θ) H is Toeplitz for any θ. It may be noted that for any positive semi-definite Toeplitz matrix, R, there always exists at least one consistent power spectrum; in fact, if R is positive definite, there is an infinite family of consistent power spectra [37]. It may be noted that for singular Toeplitz covariance matrices, the spectral representation is unique. This fact has recently been successfully utilized in atomic norm minimization problems for grid-less compressed sensing of sinusoidal signals [38]. In this work, we are mainly interested in the non-singular case, where several power spectra are consistent with given covariance matrices. In Section III, we will utilize such spectral representations in order to define a notion of distance between pairs of Toeplitz matrices. This distance will be defined in terms of the minimum optimal mass transport cost between the sets of power spectra consistent with the matrices.

B. Optimal mass transport
The Monge-Kantorovich problem of optimal mass transport is the problem of finding an optimal transport plan between two given mass distributions [22], [23]. The cost of moving a unit mass is defined on the underlying space, and the optimal transport plan is defined as the plan minimizing the total cost. The resulting minimal cost, associated with the optimal transport plan, can then be used as a measure of similarity, or distance, between the two mass distributions. The idea of utilizing the optimal mass transport cost as a distance measure has been used, e.g., for defining metrics on the space of power spectra [27], whereas the optimal transport plan has been used for tracking stochastic processes with smoothly varying spectral content and for spectral morphing for speech signals [28]. Recently, the interpretation of the optimal transport plan as providing an optimal association between elements in two mass distributions has been used as a means of clustering in fundamental frequency estimation algorithms [39].
As in [27], we consider the following distance between two spectra Φ 0 and Φ 1 : where T 2 = T × T denotes the 2-D frequency space. Here, the cost function, c(θ, ϕ), details the cost of moving one unit of mass between the frequencies θ and ϕ. The transport plan, M (θ, ϕ), specifies the amount of mass moved from frequency θ to frequency ϕ. The objective in (5a) is the total cost associated with the transport plan M and the constraints (5b) and (5c) ensure that M is a valid transport plan from Φ 0 to Φ 1 , i.e., the integration marginals of M coincide with the spectra Φ 0 and Φ 1 . It may be noted that, due to these marginal constraints, the distance measure S is only defined for spectra of the same mass, or total power. However, S may be generalized in order to allow for mass differences by including a cost for adding and subtracting mass by postulating that the spectra Φ 0 and Φ 1 are perturbations of functions Ψ 0 and Ψ 1 that have equal mass. As in [27], this may be formulated as where κ > 0 is a used-defined parameter detailing the cost of adding or subtracting mass. Also, if the cost function in the optimal mass transport problem in (5) is chosen as c(θ, φ) = d(θ, φ) p , p ≥ 1 for any metric d(θ, φ) on T, then W (Φ 0 , Φ 1 ) = S(Φ 0 , Φ 1 ) 1/p is the so-called Wasserstein metric on M + (T). Similarly, for S κ , the following theorem holds.
Consider a situation where we need to discriminate between two signals on the basis of their statistics or of their power spectra. In such cases, additive or multiplicative noise typically impede our ability to differentiate between the two. This was considered in [27] and it was shown that the transportation distance respects this property in the sense that corrupting two signals with additive and (normalized) multiplicative noise decreases their transportation distance. Specifically, the following theorem holds.

Theorem 2 ([27]
). Let p ≥ 1, κ > 0, and let the cost function be c(θ, φ) = |θ − φ| p , and let W κ (Φ 0 , Φ 1 ) be defined by (7) Then, W κ (Φ 0 , Φ 1 ) is contractive with respect to the additive and normalized multiplicative noise, i.e., with respect to As we shall see, it is possible to construct notions of distance on the space of positive semi-definite Toeplitz matrices that have properties similar to those stated in Theorem 2.

III. A NOTION OF DISTANCE FOR TOEPLITZ MATRICES
As noted above, any positive semi-definite Toeplitz matrix R has at least one spectral representation, i.e., there exists at least one spectrum Φ that is consistent with it. Thus, we define the distance, T , between two positive semi-definite Toeplitz matrices, R 0 and R 1 , as the minimum transportation cost, as measured by S, between spectra consistent with the respective matrices, i.e., Considering the definition of S in (5), this can equivalently be formulated as the convex optimization problem Note that the formulation in (9) is only defined for covariance matrices with the same zeroth covariances, or, equivalently, covariance matrices whose consistent spectra have the same mass. However, in order to allow for mass differences, T (R 0 , R 1 ) can be generalized in analog with (6) as or, equivalently, Typically, the cost function c would be selected to be symmetric in its arguments, in which case T κ would be symmetric as well, i.e., T κ (R 0 , R 1 ) = T κ (R 1 , R 0 ) for any Toeplitz matrices R 0 , R 1 ∈ M + . Although many possible choices of such cost functions exist, we will in most examples presented herein consider c(θ, ϕ) = e iθ − e iϕ 2 , i.e., the cost function quantifies distances as the square of the distance between the corresponding points on the unit circle. As we show in Section V, this particular choice of cost function allows for a sum-ofsquares relaxation of the dual formulation of (11). This dual formulation is presented next.

A. The dual formulation
In order to study properties of the distance notion T κ , we consider the dual formulation of (11), where we assume that the cost function, c, is non-negative on T 2 . In order to address this problem, we first note that the adjoint operator 2 where, in the first line, the inner product is the one associated with M n , and the second line is the bilinear form with arguments Γ * (R) ∈ C perio (T) and Φ ∈ M + (T). With this result, we can derive an expression of the dual problem by considering the Lagrangian relaxation of (11). The Lagrangian is given by Note that the Lagrange multiplier matrices Λ 0 and Λ 1 may be taken as Hermitian matrices, as R 0 − Γ (Ψ 0 ) and R 1 − Γ (Ψ 1 ) are Hermitian, and thus all inner products are real. Considering the infimum of L κ with respect to Ψ 0 and Ψ 1 , it may be noted that this is only finite if Γ * (Λ 0 )(θ) ≤ κ and Γ * (Λ 1 )(ϕ) ≤ κ, for all θ, ϕ ∈ T. If this is satisfied, the Lagrangian is, for any fixed non-negative M , minimized by Ψ 0 and Ψ 1 given by Strictly speaking, this is the pre-adjoint operator of Γ. where where we have used that the cost function c is non-negative. This yields the dual problem.
Proposition 1. For any cost function c(θ, ϕ), such that c(θ, ϕ) ≥ 0, for all θ, ϕ ∈ T, the dual problem of (11) is where Since the primal and dual problems are convex and the set of feasible points Ω κ c has non-empty interior for any κ ≥ 0, Slater's condition (see, e.g., [40]) gives that strong duality holds and hence the duality gap between (11) and (12) is zero. Also, note the strong resemblance in form compared to the dual of (6), which is given by [27] maximize Similarly, for the case when the zeroth covariances of R 0 and R 1 are required to be equal, the dual problem is where This is the dual problem of (8), where adding and subtracting mass is not allowed in the transport problem.

B. Properties of distance notion T κ
For the distance notion T κ in (10), the following proposition holds.
Proposition 2. Let the cost function c be non-negative and such that c(θ, θ) = 0, for all θ ∈ T. Then, the distance measure T κ defined in (11) is contractive with respect to additive noise. If c is also shift-invariant, i.e., c(θ − φ, ϕ − φ) = c(θ, ϕ), for all θ, ϕ, φ ∈ T, then T κ is also contractive with respect to multiplicative noise whose zeroth covariance is smaller than or equal to unity.
The statement of this proposition is that when two stochastic processes become contaminated by noise, the distance notion T κ decreases and hence the processes become harder to distinguish. Intuitively, this is a desirable property of T κ , as additive or multiplicative noise should indeed impede ones ability to discriminate between two processes. Proposition 2 may be proven by utilizing results in [27]. However, in the interest of making the exposition self-contained, we will here provide a direct proof based on the dual formulation in (12). Here, we prove the first part of the proposition, i.e., the case of additive noise. For the case of multiplicative noise, the proof is deferred to Appendix A.
Proof. Consider two processes with covariance matrices R 0 and R 1 and assume that they are both additively corrupted by an independent noise process with covariance R w . This results in processes with covariances R 0 = R 0 + R w and R 1 = R 1 + R w , respectively. From the dual formulation in (12), the distance between these covariance matrices is given by Let Φ w be any spectrum consistent with the noise covariance R w , i.e., with Γ(Φ w ) = R w . We then have that Hence, it follows that It may be noted that the assumptions made in Proposition 2 regarding the cost function c are quite mild, allowing for a large class of potential cost functions. In particular, c(θ, ϕ) = e iθ − e iϕ 2 satisfies the conditions of Proposition 2 for both additive and multiplicative noise.

IV. INTERPOLATION, EXTRAPOLATION, AND TRACKING
The formulation in (8) does not only define a notion of distance between two Toeplitz covariance matrices, R 0 and R 1 ; it also provides a means of forming interpolating matrices, i.e., defining intermediate covariance matrices R τ , for τ ∈ (0, 1). In order to define the interpolating matrices, we again utilize the spectral representation of positive semi-definite Toeplitz matrices. To this end, we note that, given an optimal transport plan, M , found as the functional minimizing (8), one may define spectra intermediate to the marginals T M (θ, ϕ)dϕ and T M (θ, ϕ)dθ by linearly shifting the frequency locations of the spectral mass, as dictated by M . That is, any mass transferred from φ to φ + ϕ is defined to, at τ ∈ [0, 1], be located at frequency φ + τ ϕ. Using this, the intermediate spectrum is given by Here, δ θ is the Dirac delta function localized at θ, i.e., δ θ (φ) δ(φ − θ) and the integrands are extended periodically with period 2π outside the domain of integration. Also, we denote by {x} mod T the value in T that is congruent with x modulo 2π. Based on this definition of the interpolating spectrum Φ M τ , the corresponding interpolating covariance matrix, R τ is defined as the unique Toeplitz matrix consistent with this spectrum, i.e., To simplify the following exposition, let I τ (M ) Γ(Φ M τ ) denote the linear operator in (15) that maps a transport plan to an interpolating covariance matrix, i.e., R τ = I τ (M ). It may be noted from (15) that I τ (M ) is welldefined also for τ = [0, 1], i.e., the formulation allows also for extrapolation. The following proposition follows directly from the definition in (15).
Proposition 3. For any τ ∈ R, the following basic properties hold for R τ = I τ (M ): a) If R 0 and R 1 have the same zeroth covariance, then it is also the zeroth covariance of R τ . b) The matrix R τ is a Toeplitz matrix. c) The matrix R τ is positive semi-definite.
Due to these properties, the proposed method offers a way of interpolating the covariances of, e.g., slowly varying time series as the interpolant R τ allows for modeling linear changes in the spectrum of the process. Remark 1. The interpolation approach generalizes trivially to the formulation in (11) between the covariances Γ(Φ M 0 ) and Γ(Φ M 1 ) under the assumption that Ψ 0 and Ψ 1 are perturbations of Φ M 0 and Φ M 1 , respectively. In order to define an interpolation and extrapolation procedure from R 0 to R 1 where there is a cost κ for adding and subtracting mass, one may, along with the interpolation path I τ (M ), linearly add the part corresponding to the added and subtracted mass, i.e., Note that in this scenario, positive semi-definiteness cannot be guaranteed for the extrapolation case.

A. Comparison with other methods
The properties in Proposition 3 distinguish the proposed interpolant R τ from other proposed matrix geodesics. As an example, consider the basic method of forming interpolants using convex combinations of R 0 and R 1 , i.e., R conv τ = τ R 0 + (1 − τ )R 1 , for τ ∈ [0, 1]. This preserves the Toeplitz structure, as well as the zeroth covariance and the positive semi-definiteness. However, from a spectral representation point of view, the convex combination gives rise to fade-in fade-out effects, i.e., only spectral modes directly related to R 0 and R 1 can be represented, and there can be no shift in the location of these modes (see also Example VI-B and Figure 3). Other more sophisticated options include, e.g., the geodesic with respect to g-convexity [9] and the geodesic in [20], [21], which builds on optimal mass transport of Gaussian distributions and can be expressed as Although both of these geodesics preserve positive semidefiniteness, they neither preserve the Toeplitz structure nor the zeroth covariance. As noted above, the three properties in Proposition 3 hold for any τ ∈ R for the proposed approach, and thus directly allows for extrapolating using (15). In contrast, it may be noted that for the linear combination R conv τ there are no guarantees that the resulting matrix is positive semi-definite if τ / ∈ [0, 1]. Also, note that the alternative geodesics,R τ andȒ τ , do not naturally generalize to extrapolation.

B. Tracking of slowly varying processes
The proposed interpolant R τ = I τ (M ) may also be readily used for tracking slowly varying stochastic process. As noted above, I τ (M ) allows for the modeling of slow, i.e., locally linear, shifts in the location of spectral power. Building on this property, we can extend the optimal transport problem in (8), in order to fit a covariance path R τ to a sequence of J covariance matrix estimates,R τj , for j = 1, . . . , J. As R τ is unambiguously determined from a transport plan M via I τ (M ), this tracking problem may be formulated as the convex optimization problem where λ > 0 is a user-specified regularization parameter. As may be noted from (19), the optimal transport plan M is here determined as the one that minimizes not only the transport cost, but also takes into account the deviations of the interpolant R τj = I τj (M ) from the available covariance matrix estimatesR τj . The behavior of this construction is illustrated in the numerical section.

C. Clustering: Barycenter computation
As a further example, we will see that the barycenter with respect to the distance notion T κ may be formulated as a convex optimization problem. This might be desirable in clustering or classification applications, where one is interested in either identifying classes of signals or processes based on their covariance matrices or associate a given covariance matrix with such a signal class. Considering the case of clustering, assume that L covariance matrices, R , = 1, . . . , L, are available. Then, we may define their barycenter via T κ according to i.e., as the covariance matrix that minimizes the sum of T κ for the set of covariance matrices R . Explicitly, R bary solves the convex optimization problem This formulation allows for using, e.g., K-means clustering (see, e.g., [41]) in order to identify classes of covariance matrices, as well as classify a given covariance matrix according to these classes. Classification of a covariance matrix R according to classes defined by a set of barycenters R In Section VI-D, we present a simple illustration of this potential application, considering unsupervised clustering of phonemes.
V. SUM-OF-SQUARES RELAXATION In order to solve the dual problem in (12) in practice, it has to be implemented as a finite dimensional problem, e.g., by gridding the space T 2 and thereby approximating the set Ω κ c using a finite number of constraints. However, for the special case of c(θ, ϕ) = e iθ − e iϕ 2 , this can also be done using a sum-of-squares (SOS) relaxation (see, e.g., [42]) of (12). For simplicity of notation, we present only the case without κ (κ = ∞), i.e., the dual of (8) as formulated in (13). To this end, identify z = e iθ , w = e iϕ which allows us to write Thus, the set of constraints defining the feasible set Ω c is given by Γ * (Λ 0 )(θ) + Γ * (Λ 1 )(ϕ) ≤ c(θ, ϕ), for all θ, ϕ ∈ T, or, equivalently, 23) Note that in the two-dimensional trigonometric polynomial (23), the coefficient for z −k1 w −k2 is equal to T 1 for k 1 = 0 and k 2 ∈ Z n\0 −1 for k 1 = 1 and k 2 = −1 −1 for k 1 = −1 and k 2 = 1 0 otherwise where Z n\0 = {k ∈ Z | |k| ≤ n, k = 0}. Here, diag (X, k) denotes the column vector containing the elements on the kth super-diagonal of the matrix X, if k > 0, and the elements on the kth sub-diagonal, if k < 0, with 1 denoting a column vector of ones of appropriate dimension.
In order to formulate a computationally feasible optimization problem, we remove the non-negativity constraint for the two-dimensional polynomial in (23) and instead impose that it should have a sum-of-squares representation [42], [43]. In particular, we impose that the polynomial in (23) should be of the form where Q ∈ M m 2 is positive semi-definite, and ⊗ is the Kronecker product. Note that any polynomial on this form is non-negative by definition, and, furthermore, for any non-negative polynomal P * , there is a sequence of polynomials, P m , on the form (24) such that P m − P * ∞ → 0 as m → ∞.
Next, note that the coefficients of P are associated to the elements of Q in (24) according to and T k is the matrix with ones on the kth diagonal and zeros elsewhere. Putting these facts together by requiring that the polynomial (23) can be written as a sum-of-square (24), we subject to Q 0 It is worth noting that the P defined by the optimal Q in (24) is non-negative (actually, strictly positive) on the unit torus. Thus, for the solution of the SDP problem in (25), there will be a corresponding feasible solution of (13). Therefore, any solution to (25) will give a lower bound for the optimal objective value of (13). However, any non-negative polynomial may be arbitrarily well approximated by P of the form (24) by a suitable choice the degree m, and thus the maximal objective value of (25) will converge to the maximal objective value of (13) as m grows. By comparison, directly discretizing T, thereby approximating Ω c using a finite number of constraints, yields an upper bound for the maximal objective value of (13). In this case, the maximal objective value will also converge to that of (13), this time from above, as the spacing of the discretization of the grid becomes finer.
VI. NUMERICAL EXAMPLES In this section, we present some numerical examples illustrating different aspects and application areas of the proposed distance notion T κ , as well as the interpolant R τ = I τ (M ).

A. Trajectory example
In order to illustrate the behavior of the proposed interpolation method, we consider a simple scenario with covariance matrices of dimension n = 2. Consider covariance matrices of the form where r ∈ [−1, 1]. Thus, considering interpolating paths R τ , these will be on the form Figure 1 displays the real and imaginary part of r 1,τ for τ ∈ [0, 1] when varying the magnitude r of the off-diagonal element of R 0 and R 1 between 0 and 1. As can be seen, the trajectories of r 1,τ approximately correspond to convex combinations (1 − τ )r + τ re i 5π 6 when r is close to zero, whereas they are considerably curved for r closer to 1. It may be noted that for the singular case, i.e., r = 1, the trajectory of r 1,τ coincides with the unit circle. Thus, we see that for covariance matrices that have consistent spectra that are essentially flat, the interpolant R τ will be approximately equal to the convex combination (1−τ )R 0 +τ R 1 . At the other extreme, the interpolant corresponding to covariance matrices R 0 and R 1 that have consistent spectra that are close to being singular will also have almost singular consistent spectra.

B. Interpolation and extrapolation for DOA
Next, we illustrate the proposed methods ability to produce interpolants R τ that are consistent with locally linear changes in the frequency location of spectral power. This is illustrated using a direction-of-arrival (DOA) estimation problem. Consider a uniform linear array (ULA) with 15 sensors with half-wavelength sensor spacing and a scenario where two covariance matrices  (9), as R 0 and R 1 have the same zeroth covariance, in order to find the optimal transport map M . Then, using (15), we compute covariance matrices R τ , for τ ∈ [0, 2], i.e., we both interpolate on τ ∈ [0, 1] and extrapolate on τ ∈ (1,2]. This is then compared to the basic interpolant R conv τ based on convex and linear combinations of R 0 and R 1 , as well as the more sophisticated covariance matrix geodesicsR τ andȒ τ , as defined in (17) and (18), respectively. For these four cases, we then estimate the corresponding inter-and extrapolated spectra using the correlogram, i.e., as Φ corr (X, θ) = a(θ) H Xa(θ) (28) where X is substituted for the four different covariance interpolants. The results for the proposed interpolant R τ are shown in Figure 2. As can be seen, R τ indeed models a scenario where one of the targets has a constant location, whereas the second target splits upp into two smaller targets. Note also that R τ implies that the smaller target continue linearly with respect to the look-angle, θ, also for the extrapolation case, i.e., for τ > 1. In contrast, the convex (linear for τ > 1) combination R conv , as shown in Figure 3, display undesirable behavior; clear fade-in fade-out effects are visible, and non-negativity is violated as τ → 2 due to R conv becoming indefinite. Similar objections may be raised against the geodesicsR τ andȒ τ , shown in Figures 4 and 5, both displaying fade-in fade-out effects and thereby fail to model any displacement of the targets. Also, note that the total power of the signal varies greatly as τ goes from 0 to 2, especially forȒ τ .  (19), fitted to a sequence of five covariance estimates.

C. Tracking of an AR-process
Next, we illustrate the approach in (19) for the tracking of signals with slowly varying spectral content. To this end, consider a complex autoregressive (AR) process with one complex, time-varying pole. The pole is placed at a constant radius of 0.9, and moves from the frequency 0.4π to 0.6π. Spectral estimates based directly on covariance matrix estimatesR are shown in Figure 6. These covariance matrix estimates are obtained as the outer product estimate, based on 100 samples each, where the overlap between each estimate is 80 samples. As can be seen, the spectral estimates are very noisy and vary greatly in power. Using five of these covariance matrix estimatesR, evenly spaced throughout the signal, we solve (19) in order to obtain an estimated covariance path, R τ . The resulting spectra, estimated using (28), are shown in Figure 7. As can be seen, the path resulting from the proposed method allows for a smooth tracking of the shift in spectral content.

D. K-means clustering
As a simple illustration of how to utilize the barycenter formulation in Section IV-C, we consider the application of unsupervised clustering of phonemes. Specifically, we consider 7 utterances; 3 utterances of the phoneme /ae/, 2 utterances of /oy/, and 2 utterances of /n/. For each utterance, we estimate a covariance matrix of dimension n = 10 and then run a K-means algorithm that alternates between classifying each estimated covariance matrix according to (22) and computing new barycenters according to (20). We initiate the algorithm by choosing initial barycenters as convex combinations of the available covariance matrix estimates, and demand separation into three clusters. The algorithm is then run until convergence, i.e., until the classification has stabilized. The results are shown in Table I. As can be seen, the K-means algorithm utilizing (20) and (22) is indeed able to correctly identify three clusters corresponding to the three different phonemes. Table II presents the corresponding distance matrix, i.e., the matrix detailing the distance T κ between each covariance  Fig. 8. Interpolated spatial spectrum estimated as a(θ) H Rτ a(θ), where Rτ is obtained by solving (9), for the case of one moving source and one static interferer, using with the cost function c(θ, ϕ) = e iθ − e iϕ 2 . matrix estimate and each barycenter. As can be seen, the clusters are quite well-separated.

E. Fixed cost for moving mass
In some scenarios, there might be a relatively large noise component present in both R 0 and R 1 , where the noise power is localized in frequency and this localization is the same in both R 0 and R 1 . This might be the case in, e.g., DOA estimation scenarios where a source of interest is moving in the presence of a stationary interferer. Then, if the source for example moves past the location of the interferer, the optimal transport problem in (9) may result in a power association such that the source and the interferer are mixed together. Such a scenario is illustrated in Figure 8, showing estimated spectra obtained from the interpolant R τ , as given by (9) with the cost function c(θ, ϕ) = e iθ − e iϕ 2 . Here, the source is moving from look-angle θ = 30 • to θ = −20 • , whereas there is a fixed interferer located at θ = 0 • having a third of the power of the source. In order to avoid this type of problem and promoting transport plans that avoid transporting stationary masses, the cost function may be modified according to Thus, the cost function is here formulated such that there is a fixed, baseline cost of moving any mass. It may be noted that this cost function satisfies the assumptions of Proposition 2,   Fig. 9. Interpolated spatial spectrum estimated as a(θ) H Rτ a(θ), where Rτ is obtained by solving (9), for the case of one moving source and one static interferer, using with the modified cost function in (29). and thus the proposition holds also for this choice of c(θ, ϕ). The resulting estimated spectra obtained from the interpolant R τ resulting from solving (9) using this modified cost function is shown in Figure 9. As can be seen, the source and the interferer are now well separated throughout the interpolated path.

VII. CONCLUSIONS AND FURTHER DIRECTIONS
In this work, we have proposed a notion of distance for positive semi-definite Toeplitz matrices. By considering spectral representations of such matrices, the proposed measure is based on distances, in an optimal mass transport sense, between families of spectra consistent with the Toeplitz matrices. We have shown that the proposed distance measure, under some mild assumptions, is contractive with respect to additive and multiplicative noise. The proposed measure may be used to, for example, define inter-and extrapolation of Toeplitz matrices, being of interest in applications such as tracking of slowly varying signals.
A future direction of this methodology is to generalize the distance measure for structured matrices corresponding to input-to-state covariances in the THREE framework for spectral analysis [12]. This generalizes ideas from beamspace processing, enabling the user to improve resolution and robustness in power spectral estimation over selected frequency bands [44], [45]. This will be the subject of further research.

A. Proof of Proposition 2 for multiplicative noise
Here, we present the proof for T κ being contractive with respect to multiplicative noise.
Proof. Let the noise covariance matrix be R w , implying that the covariances of the contaminated processes are R 0 = R 0 R w and R 1 = R 1 R w , respectively. Let Φ w be any spectrum consistent with R w , i.e., Γ(Φ w ) = R w . Also, let the zeroth covariance of the noise process be smaller than or equal to 1, so that