Vehicle Tracking on Satellite Video Based on Historical Model

Vehicle tracking on satellite videos poses a challenge for the existing object tracking algorithms due to the few features, object occlusion, and similar objects appearance. To improve the performance of the object tracking algorithm, a historical-model-based tracker intended for satellite videos is proposed in this study. It updates the tracker by using the historical model of each frame in the video, which contains plenty of object information and background information, so as to improve tracking ability on few-feature objects. Furthermore, a historical model evaluation scheme is designed to obtain reliable historical models, which ensures that the tracker is sensitive to the object in the current frame, thus avoiding the impact caused by changes in object appearance and background. Besides, to solve the drift issue of the tracker caused by object occlusion and the appearance of similar objects, an antidrift tracker correction scheme is proposed as well. According to the comparative experiments conducted on satellite videos dataset SatSOT, our tracker produces an excellent performance. Moreover, sensitivity analysis, varying criteria comparative experiments, and ablation experiments are conducted to demonstrate that the proposed schemes are effective in improving the precision and success rate of the tracker.

the computer. Due to the diversity of moving objects and the video equipment, it is necessary to adapt the object tracking algorithm to different data sources. For example, the ordinary video [1]- [3], the unmanned aerial vehicle (UAV) video [4], [5], the thermal infrared video [6], the synthetic aperture radar video [7], [8], and the satellite video [9]- [21], which is an emerging type of space-based video data in recent years. The satellite video has demonstrated such advantages as wide shooting range, high resolution, and the capability of continuously monitoring the target objects on land.
The object tracking on satellite videos has a wide range of applications in national defense, environmental protection, disaster prevention, and traffic monitoring [12], [22]. Our study focuses on the vehicle tracking, which is vital for traffic and military spying. Besides, due to the wide shooting range of satellite videos, the long-distance tracking is achievable, which is conducive to analyzing the driving motives and trajectories.
There are some issues when the vehicles are tracked on satellite videos. The first one is a small number of features and textures due to the small size of the object. Then, the vehicle is often occluded by trees, bridges, and other obstacles. Besides, similar objects often appear around the tracked object. Thus, it is a challenge to apply the existing object tracking algorithms to satellite videos.
According to the principle followed by the algorithm, the existing object tracking algorithms can be divided into two categories: generative models [23]- [26] and discriminative models [27]- [33]. Generative models focus on the characteristics of the object itself and track the object by conducting iterative search for the similar object frame by frame. Some classical tracking algorithms, such as MeanShift [23], CAMShift [24], particle filter [25], and optical flow [26] are classed as generative models. Differently, discriminative models focus on the difference between the object and the background. It takes background information into consideration, e.g., MOSSE [27], CSK [28], KCF [29], and STC [30]. With the widespread application of deep learning, the object tracking algorithms combined with deep learning have been gradually developed. Some algorithms extract deep convolutional features, e.g., C-COT [34] and ECO [35], both of which use the VGGNet [44], whereas others use end-to-end object tracking methods, e.g., CFNet [45], SiamFC [46], and MDNet [47]. Given the lack of background information, the generative algorithm only models the object itself, which is unreliable for the few-feature objects on satellite videos. By contrast, the discriminative algorithm introduces background information, which makes it advantageous over the generative algorithm. In particular, the correlation filter (CF), which is one of the discriminative algorithms, is widely used in object tracking on satellite videos due to its high speed and high precision. Although the algorithm based on deep learning is equally effective, it is disadvantaged by heavy computational workload and slow speed. In some cases, the few and weak features of objects make it inappropriate for neural networks to learn. In general, CF-based algorithms are superior to generative algorithms in accuracy and are superior to deep-learning-based algorithms in speed, which makes them have greater advantages on the satellite videos with larger image sizes and fewer features. For this reason, the CF is applied in this study.
For most of the current CF algorithms, the focus is on ridge regression. By learning the features of the region of interest (ROI), they can classify the object and background to locate the object. When the filter is updated, most algorithms use a single ROI in the current frame and ignore the historical information of the previous frames, which leads to the waste of information. Nevertheless, historical information contains a great deal of the object information and background information required for modeling the object. Especially, background information could be inputted into CF as negative samples, which could improve the tracking performance for the few-features object. For this reason, a historical-model-based tracker for satellite (HMTS) videos is proposed in this study, as illustrated in the flow diagram (see Fig. 1). Historical models (HMs) are applied to update the tracker. However, not all HMs are suitable due to changes in object appearance and background during the tracking process. Hence, a scheme is used to find the HMs that are best suited to the current frame. In addition, since the vehicle is often occluded by obstacles or there are similar objects appearing around the vehicle, a strategy is adopted to detect vehicle state and prevent the drift of the tracker. The peak and kurtosis (PK) of the response map of the CF [50] are employed to monitor the vehicle state. Additionally, Kalman filter (KF) is adopted to correct the tracking result to avoid drift. To sum up, the contributions of this study are as follows: 1) In this article, we designed a novel HMTS videos to improve the tracking performance on small-size objects, thus addressing the trace failure caused by a lack of features.
It improves the acquisition of object information by utilizing historical information, ensures the reliability of tracking, and reduces historical information waste. 2) We proposed an HM evaluation scheme to avoid the impact caused by changes in object appearance and background. The scheme relies on cross-correlation functions to measure the similarity between the response map and the ideal regression targets. HMs are scored by this scheme to identify the best-suited ones for the current frame. 3) We constructed a tracker correction scheme to prevent the tracker drift caused by object occlusion and similar objects appearance. The scheme uses KF to estimate the motion of the object and PK for the response map evaluation. Based on PK, the trajectory of the tracker is corrected by the prediction function of KF. The rest of this article is organized as follows. In Section II, the relevant work is summarized. Section III presents the design of the proposed tracker HMTS. In Section IV, the experiments are detailed, including configurations, experimental design, experimental analysis, and experimental results. Finally, Section V concludes this article.

A. CF-Based Object Tracking
The emergence of CF marks a key milestone in the development of object tracking algorithm. By introducing the Fourier transform method to convert the originally complex matrix operations into calculations in the frequency domain, it accelerates the calculation process. At the same time, CF increases the use of background information, which improves the accuracy of tracking. In the earliest CF algorithm MOSSE [27], the least squares method is used to distinguish between the background and the object. Then, kernel CF is used by CSK [28] and KCF [29] to introduce circulant matrix, which increases the number of samples, thus improving the accuracy of training. Following the KCF, scholars have successively improved the KCF framework to develop many algorithms based on CF. There are three main approaches to improving the CF. The first one is to apply a suitable feature extraction algorithm [34]- [37]. In [36], color names (CNs) are used, whereas in [34], deep convolutional features are used. The second one is to adapt to the object scale change by introducing a scale-adaptive mechanism [38], [39]. The mechanism is effective in increasing the success rate of the tracker. The last one is to improve the structure of the filter [40]- [43], e.g., BACF [41] adds the spatial regularization term, and AutoTrack [43] adds the temporal regularization term.

B. Object Tracking on Satellite Video
The existing object tracking algorithms have been widely used for the tracking of objects on satellite videos. The traditional generative model is applied in [9] and [10]. Besides, CF algorithm is used in [11]- [17]. Moreover, deep learning is adopted in [18]- [21]. At present, traditional generative model is rarely adopted for satellite videos, which is attributed mainly to its low accuracy and slow speed. By contrast, the discriminative model has been more commonly used. A major solution is to improve CF, which can be achieved from three perspectives. The first one is to change the feature extraction algorithm, the second one is to use motion model, and the last one is to introduce tracker state monitoring scheme. To change the feature extraction algorithm, Shao et al. [12] combined LK optical flow method and histogram of oriented gradient (HOG) feature extraction algorithm to extract features, whereas Wu et al. [15] combined Hu CF and median filter to extract features. Due to the rotation invariance of Hu invariant moment, this algorithm performs well in tracking rotation objects. In [17], Gabor filter is employed to extract features, which leads to an excellent performance for the objects with textures. In terms of motion models, KF is frequently used to predict the position of the moving object [14]- [16], which improves the robustness of the algorithm. As for the tracker state monitoring scheme, Wang et al. [17] proposed the tracking status monitoring indicators (TCMIs) based on the Bayesian framework. If the object is occluded, TCMI can be relied on to guide the tracker, with updating terminated to prevent model drift. In addition, Xuan et al. [16] used the peak value of response map to monitor tracker states and guide the updating of the tracker.

C. Benchmark Dataset Based on Satellite Video
In addition to the improvement of algorithm, the benchmark dataset based on satellite video has also been proposed to make the verification of algorithms more reliable. Yin et al. [48] constructed a large-scale satellite video dataset with a wide range of annotations. At the same time, a benchmark was proposed to evaluate algorithms for their performance, e.g., the multiobject tracking algorithm and the single-object tracking algorithm on satellite videos. Zhao et al. [49] put forward another densely annotated satellite video dataset, which is purposed to evaluate the single-object tracking algorithm.

III. PROPOSED METHOD
In this section, it will be explained how to apply the HM, cross-correlation function, PK, and KF for CF-based object tracking. First, a tracker HMTS based on history model is proposed to track the object. Then, an HM evaluation scheme is introduced to validate each HM, so as to choose the HMs that are sensitive to the object in the current frame. Finally, a scheme is introduced to correct the tracker by using KF and PK.

A. Overall Architecture
In our tracker, kernel CF (KCF) [29] is taken as the baseline. Suppose there is a vectorized feature map X, X ∈ R D . KCF solves the following ridge regression problem: where w represents a filter; X(i) refers to the ith circular shift of X; y(i) indicates the regression target for X(i); λw 2 2 denotes a shrinkage penalty, which controls overfitting; λ is referred to as the penalty parameter; and the superscript T refers to transposition.
By mapping the linear inputs X(i) to nonlinear space ϕ(X(i)), the filter w can be expressed as a linear combination of ϕ(X(i)). In nonlinear space, as opposed to linear space, the feature map can be simply split into the object and background. Additionally, the kernel trick can be used to express w T X(i) as where K(·) refers to the kernel function (e.g., radial basis function) used to compute the dot-product of nonlinear space; K XX represents the kernel matrix with elements K XX (i, j) = K(ϕ(X(j)), ϕ(X(i))); the ith rows of K XX are referred as ; and α is the vector of coefficients α(i) (i = 1, · · · , D). Similarly, w 2 2 can be expressed as Thus, (1) can be reformulated as The solution of (4) is expressed as where y is defined as y = [y(1), · · · , y(D)] T . I refers to a D × D identity matrix. The filter w is replaced by α. K XX + λI can be reformulated as circulant matrix C(K XX + λδ)(δ = [1, 0, · · · , 0] ). K XX refers to the first row of K XX . By using the discrete Fourier transform (DFT), (5) can be expressed asα where the hat denotes the DFT of a vector. When KCF filter is used, the filterα T −1 is updated by (6) refers to a target-centered patch of last frame T − 1 whose target is located in the center. Suppose that a new patch x T is sampled in the current frame T , and the location of the target is detected by where F −1 denotes the IDFT operator.Kz T x T is computed by x T and sample modelz T . The position of the target is determined by the position of the maximum value of R T . Typically, the classical KCF does not use the filterα derived from (6) directly when an object is detected. Instead, linear interpolation is performed to update the filterα T as follows: where β represents the coefficient of interpolation. Similarly, the sample modelz T is also derived from linear interpolation where γ is also referred to as the coefficient of interpolation. By using (8) and (9), the tracker can be updated at a fixed learning rate.
Updating the filter with a fixed proportion can make the filter retain part of the historical information and avoid the interference from new information. However, the fixed learning rate also causes the tracker to treat all the historical information indiscriminately, which causes the waste of historical information. Meanwhile, since the initial value of the tracker is obtained from the first frame, the information contained in the first frame will be exponentially attenuated when the tracker is updated [see (8) and (9)]. In case of a low learning rate,α 1 and z 1 may account for a larger proportion in the filter and sample model, respectively, in the early stage of tracking, which has a significant impact on the tracker. Therefore, the tracker is capable to recognize the object, which is similar to the first frame. Hence, when the feature of the object or the background changes, the tracking may fail.
To make full use of historical information, our tracker retains the filter and the target-centered patch for each frame as HMs , where s i represents a score of each HM. As for the approach to score learning, it will be detailed in Section III-B. Suppose that a new patch x T is sampled in the current frame T , then the filterα T is calculated bŷ The sample modelz T is calculated bȳ where linear interpolation is still performed to update our filter and sample model at the learning rate β and γ, respectively. At each frame,α T −1 andż T −1 are calculated. They both consist of the HMs. The score s i (i = 1, · · · , T − 1) leads toα i (i = 1, · · · , T − 1), which is more sensitive to x T accounting for a larger proportion inα T −1 . In the same way, s i (i = 1, · · · , T − 1) leads to z i (i = 1, · · · , T − 1), which is more similar as x T accounts for a larger proportion iṅ z T −1 . In these operations,α T −1 andż T −1 classify the object and background in x T with high precision. In the meantime, it makes full use of the HMs, with the historical information treated differently. In addition, when the feature of the object or background changes, the scores can be used to reduce the impact of the initial value by adjusting the proportion of the initial value dynamically, because the scores will change in each frame. After the filterα T and the sample modelz T are obtained, the response of x T can be calculated by The position of the maximum value of R T is that of the object.

B. HM Evaluation Scheme
In Section III-A, it is proposed to update the tracker by using HMs. In order to improve the capability of the tracker to distinguish between the object and the background, it is necessary to improve the sensitivity of the selected filters to the object. Meanwhile, the selected sample models must be similar to the object. In this section, an HM evaluation scheme will be introduced to obtain the appropriate filters and sample models by scoring the HM.
The cross-correlation function is applied to measure the correlation of two signals in signal analysis. The discrete cross-correlation function is expressed as follows: where φ 1 and φ 2 represent signals. S(τ ) is referred to as the similarity of two signals when one of them is applied to a τ -step discrete circular shift. The higher the value of S(τ ), the higher the similarity of the two signals. The maximum value of S(τ ) (τ ∈ (−∞, +∞)) can be taken to represent global similarity. Meanwhile, the DFT is used to define the discrete cross-correlation function as where S represents the vector of S(τ ) (τ ∈ (−∞, +∞)) and the superscript H refers to conjugate transposition.
By taking the regression targets and the response as signals, the global similarity of them can be calculated. In theory, the better the performance of the tracker, the higher the global similarity between the response and the regression targets. According to this theory, an HM evaluation scheme is designed to score the HMs.
In frame T , after position detection is completed, a new targetcentered patch z T can be obtained and a new filterα T can be calculated by (6). HMs are updated to . In order to obtain s i (i = 1, · · · , T ), each HM is combined as a tracker in the first place. Subsequently, the responses can be received on z T Then, the discrete cross-correlation function is applied on R i (i = 1, · · · , T ) and regression targets y The score s i (i = 1, · · · , T ) for the ith HM in frame T is given by where max(·) denotes the maximum of the vector; Rank(max(S i )) indicates the index of max(S i ) in the set {max(S i )} T i = 1 , which is ranked in descending order; and θ refers to the threshold.
By using this scheme, the HM, which performs better in distinguishing between the object and background in z T , is given a higher score. Meanwhile, since the HM with low global similarity is irrelevant to confirming the object location, the scores of this kind of HMs are set to 0. All scores are used to update the tracker [see (10) and (11)].
Under this scheme, the historical information can be applied flexibly. Besides, the filter and the sample model are made more robust for object tracking.

C. Antidrift Tracker Correction Scheme
It is a challenge to prevent tracker drift for object tracking. In general, tracker drift occurs when object occlusion happens and similar objects appear. In these cases, a significant deviation can arise between the maximum position of the response map and the ground truth of the object. Besides, the response map fluctuates significantly or shows multiple peaks. Therefore, it is necessary to monitor the tracking state and correct trajectory in the process of tracking. In this section, it will be explained how the KF and the PK of the response map can be applied to correct the tracker.
Assume that the ground truth of the vehicle is represented by {p t } T t = 1 . The motion of the vehicle can be expressed as dynamic equation p t = f t (p t−1 , u t ). In the process of tracking, the position detected by our tracker can be expressed as {d t } T t = 1 . The relationship between detected position and ground truth is expressed as measurement equation d t = h t (p t , v t ), where both v t and u t are the noises that are independent and uniformly distributed. If u t and v t are Gaussian while dynamic equation and measurement equation are linear, KF is applicable to estimate the motion of the vehicle [51].
The dynamic equation and measurement equation of KF are expressed as where A denotes the state transition matrix and B denotes the measurement matrix. The aim of KF is two-fold. One is correction and the other is prediction. Correction means correcting the measured position, in our case is the detected position d t , to an estimate of the ground truthp t . Prediction means usingp t to predict the ground truthp t+1 in the next frame. Herein, the hat denotes the estimated value, and the bar¯denotes the predicted value. The way of correction is defined as The way of prediction is defined as where G t denotes the Kalman gain; C t indicates the covariance matrix of error between ground truth and estimation (p t −p t ); C t represents the covariance matrix of error between ground truth and prediction (p t −p t ); V t refers to the covariance matrix of v t ; U t stands for the covariance matrix of u t ; and I denotes an identity matrix.
In the current frame, the detected position d t is inputted into the KF to obtain a predicted positionp t+1 for the next frame. If the tracker is found to have drifted in the next frame, the detected position d t+1 will be replaced by the predicted positionp t+1 . To measure the reliability of tracker, the high-confidence criterion proposed by Han et al. is adopted [50]. It relies on PK of the response map to detect whether the tracker has drifted. The two measurement thresholds are defined as where P tr and K tr represent peak threshold and kurtosis threshold derived from the historical average values with certain ratio δ 1 and δ 1 , respectively; max(R i ) refers to the maximum of response that denotes the peak value; and BK i denotes the kurtosis value. Ideally, the response map resembles a Gaussian distribution, with a single peak and smooth surroundings. When the drift occurs, the response map exhibits multiple peaks and fluctuates significantly, with the PK in decline (see Fig. 2). Therefore, if the PK of the response map at the current frame reaches above the PK threshold, it can be considered that the response map at the current frame is reliable and no drift occurs to the tracker. If one indicator falls below the threshold, the position detected by the tracker is unreliable and it is necessary to replace it with the position predicted by KF. Meanwhile, the correction of KF and the update of tracker will be terminated until the tracker returns to normal. Fig. 2 shows the effect of antidrift tracker correction scheme. When PK decreases, the tracker with the scheme ceases to be updated and KF predicts the position of the object. Since the incorrect information is not accepted by tracker, the tracking returns to normal when the object is no longer obscured. On the contrary, the tracker without the scheme keeps updating and cannot find the position of the object. When the object appears again, the tracker becomes ineffective. This scheme is effective in avoiding tracker drift through the combination of KF and PK. Under the strategy of updating termination, the tracker and KF avoid receiving false information, which ensures the reliability of the tracker and the accuracy of KF. When the tracker stops working, the prediction function of KF can be used to predict the trajectory reasonably and ensure the continuity of the trajectory. The full algorithm is detailed in Algorithm I.

IV. EXPERIMENTAL RESULTS AND ANALYSIS
A. Experimental Setups 1) SatSOT Dataset: As a densely annotated satellite video single-object tracking benchmark dataset, the SatSOT dataset [49] covers four categories of objects, including vehicle, train, plane, and ship. Herein, the videos of vehicle are used, involving 65 objects, 19 948 frames, and 9 challenges. The object is annotated by a bounding box and recorded with upper-left coordinate, the width, and height of the bounding box. Table I lists the description of nine challenging attributes given by the article [49].
2) Evaluation Metric: The evaluation metric used in the experiments is one-pass evaluation (OPE) [1], [2], with two evaluation indicators involved to perform quantitative evaluation, i.e., precision and successful rate. By plotting precision plot and success plot, the performance of algorithms can be visualized.
Precision is an indicator used to evaluate the error between the detected position of the tracker and the ground truth. The location error of the tracker can be obtained by calculating the Euclidean distance between the center point of the tracking box and the center point of the ground truth. Typically, there is a threshold of location error (e.g., 20), and the percentage of frames whose location error is smaller than this threshold is referred to as the tracker's precision. In the meantime, precision plot is created. The horizontal axis of the precision plot is the location error threshold, and the vertical axis is the precision. Success rate is an indicator used to evaluate the overlap between the tracking box and the ground truth. Suppose that the tracking box is R t and the ground truth is R 0 , then the overlap of tracker is expressed as S = |R t ∩ R 0 |/|R t ∪ R 0 | , where |R t ∩ R 0 | refers to the number of pixels at the intersection of R t and R 0 , and |R t ∪ R 0 | refers to the number of pixels at the union of R t and R 0 . The success rate is expressed as the percentage of frames whose overlap exceeds the overlap threshold. Also, the success plot is created. The horizontal axis of the success plot represents the overlap threshold, and the vertical axis indicates the success rate.
Based on this evaluation metric, the location error threshold is set to 5 pixels to sort the algorithms on precision. As for success rate, the average overlap score (AOE) is introduced as the basis for sorting, which is an average of overlap. Hence, all algorithms are evaluated from two perspectives.
3) Implementation Details: Our experiments are performed on a PC with 2.40-GHz CPU, and our tracker is implemented in MATLAB 2018a. The sample patch is 3.5 times the size of the object's ground truth. Besides, the standard deviation of the Gaussian surface, which is used to fit the regression targets, is 0.125. By comparison, the standard deviation of the kernel function is 0.6. The penalty parameter λ is 1e-4, whereas the learning rates β and γ are both 0.02.
The features extracted in the tracker are grayscale feature and CN [52]. Additionally, the principal component analysis is conducted to simplify the features to ensure the computational speed of the tracker. The features of the vehicle on satellite videos are weak and the size of the object is tiny. Normally, the texture feature is not suitable for representing the vehicle, e.g., HOG and local binary patterns [53]. This is mainly because these features must be extracted from the minimal cells, e.g., 4 pixels. Whereas, the object is too tiny to be divided into cells, so that it may make the loss of information and be meaningless. As stated by Shao et al. [12], the texture feature is ineffective for object tracking on satellite videos. As a result, simple features, such as grayscale and CN, are much more reliable.

B. Sensitivity Analysis and Varying Criteria Comparison
In this section, we first conduct a sensitivity analysis to determine the most effective usage of the HM evaluation scheme. Then, we conduct varying criteria comparative experiments to compare PK with other high-confidence criteria.
In the sensitivity analysis, we vary the start frame ϑ of the HM evaluation scheme and test the sensitivity of ϑ. Besides, we test the sensitivity of θ in (17). When test ϑ, θ and the other parameters in the tracker are fixed. When test θ, ϑ is set with the optimal one received from the last test and the other parameters in the tracker are also fixed. Table II shows the OPE results of two tests. The tracker works best when ϑ = 10 and θ = 70.
In the varying criteria comparative experiments, we compare the PK with other high-confidence criteria, including the average peak-to-correlation energy (APCE) [54], the peak, and the kurtosis. We designed four trackers, and Table III lists the OPE results of them. HMTS_APCE uses APCE as the high-confidence criterion. HMTS_PEAK and HMTS_KURTOSIS employ the peak and the kurtosis of the response map as the criterion individually. They are compared with our tracker HMTS, which employs both the peak and the kurtosis of the response map as the criterion. From the OPE results, it can be found that high-confidence criterion can indeed improve the precision and robustness, and the PK is more suitable for our application.

C. Ablation Experiments
As mentioned earlier, our tracker treats KCF as the baseline, with the HM used to improve the update strategy of KCF. At the same time, KF is introduced for motion estimation, and the PK of the response map is used to detect whether the tracker drifts. To verify whether the HM, KF, and PK are effective in improving tracker performance, ablation experiments were conducted.
We design four trackers. Table IV lists the composition of them. Among these trackers, KCF_CN uses only CF, HMTS_NKFPK does not rely on KF and PK to correct the trajectory, and HMTS_NHM does not require HM for updating. They are compared with our tracker HMTS. The OPE results of these four groups of trackers are shown in columns 5 and 6 of Table IV.
From the performance of HMTS_NHM, it can be found out that the introduction of KF and PK improves the precision and success rate of the algorithm, which is because KF and PK terminate the tracker updating and predict the trajectory of the object when the tracker drifts. From the performance of HMTS_NKFPK, it can be discovered that using HM plays no role in improving the precision and success rate, which is mainly because the baseline tracker tends to drift when the object is occluded and similar objects appear. Since the antidrift scheme is not used, the tracker does not stop updating. At the same time, the HM evaluation scheme scores HMs based on the drifting frame, and these HMs render new filter unreliable. Consequently, the tracker is less likely to return to normal, which reduces the performance of HMTS_NKFPK without KF and PK control.
Our tracker HMTS combines HM, KF, and PK. Also, the HM evaluation scheme keeps scoring HMs based on the correct sample. In this way, the updated filter will become more reliable. Hence, our tracker can take full advantage of HM.
In summary, the introduction of HM, KF, and PK is verified as effective. Our tracker shows a significant improvement over the baseline on the precision and success rate.
1) Overall Results: The overall OPE results are shown in Fig. 3. Our tracker HMTS ranks top on both indicators. The precision of HMTS is 0.7253, and the AOE of HMTS is 0.4344. Compared with the classical CF-based tracker KCF, our tracker HMTS improves the precision and AOE by 0.05 and 0.0331, respectively.
Meanwhile, as shown in the ranking list, CF-based trackers, e.g., KCF, MCCT, ECO-HC, STAPLE, and CFME, are advantageous over end-to-end networks based tracker, e.g., SiamFC and transformer-based tracker, e.g., STARK. That is to say, CF-based trackers are more suitable for vehicle tracking on satellite videos.
2) Attributes-Based Results: Overall, HMTS has excellent performance on most of challenges. Since the antidrift tracker correction scheme are effective on solving the object occlusion  Tables V and  VI). When these challenges arise, the antidrift tracker correction scheme plays a dominant role in preventing drift. When these challenges are addressed, the HM supports the tracker in returning to normal rapidly and precisely.
Meanwhile, tiny object is the object with a small size. There are few features on itself, so it is better to use background information as negative sample. The HMs make HMTS collect rich background information, which is favorable to tracking tiny object. Hence, our tracker has an excellent performance on TO (see Tables V and VI).

E. Qualitative Comparison
In order to better demonstrate the performance of each tracker, qualitative experiments are carried out, with the tracking box of each tracker at each frame displays the video frame. Partial performance is shown in Fig. 4. 1) Qualitative Analysis: From video car_01 and car_24, it can be seen that HMTS is effective in solving the occlusion issue while most of trackers drift when hitting the obstacle. In video car_01, the object is occluded at frame 173, only CFME and HMTS overcome this challenge. The reason is CFME and HMTS use the motion estimation scheme, whereas other trackers do not perceive the disappearance of the object and collect false information when object is occluded. Meanwhile, the time of occlusion is long, and the object moves a long distance. Hence, when the object occurs again, the object is not in the ROI of them, and it is difficult for them to detect the object after collecting so much false information. On the contrary, the antidrift tracker correction scheme makes HMTS stop collecting information and predict the object's position when object is being occluded. Even with long occlusions, the object will be in our ROI. This scheme not only makes HMTS reliable but also leads to a continuous trajectory.
Meanwhile, from video car_08 and car_26, it can be found out that the performance of our tracker is superior in case of SOB, TO, IV, etc. Although the feature of object itself is insufficient to identify it, HM is beneficial for HMTS to collect the background information. As a result, HMTS can ensure accurate tracking.
2) Inferior Cases: Although HMTS is verified as the optimal tracker, there are some inferior cases within the tracking process (see Fig. 5).
HMTS does not use the scale-adaptive mechanism, which leads to the tracking box cannot be adjusted along with the change of the object size, as shown in Fig. 5(a). It results in a lower success rate. Besides, the prediction function of KF is limited. If the velocity and orientation of the object change during occlusion, the KF prediction is likely to be unreliable and may result in drift, as illustrated in Fig. 5(b).

V. CONCLUSION
In this article, a tracker (HMTS) is proposed for vehicle tracking on satellite videos. In addition, an HM evaluation scheme is proposed for the evaluation of HMs by using cross-correlation function. Besides, to prevent tracker drift, an antidrift tracker correction scheme is proposed. Our tracker is compared with 16 state-of-the-art trackers. As demonstrated by the quantitative and qualitative experiments, HMTS produces an excellent performance. Additionally, sensitivity analysis, varying criteria comparative experiments, and ablation experiments reveal that HM, KF, and PK are effective in improving the accuracy and robustness of the algorithm.


Uncertainty-Aware Manipulation Planning Using Gravity and Environment Geometry

Factory automation robot systems often depend on specially-made jigs that precisely position each part, which increases the system's cost and limits flexibility. We propose a method to determine the 3D pose of an object with high precision and confidence, using only parallel robotic grippers and no parts-specific jigs. Our method automatically generates a sequence of actions that ensures that the real-world position of the physical object matches the system's assumed pose to sub-mm precision. Furthermore, we propose the use of “extrinsic” actions, which use gravity, the environment and the gripper geometry to significantly reduce or even eliminate the uncertainty about an object's pose. We show in simulated and real-robot experiments that our method outperforms our previous work, at success rates over 95%.

However, these jigs are expensive, to the point that their cost including engineering and system integration can make up over 50% of the total cost [3]. Their versatility is also limited, since only few jigs can be placed within reach of each robot arm and the changeover to a different product line incurs additional costs. It would be highly desirable for robot systems to handle a large variety of parts without the need to change jigs.
Programming robots without jigs poses difficulties, such as noisy pose estimation from parts detection which can introduce misalignment of the grasped parts and result in task failure. For the World Robot Summit Assembly Challenge 2018 and 2020, we developed a robotic system composed of multiple robot arms and generic grippers aiming to realize jig-less assembly [4]. Although our system won 4th and 3rd place in 2018 and 2021 respectively, our robots (and those of the other winning teams) frequently failed insertion tasks due to position uncertainty.
Previously, we developed an in-hand pose estimation method in which the robot touches the environment with the grasped part to reduce the uncertainty of the part's pose to sub-mm precision [5] (the Touch action). We extended this by also evaluating the effect of a camera view on the object pose uncertainty (the Look action) and allowing a combination of both actions [6]. However, these methods were limited: the user had to define the action sequence manually, and the actions required either slow motions and a force sensor (for the Touch action), or careful camera calibration (for the Look action).
In this work, we solve both of those problems by first extending our method with three new actions (Grasp, Place, Push) which require no force sensor or camera calibration, and which are efficient and fast as they take advantage of both gravity and the environment geometry. Furthermore, we formulate a planner that selects action sequences which reduce or even eliminate object pose uncertainty, thus enabling reliable manipulation planning for tasks that require high precision, such as industrial assembly.
We show that our proposed method performs reliably in simulation for 13 parts used in a model assembly task, as well as in a real robot system. Using the three new actions, the uncertainty of many objects can be eliminated without the use of any sensors. The only constraints on the object are that the robot needs to be able to grasp it stably, that it may not topple or roll during pushing, and that it may not roll away when placed (either by using a placing surface of appropriate softness or by disallowing spheres and cylinders). r The formulation of five actions that reduce the pose uncertainty of an object r The calculation of the uncertainty that results from each action r The planner that minimizes uncertainty using these actions.

A. Uncertainty-Aware Manipulation
Creating "robust" manipulation plans that succeed even under uncertainty (also called "open-loop," "oblivious" or "sensorless" manipulation plans) has a long history. Mason described the "Place" action as early as 1982 [7], and investigated the effect of pushing actions. Erdmann and Mason [8] showed that sensorless manipulation can arrange planar objects into a known state. Goldberg [9] described an algorithm to obtain a sequence of grasps that orients a planar polygonal object without knowledge of the initial orientation. They assume infinitely large grippers and no friction. Zhou et al. [10] formulated a planner for planar grasping sequences that allow arbitrary shapes with bounded uncertainty to be arranged into a known state, by modelling friction and the contact between gripper and object. In contrast to the works above, our method is not restricted to 2D, but works with 3D meshes that are widely used in motion planning and collision checking libraries. Allowed grasps can be defined anywhere on the object, and the uncertainty is minimized for 3D poses.
Lozano-Pérez et al. [11] proposed "pre-image backchaining," an algorithm to obtain a robust manipulation ("fine-motion") plan by defining permissible regions of space backwards from the goal position, from which a motion will certainly achieve its goal. To avoid having to quantify these regions and because the search space is small and can be exhausted easily, our method simply applies a breadth-first search and samples actions from the initial pose and uncertainty.
More recent approaches propose machine learning methods to handle uncertainty. Kahn et al. proposed learning a collision probability using a predictive model for collision avoidance of mobile robots [12]. Kou et al. developed a model-based reinforcement learning approach, where the maximum control input value was scaled down with high model uncertainty [13]. Lee et al. proposed to learn switching a model-based and learningbased controller with low and high uncertainty, respectively [14]. Unlike these approaches, our method does not need to collect data or train machine learning models. This is very desirable for users.

B. Active Perception
Many studies have explored active perception approaches where robots interact with objects to determine the objects' location or properties [15]. In tactile exploration, a probabilistic model has been learned for the target object's surface which efficiently identified the shape by touching to reduce the model uncertainty [16], [17], [18]. Koval et al. and Wirnshofer et al. proposed particle filter-based approaches to estimate object pose by active contacts [19], [20]. We also use particle sampling in some of our actions.
Paolini and Mason learned a statistical model for an object by pick and place regrasping [21]. However, our approach requires no data collection. Chavan-Dafle and Rodriguez demonstrated other actions utilizing the environment [22] as well as stable regrasping without fixtures by sliding the object in the gripper [23]. These approaches are the most similar to our actions, but they do not consider the uncertainty around the object, and they assume that an object can slide within the gripper. Our method explicitly tracks uncertainty, and uses gravity to obtain a stable pose for the object, which is more generally applicable than in-hand sliding.
While most approaches have focused on a single action or several objects, our method combines 5 different actions that apply to a wide range of parts.

C. Object Pose Estimation Using Multi-Modal Sensors
Numerous computer vision and pattern recognition works have developed not only 2D [24] but also 6D pose estimation techniques using RGB-D information [25]. CAD information and deep learning have been used to improve robustness [26], [27]. Although vision-based methods are helpful to determine object poses, they are affected by lighting conditions and materials. Our method is particularly robust in the context of assembly tasks, where metallic parts with reflection may degrade the estimation performance.
Other sensors or sensor fusion for object pose estimation have also been used. For example, tactile sensors have been shown to improve grasped object localization [28], and in-hand object localization approaches have combined vision, force-torque, and joint angles sensors to improve estimation performance [29], [30]. However, while these multi-modal methods improve pose estimation, they do not quantify the measurement uncertainty of the observations. Our method reduces the remaining uncertainty by tracking and minimizing it explicitly.

III. PROPOSED METHOD
In our proposed method, the robot reduces the object pose uncertainty by executing a sequence of actions obtained by a planner (Subsection E). We use five type of actions: Touch, Look, Place, Grasp and Push, of which the last three are newly added in this letter. The pose belief is represented by a Gauss distribution (Subsection A). The planner's objective is to minimize its size (Subsection D). To produce effective sequences, the planner must calculate the effect of each action on the pose belief. To calculate the last three types of actions, we sample particles from the distribution to obtain discrete poses without uncertainties (Subsection B), for which the effect can be calculated geometrically (Subsection C).
The first two actions (Touch and Look) are described in [5], [6]. The Look action requires either a calibrated camera or a calibration geometry in the image, while the Touch action requires the position of at least one calibrated support surface (and optionally one edge) in the environment. The remaining three are extrinsic manipulations using the interaction between the object and the gripper geometry and/or supporting surface. They require only the support surface and no sensor. The Place action consists of placing the grasped object on a support surface with known height. Grasp actions consist of grasping an object resting on a support surface, as shown in Fig. 1. Push actions consist of pushing the object using the gripper geometry.
The method assumes that 1) the geometry is provided as a mesh or collection of bodies, 2) friction can be neglected until the object is grasped, and 3) the object does not topple or roll during pushing.

A. Representation of Pose Beliefs
To model the belief about the object pose, we follow the method described in [31]. We represent a pose T ∈ SO(3) with uncertainty as a mean poseT ∈ SO(3) with a small perturbation by a random variable ξ ∈ R 6 : where ∧ : Let ∨ : se(3) → R 6 be the inverse of ∧ . We assume that the variable ξ follows the zero-mean Gaussian distribution N(0, Σ) where Σ is the covariance matrix. Thus, we represent a pose belief using the mean poseT and the covariance matrix Σ of the perturbation variable.

B. Calculating Pose Uncertainty
Each action reduces the pose uncertainty in different ways. The calculation of the Touch and Look action is described in our previous work [5], [6]. For the three extrinsic manipulations, we assume that the pose distribution after the action can be approximated by the representation in §III-A. For this assumption, it is required that the function from poses before the action to poses after the action is continuous in the region of uncertainty, but not necessarily C 1 -class. Note that it may be discontinuous when an extrinsic action is performed badly (for example, placing a pole on its end in an unstable position), but our planner ( §III-E) avoids such actions.
We calculate the effect of the extrinsic manipulations on pose beliefs as follows. We first generate the pose particles T 1 , . . . , T N by (1) where N is the number of particles. Then, we calculate the poses T 1 , . . . , T N after the manipulation for each particle using the methods described in the next subsection. Finally, we estimate the new pose belief (T , Σ ) as a distribution fitting the samples T 1 , . . . , T N .
To calculateT , we solve the equation using the Newton method. More precisely, letT s be the current approximation of the solution of (2) and let and J s is the Jacobian of the map which can be calculated using (33) and (34) in [31]. After calculation ofT , Σ is calculated by where

C. Calculating Poses After Extrinsic Manipulations
Apart from the uncertainty, we must also calculate the mean pose resulting from the manipulation. For all extrinsic actions, we assume that the grasped object is a polyhedron and its center of gravity is known. We ignore friction between both the object and the gripper, and the object and the support surface.
In the following sections, the object's phases of rotation correspond to Fig. 2. 1) Place Action: Let C be the projected point of the center of gravity of the object onto the support surface. Since we ignore friction between the object and the support surface, the object cannot topple and the position of C in the support polygon does not change. Thus, it is sufficient to calculate the rotation of the object as it is being placed.
The rotation of the object while it is being placed on the support surface is divided into the following four phases: 1) When no vertex of the object is on the support surface, it descends vertically without rotating. 2) When only one vertex of the object is on the support surface, the object rotates around the line on the support surface which is orthogonal to the line connecting the vertex and C. 3) When two vertices of the object are on the support surface and there exists a perpendicular line from C to the segment connecting the two vertices, the object is rotated around the line connecting the two vertices. 4) When more than three vertices are on the support surface, and if C is inside of the convex hull of the vertices of the object on the surface, the object is placed stably. Otherwise, this place action is regarded as unstable. Note that the condition for the third case is always satisfied when the object is placed on the face of an acute-angle triangle or a rectangle. As the rotation can be approximated by the described method even when the condition does not apply, we use the calculation for all cases.
2) Grasp Action: We assume a parallel gripper with two flat gripper pads ("fingers"). We consider the fingers as parts of two parallel planes which are orthogonal to the support surface and move with the direction of their normal lines. We define three orthogonal axes x g , y g , z g as displayed in Fig. 2, where x g is parallel to both the support surface and the fingers, y g is orthogonal to the fingers, and z g points upwards, away from the support surface. When the initial pose of the object is near the stable pose and the center of gravity is not very far from the grippers, the object's change of position with respect to x g is negligible. The position in y g is fixed by the fingers and the position in z g is determined by the support surface. Thus, we only need to determine the rotation of the object during grasping.
The rotation is divided into the following phases, as the gripper closes: 1) When no finger touches the object, it does not rotate.
2) When one vertex touches one of the fingers or two vertices touch one finger each, it rotates around z g . 3) When two vertices touch one of the fingers but no vertex touches the other pad, it is only pushed and does not rotate. 4) When two vertices touch one of the fingers and one vertex touches the other finger, it rotates around the line connecting the two vertices on the same side. 5) When three or more vertices touch one finger and at least one vertex touches the other finger or two vertices touch each finger, the object is grasped stably if there is an intersection of the convex hull of the vertices touching one finger and that of the other finger. Otherwise the grasp is considered unstable. Note that the third phase is not always passed. Note also that, in the conditions of the second, third, and fourth phases, vertices with the same x-and y-coordinate are considered as one vertex.
3) Push Action: We use the same coordinates as for the grasp action. When the push action is executed, the fingers are closed and the object is pushed within the direction of the x g -axis. The reasoning of the grasping case applies analogously, so we only need to calculate the rotation.
The rotation is divided into the following phases: 1) When the fingers do not touch the object, it does not rotate.
2) When at most one vertex touches the fingers, the object rotates around z g . 3) When two vertices touch the fingers, the object is only pushed and does not rotate. Note that, as in the grasp action, vertices with the same xcoordinate and y-coordinate are considered as one vertex.

D. Quantification of Uncertainty
To allow the planner to minimize the uncertainty, we need to quantify it as a scalar. We use the alternative form of (1): Let Σ be the covariance matrix of ξ . The transformation from Σ to Σ can be calculated by the coordinate change formula (26) in [31]. We define the amount of uncertainty by where c i,j are the fixed coefficients. The reason we transform Σ to Σ is that Σ is associated with the object's own coordinate system. The first three diagonal components of Σ correspond to the object's position uncertainty in x, y and z, the second three to the orientation uncertainty. The non-diagonal components correspond to the covariances between them.
Some objects have rotational symmetry around an axis, around which uncertainty cannot be eliminated. In these cases, we set the coefficients c i,j such that unavoidable uncertainty is ignored and the plan succeeds in minimizing the rest. If we did not use Σ instead of Σ, this would be impossible, as the symmetry axis depends onT .

E. Planning
To obtain the action sequence from the planner, we first enumerate the possible actions.
For the Place action, we enumerate all faces of the three dimensional convex hull of the object as possible placement orientations. For the Grasp action, we define the grasp poses for each object by hand for simplicity and calculation speed, but any grasp-generating algorithm could be applied. For the Push action, we enumerate the edges of the projected object on the support surface as pushing directions. For round objects, we randomly choose a lower number of directions when there are too many.
We limit the number of possible actions to around 30, and the lengths of action sequences to 3 because this is sufficient for the majority of objects in our experiments. As the size of search space is O(N d ) where N is the upper bound of the number of possible actions and d is the depth of search, it is small enough that we can use breadth first search. At each step, we select the next action and calculate the resulting uncertainty. For the Touch and Look action, we assume that the object is observed at the expected mean pose, and that the uncertainty decreases according to the selected action.

IV. PLANNING EXPERIMENTS
To verify that A) our planner gives valid results and B) the performance of our newly proposed actions (Place, Grasp and Push) compared to the previous ones (Touch and Look), we generated action sequences for a variety of objects and evaluated the resulting uncertainty.

A. Setup
We evaluated the resulting plans for the following conditions: r We created plans for 13 different objects that were used in the WRS 2020. See Fig. 3, which shows the parts with the initial uncertainties. These values represent a standard deviation of Σ i,j = 3.16 mm in translation and 0.0316 radians in rotation. This means that 99% (three standard deviations) of sampled particles will be within 9.5 mm and 0.95 radians (5.4 degrees) of the reference pose. Then, the object is grasped and the uncertainty is calculated. This grasped pose with uncertainty is the initial pose.
For objects which have no rotation symmetry (base, panel motor, panel bearing, motor), the amount of uncertainty is defined according to (5) with the coefficients: For objects which are symmetric around the x-axis (motor pulley, bearing, shaft, end cap, bearing spacer, output pulley, idler spacer, idler pulley and idler pin), we ignore the rotation uncertainty around the x-axis by defining the uncertainty as in (5) with the coefficients: Note that although the output pulley (Fig. 3(e) is not strictly symmetric, it is close enough that we choose to disregard its rotation uncertainty.

B. Results
The result of the planning experiments are displayed in Table I. The box-plots showing the amount of remaining uncertainty (after 3 actions) are displayed in Fig. 4 When planning with all actions, the results are best, and uncertainty is eliminated or severely reduced, generally to sub-mm levels. The results show that the extrinsic manipulations are responsible for reducing the uncertainty to zero, as this occurs for 8 out of 13 parts (panel motor, panel bearing, motor pulley, bearing, end cap, bearing spacer, idler spacer, idler pulley and idler pin) even when planning only with extrinsic manipulations. For two parts (motor and shaft) the planner could not find a valid plan using only the extrinsic manipulations, as the parts have no faces to place them stably. For the base, the reason why uncertainty remains is that while a face to place it exists, it is so thin that the placing can fail. For the output pulley, the push action can sometimes fail because due to its protruding part, the center of mass and the center of the circle in the outline do not coincide, so some uncertainty remains.
When planning only with static manipulations, uncertainty is reduced for all parts, but never reduced to zero.
In conclusion, the extrinsic manipulations are powerful but not usable for all objects, while the static manipulations are effective regardless of the grasped object.

V. SIMULATION
To verify that the action sequences are valid and can be executed, we first conducted simulation experiments.
We used Gazebo as the simulation environment and simulated the same robot that was used in the real-robot experiment. We choose the bearing plate (see Fig. 3(c) as the target object. Friction parameters mu and mu2 were set to 0.5 for the bearing plate and 1.0 for the gripper fingers and support surface.
Both in simulation and the real robot experiments, we use gripper fingers with a width of 17 mm and a thickness of 6 mm.

A. Setup
For this experiment, we obtained an action sequence from our planner which satisfies the following conditions: r At the initial state, the panel is placed with uncertainty as shown in Fig. 5(a). The uncertainty is represented by a covariance matrix Σ such that r At the final state, the panel must be placed in the same orientation and have no uncertainty. The sequence consists of the following actions: 1) Grasp the plate as shown in Fig. 5(b), 2) Place the plate at the same position, and 3) Push the plate horizontally as shown in Fig. 5(c).
As a control, we also simulate a manually created action sequence that leaves a dimension uncertain: 1) Same initial state as above, 2) Grasp the plate, and 3) Place the plate at the same position as the above sequence.
After the above action sequences, we evaluate the resulting object positions.  of the panel is considered acceptable. For the planned action sequence, placement was successful 97% of the time, while the control sequence only succeeded in 15% of cases.

B. Results
As expected, both the planned sequence and the control sequence decrease the scattering in x to a few millimeters, but only the planned sequence removes it in y.
Some noise remains even after executing the planned sequence, especially for the right hole. This is likely due to numerical instabilities in the gripper model and the contact between the gripper and the plate. Without this effect, we would expect a 100% success rate in simulation.

VI. REAL ROBOT EXPERIMENT
To confirm the effectiveness of our method and the new actions (Push, Place, Grasp) that have been added to the Touch and Look actions, we conducted a real-robot experiment based on the actual assembly task at the World Robot Summit Assembly Challenge 2020. The experiment recreates the placement of the two L-shaped panels, the bearing, and the shaft (shown in Fig. 3(b), 3(c), 3(f), and 3(g), respectively). We evaluated the precise placement of the through-holes of the L-shaped plates with the screw holes of the base, and the alignment of the bearing and shaft.

A. Experimental Setup
We recorded a total of 480 trials. At the start of each trial, the parts were placed at an initial position, defined by a set position P , translated by a random offset within ±10 mm in each horizontal axis and rotated by a random angle within ±15 degrees around the vertical axis. The random offset and rotation represent a significant amount of noise in an object detection pipeline.
As seen in Fig. 5, for the motor panel and bearing panel, the robot first performs a grasp action centering and then a push action. For the bearing and shaft, the robot performs several grasp actions. As this sequence should leave the object pose well defined without any remaining uncertainty, the robot then picks up the part and places it at the target location. The hole alignments are then recorded with an RGB camera, either from above or from the side, and the distance from the reference position estimated manually from the image. The robot then places the parts at a new random location to prepare for the next trial.

B. Results
Of 120 trials each, the motor panel, bearing panel, bearing, and shaft were placed successfully within 1 mm of the target location 106, 114, 115, and 120 times, respectively. The bearing panel and part were placed three times, the motor panel 15 times within 3 mm of the target location but at more than 1 mm distance. We recorded these offsets as a failure due to the high precision requirements for these applications, although this error can be permissible. Three and two trials for the bearing panel and bearing were not counted due to unrelated issues (e.g., motion planning issues or latency). In total, the parts were successfully placed at the target location in 95.7% of the trials (455 of 475).
All of the failures were due to mechanical problems that can be mitigated, such as the gripper surfaces sticking to the part and pulling it away from the assumed position as the gripper opened.

VII. DISCUSSION
The simulated and experimental results show that our method produces feasible action sequences that reduce or eliminate the uncertainty of an object pose, and that they lead to repeatable results on a real robot preparing a representative model assembly task.
In almost all cases, the Extrinsic actions we proposed (Grasp, Place, Push) reduce uncertainty more than the Touch and Look actions. This makes intuitive sense, as aligning an object with a surface (table or gripper pads) constrains multiple degrees of freedom, while the Touch action constrains only a single point. They are also easier to calibrate than the Look action, as they require only the position of a flat surface instead of complete camera calibration.
Due to size and measurement hardware constraints, we could only perform real robot experiments for a limited amount of objects and action combinations. However, previous experiments [5], [6] have shown that even using only Touch and Look actions, uncertainty can be reduced to sub-mm levels, thus confirming the simulation results.
Our method has some limitations. We make a number of simplifications, most importantly that friction has no significant effect. For the gripper, this can be reproduced by adding a linear guide to one of the fingers, or simply starting close to a stable grasp's position. For the Push action, the motion needs to be slow enough to be quasi-static, and the part may not topple. The Push action also assumes that the pushing part of the gripper is a flat surface, which is permissible for many but not all industrial grippers.
We have omitted grasp generation and stability evaluation for the Push and Grasp actions, as for industrial use it is often preferable to explicitly define permissible interactions (predictable behavior, easier tuning than e.g. friction parameters). However, automatic grasp generation and push stability evaluation can easily be added on top of our method, by using them to define the permissible grasps and pushes.
A limitation of the Place action is that it assumes a stable placement pose, which does not apply for meta-stable objects which can roll, such as cylinders or balls. In practice, a layer of foam can inhibit the rolling enough for the action to work.

VIII. CONCLUSION
We presented a method to generate a sequence of actions that reduces or eliminates the uncertainty about an object pose. In real-robot experiments, the sequence resulted in sub-mm precision over 95% of the time, and in <3 mm precision 99% of the time. In simulation, we showed that the "extrinsic" actions proposed in this letter (which use gravity, the environment and the gripper geometry) improve the uncertainty reduction significantly over our previous work, and can eliminate uncertainty entirely. Our method can be used to precisely position objects without custom jigs, or it can be implemented into existing Task and Motion Planning frameworks to increase repeatability and resilience to noise.


Tensor Program Optimization with Probabilistic Programs

Automatic optimization for tensor programs becomes increasingly important as we deploy deep learning in various environments, and efficient optimization relies on a rich search space and effective search. Most existing efforts adopt a search space which lacks the ability to efficiently enable domain experts to grow the search space. This paper introduces MetaSchedule, a domain-specific probabilistic programming language abstraction to construct a rich search space of tensor programs. Our abstraction allows domain experts to analyze the program, and easily propose stochastic choices in a modular way to compose program transformation accordingly. We also build an end-to-end learning-driven framework to find an optimized program for a given search space. Experimental results show that MetaSchedule can cover the search space used in the state-of-the-art tensor program optimization frameworks in a modular way. Additionally, it empowers domain experts to conveniently grow the search space and modularly enhance the system, which brings 48% speedup on end-to-end deep learning workloads.

Introduction
Deep learning has become pervasive in daily life. From video understanding [28], natural language understanding [15], and recommendation system [29] to autonomous driving [21], different deep learning models are deployed on different hardware platforms and devices. Deep learning frameworks usually rely on manually optimized libraries [13,22] to accelerate deployment. Engineers need to choose from many tensor programs that are logically equivalent but differ significantly in performance due to memory access, threading, and the use of specialized hardware primitives. The engineering effort required for tensor program optimization has become a significant bottleneck for machine learning deployment with the growing number of models and hardware backends. Figure 1: Automatic tensor program optimization contains two key elements: the search space S(e 0 ) and the search algorithm that finds the optimal tensor program e . The search space usually incorporates choices over loop transformation, vectorization, threading patterns, and hardware acceleration.
Automatic program optimization [10,43,1] is a recent sequence of efforts that aims to use machine learning to solve this problem. There are two vital elements of automatic tensor program optimizations. First, a search space is defined to provide a possible set of equivalent tensor programs. Then existing systems use learning-based search algorithms to find an optimized tensor program in the search space with feedback from the deployment environment. Most of the current approaches [10,43,1,25,4] use pre-defined search spaces that effectively encode the domain knowledge of the authors once and focus on developing efficient search algorithms.
While efficient search is essential, the search space itself fundamentally limits the best possible performance search algorithms can get. To construct a good search space, domain experts have to make numerous choices over loop transformation, vectorization, threading patterns, and hardware acceleration. Additionally, the best search space itself evolves as new tensor program optimization techniques [24] and hardware primitives [30] grow. As a result, there is a strong need to enable easy customization and construction of the search space at scale by taking inputs from system engineers and domain experts. Unfortunately, any change to search space construction currently requires surgical modifications to the automatic program optimization frameworks.
This research asks the following question: can we decouple the search space construction from search and provide an adequate abstraction for domain experts and the learning system to collaborate on search space construction? We give an affirmative answer to the question with two key observations. First, we can parameterize an optimization search space by the initial program followed by a sequence of transformations on the program. Next, using this parameterization, domain experts can then provide probabilistic choices that represent possible transformations after examining the program state. These two observations lead to a simple yet powerful abstraction for search space construction through a domain-specific probabilistic language. Finally, our framework composes multiple possible probabilistic transformations to form a rich search space. We make the following contributions: • We introduce a simple yet powerful probabilistic language abstraction for tensor program search space construction.
• We build a learning-driven framework to find optimized tensor programs specified by the search space constructed using our abstraction.
• We build an end-to-end system that can take prior knowledge from domain experts to construct optimization search space to optimize deep learning deployment on multiple platforms.
Our end-to-end system can easily expand search space that matches previous approaches without surgical changes and achieve comparable performance on multiple hardware backends. Experimental results show that our abstraction is expressive enough to cover the optimization space of a diverse set of tensor programs, delivering a competitive performance of popular deep learning models, and convenient to incorporate hardware-specific knowledge into the search space to outperform the state-of-the-art frameworks.  2 Background and Problem Overview Figure 1 shows a typical workflow for tensor program optimization. For a given program e 0 , a typical tensor program optimization framework will generate candidates from a pre-defined search space S(e 0 ) containing semantically-equivalent programs. Then the framework finds optimized tensor program e * ∈ S(e 0 ) with the minimum latency on the target hardware.
A typical search space S(e 0 ) contains choices over threading, loop ordering, memory access, and hardware primitives. Defining the search space S(e 0 ) for a wide range of tensor programs brings several challenges. First, S(e 0 ) is highly dependent on e 0 . For example, S(e 0 ) of a compute-intensive program (e.g., Dense) needs to consider many more possible configurations than a communicationintensive program such as ReLU. The space also differs significantly in different hardware domains. For example, S(e 0 ) on CPU involves multi-core parallelism and vectorization, while S(e 0 ) on GPU involves thread binding and tensorization. Finally, as the hardware and model settings change, we need to bring in fresh domain experts to update the S(e 0 ) to leverage the latest improvements.
This paper aims to provide a programmable abstraction to construct S(·) in a composable and modular way. Our key goals are listed as follows: Expressiveness. We need to be able to build a rich search space that covers the optimization programs that domain experts will write. Modularity. Tensor program optimization likely will involve inputs from multiple domain experts over different periods. Therefore, we need to be able to combine prior knowledge in a composable and modular way. Designed for learning. We need to build a generic learning-based framework to enable diverse variations of the cost model and search for search space specified in our abstraction. We will address the above goals in the following two sections.

Composable Search Space Parameterization
This section presents MetaSchedule, a probabilistic approach to search space parameterization.

Stochastic Search Space Construction
MetaSchedule constructs a search space S(·) with stochastic program transformations as the primitive building blocks. Traditional program optimization can usually be represented by a sequence of transformations τ , where at step i, the program e i−1 is transformed into a semantically-equivalent program e i , which finally leads to the optimized program e n . MetaSchedule generalizes this idea by allowing further parameterization of each transformation step in τ .
Taking Figure 2 as an example: e 0 is the initial program for the program B = ReLU(A) 1 . In MetaSchedule, transformation t 1 = Split is parameterized by a loop i and a sequence of integers indicating the loop extents after splitting; Similarly, transforms t 2 = Parallelize and t 3 = Vectorize are parameterized by loops respectively. As a result, an optimized program e n is obtained by applying a sequence of parameterized transformations τ to the initial program e 0 . Accordingly, the search space S(e 0 ) is composed of e 0 and all possible sequences of parameterized transformations. for i1, j1 in grid(θ1, θ3): for k in range (16) On the other hand, it could be less practical for practitioners to determine the best combination of the parameter values in transformations. For instance, in Figure 2, it is usually efficient to use 512-bit vectorization over the inner loop when AVX-512 2 vector instructions are available on Intel CPUs, or other vectorization lengths may lead to better performance, otherwise. Therefore, deep knowledge of the target hardware is mandatory to enumerate plausible parameter combinations to control the search space size while covering the optimal program.
To let practitioners efficiently define parameterized transformations without worrying about candidate values, MetaSchedule introduces random variables drawn from analysis, sampling. Parameterized by random variables, a transformation naturally becomes stochastic, and the underlying probabilistic space reflects the space of possible transformations.
As illustrated in Figure 3, when creating Split transformations to tile loop i and j in the Dense operator, the tile sizes are drawn by random variables θ 0−3 defined from Sample-Tile. In this way, the Split transformation becomes stochastic. Similarly, we use Sample-Compute-Location to enumerate valid loops in Dense after splitting for ReLU to fuse its computation. In summary, 7 lines of MetaSchedule program covers a family of possible optimized tensor programs with stochastic transformations in its search space S(e 0 ), where e 0 is Dense-ReLU.
Notably, unlike orthogonal grid space in hyperparameter tuning, MetaSchedule captures longterm structural and arithmetic dependency between random variables and the tensor program e i being transformed. As demonstrated on Step 2 in Figure 3, sampling distribution from Sample-Compute-Location depends on the latest tensor program e 5 , whose structure depends on all previous random variables.

Modular Search Space Composition
Although the search space constructed by stochastic transformations proposed in the previous subsection is efficient and is capable of covering the optimal tensor program, it is hard for other developers to learn how the search space is constructed by reading a long sequence of transformations. It makes transformations designed for a workload hard to be reused by other workloads. Meanwhile, we observe that it is straightforward to group a sub-sequence of transformations for a particular fine-grained optimization. For example, some transformations implement multi-level tiling for better memory locality in compute-intensive operators like Conv2d and Dense; some other transformations are used to fold/inline elementwise operations such as activation functions into their predecessors or successors for better memory bandwidth efficiency.
To improve the usability and make MetaSchedule more practical, we introduce transformation module. Just like the convolutional module with Conv2D, BiasAdd and ReLU in ResNet, a transformation module in MetaSchedule is defined as either atomic stochastic transformation, or composition of for i in range (512): for j in range (256): for k in range (16)  program analysis, sampling as well as smaller transformations. Each transformation module can have a meaningful name so that it can be easily adopted by many workloads to hierarchically construct a search space. Figure 4 shows hierarchical composition of transformation modules. Specifically, Multi-Level-Tiling interleaves program analysis on loops and the stochastic tiling of the loop structure and organizes the original tensor program into a 5-level tiling structure. Notably, the transformation module is generic to tensor programs and thus could be applied to a variety of operators, including conv1d, conv3d, matmul, etc. Figure 5 depicts an example of composing a search space with transformation modules. In this simple example, we select a set of transformation modules, which are implemented in advance by practitioners with prior domain knowledge, and apply them to every available location in the tensor program to form a search space. Consequently, the formed search space covers common optimizations on diverse hardware.

Relation to Existing Tensor Program Optimization Methods
In this subsection, we discuss prior approaches for automatic tensor program optimization and illustrate that many of them can be covered by the MetaSchedule framework.
Domain specific languages for program transformations used by prior frameworks [32,5,9,37] allow developers to easily optimize a program manually. When there is no random variable sampling in the program, MetaSchedule reduces to a DSL for deterministic program transformations and achieves the same functionality.
Template-guided auto-tuning [10,2,25,27] fully relies on developers to define a search space. In MetaSchedule, it means all random variables in a search space are defined ahead of the transformations, so there is no interaction between program analysis and follow-up random sampling choices conditioned on the program state.  Auto-scheduling [43,45,1,19] requires developers to implement workload agnostic transformation rules. MetaSchedule achieves the same programmability and functionality through specific probabilistic transformation modules that correspond to the search space generation rules.
Notably, all approaches mentioned above have important use-cases in tensor program optimizations, depending on how much domain knowledge we want to incorporate for a particular scenario. By decoupling the search space construction from the search, we effectively build a single framework for all the use cases and enable further customization without surgical changes to the system.

Learning-driven Search
The last section provides a modular abstraction for search space. We still need to do an effective search to find an optimized program within the search space. This section provides a generic learning-driven framework to find an optimized program.
Objective formalization. For a given probabilistic program e 0 , let us use τ to denote the transformations performed on e 0 . τ can be sampled from a prior distribution specified by the probabilistic program. We define g(e 0 , τ ) to be the tensor program after applying transformation τ to e 0 . Let f (e) be the latency of the particular program e on the hardware environment. We define a posterior probability of an optimized program as: P (τ | e 0 ) ∝ e −f (g(e0,τ )) · P (τ ).
(1) Intuitively, we want to assign a higher probability to the programs that perform well. Our final goal is to find τ = argmax τ P (τ | e 0 ) that maximizes the posterior through maximum a posteriori estimation (MAP) estimation.
Execution tracing. To enable domain experts to express their knowledge via transformations modules productively, we embed MetaSchedule in Python. We introduced execution tracing to reduce the cost of repetitive re-execution of the Python program. Figure 6 demonstrates an example tracing process. During program execution, the system records all samplings and transformations while ignoring control flow and other constructs of the host language. The resulting trace is a sequence of MetaSchedule primitives with only sampling and transformation instructions, which could be re-executed as a normal MetaSchedule program. We can then continue to explore different sampling choices for a given collection of initial traces. Conceptually, this is equivalent to dividing up our support set and then sampling the program condition on the execution sequence of the program.
End-to-end search. Figure 7 shows the overall workflow of our learning-driven framework. The search algorithm first samples the MetaSchedule program to obtain a collection of traces. Then it continues to explore the space condition on the traces. Notably, there is a significantly higher cost measuring f (e) directly on the hardware, so we also incorporated a proxy cost modelf (e), which is updated throughout the process, similar to previous efforts on tensor program optimization [10,43]. At each iteration, we adopt an evolutionary search algorithm that proposes a new variant of the trace by mutating the random variables, then accept or reject the proposal based on the cost model. While evolutionary search could be viewed as parallel chain MCMC, we also made our system modular enough to incorporate other ways to select the probabilistic choices, such as those through Bayesian optimization and reinforcement learning. Cost model. Our approach allows extensive cost models, enabling us to supply those pre-trained from existing datasets [44]. We pick a tree-boosting-based cost model inf (·) by default and leverage a common set of features that are used in previous works [43].

Multi-Level Tiling
Trace validation. Importantly, invalid traces may show up as we propose updates. Such a scenario can happen when some of the random variable choices go beyond the physical hardware limit or a variable that induces changes to the execution sequence. Instead of enforcing a conservative proposal, we introduce a validator that validates the correctness of the trace. The trace validation allows us to move around the space more freely while still ensuring the correctness of the sample outcome to be on the right support set.

Related Work
Tensor Program Transformations are proposed by many prior works, such as Halide [32], TVM [9], Tiramisu [5] and TACO [23,37]. Note that all the previous transformation languages are deterministic and cannot be directly used for search space construction, meaning that they have to introduce a separate programming model to express a search space. This paper makes a simple but powerful generalization to domain-specific probabilistic language. The resulting abstraction enables a unified approach to deterministic transformation and search space construction.
Black-box optimization has been adopted in high-performance computing libraries [16,3]. Recent advances in automatic tensor program optimization brought a rich set of techniques to accelerate search through better cost modeling [10,4,34] and learning-based search [2,25,1,19,44], which could be incorporated into MetaSchedule search. Different variations of pre-defined search spaces have also been proposed that couple with the automatic tensor program optimization frameworks [10,43,1]. Polyhedral model [40,6,39] is one useful way to construct a rich pre-defined search space. This paper focuses on modular search space construction and provides orthogonal contributions to these prior works.
Probabilistic programming language is a powerful abstraction for incorporating domain knowledge and probabilistic inference. There are many general-purpose probabilistic languages, such as Church [17], Stan [8], Pyro [7], NumPyro [31], PyMC3 [35] and Edward [38]. This paper proposes a domain-specific probabilistic language for tensor program optimization with specializations such as tensor program analysis that would otherwise be opaque to the previous systems. Our learningdriven search can be viewed as an application of previous works [42,33,41,46] that use tracing to divide programs into subprograms with fixed support. We focus on the MAP inference problem where the posterior depends on an unknown cost function. We solve the problem through a learned cost-model-driven evolutionary search over traces and with validation.
Automatic neural program synthesis [12,18,11] has seen large progress recently. Alphacode [26] builds a system that can output creative and sound solutions to problems that require deeper-level human understanding. These approaches generate abstract syntax trees (ASTs) that can be incorrect and use input-output pairs to filter out those erroring programs. Our compiler approach requires us to ensure the correctness of all transformations. However, some ideas like validation after creation might be reusable. 6 Experiments

Expressiveness to Cover Common Optimization Techniques
This section aims to answer the following question: Is MetaSchedule expressive enough to capture the search space of the state-of-the-art optimization techniques? To answer this question, we evaluate our work on a diverse set of operators and subgraphs extracted from popular deep learning models, including variants of convolution, dense, and normalization.
As baselines, PyTorch (v1.11.0) results are provided to compare performance with vendor libraries; TVM (commit: 8d4f4dd73f), which incorporates AutoTVM [10] and Ansor [43], is used as the state-of-the-art tensor program optimization system, and we pick the best among the two in each respective setups. Full operators and hardware configurations are documented in Appendix A.2. Figure 8 shows that, in all cases on CPU and GPU, MetaSchedule delivers performance comparable with or even better than TVM, from which we could infer that MetaSchedule could express optimization techniques comparable to TVM on diverse workloads. Additionally, in most of the cases, MetaSchedule outperforms PyTorch by a significant margin except for SFM, which is highly optimized manually in PyTorch.

Optimizing End-to-End Deep Learning Models
Operator performance does not always translate directly to full model optimization. Therefore, this section is dedicated to answering the following question: Can MetaSchedule deliver competitive performance with state-of-the-art works for end-to-end models?
Therefore, a series of experiments are conducted to compare MetaSchedule and TVM, including BERT-Base [14], ResNet-50 [20], and MobileNet-v2 [36] on both CPU and GPU. As shown in Figure 9, MetaSchedule performance is on parity with TVM, while surpassing PyTorch in all cases, which indicates that MetaSchedule framework delivers end-to-end performance. Additionally, tuning time is provided in Appendix A.5.

Search Space Composition and Hardware-Specific Modules
Besides performance parity with existing work, in this section, we demonstrate the extra value of modular search space composition by answering the following question: How convenient is it to compose transformation modules, and how does it translate to performance?
We design an ablation study for transformation modules composition. As indicated in Figure 10a, by progressively enriching the search space, the performance of optimized tensor programs consistently increases. Composed with a hardware-specific module Use-Tensor-Core, MetaSchedule delivers significantly better performance compared with generic search space. The performance gain, brought by search space composition with customized rules, does translate to end-to-end model performance, as shown in Figure 10b. Specifically, on BERT-large workloads, MetaSchedule with Use-Tensor-Core delivers 48% speedup over TVM.
Notably, it took a graduate student only 2 days to craft the 82-line Use-Tensor-Core module in Python (see supplementary materials), which provides strong evidence of the convenience of customization and composition. More details are in Appendix A.4.

Conclusion
This paper presents MetaSchedule, a programming model to describe search space construction in tensor program optimization. Our method abstracts search space as a probabilistic language and enables flexible incorporation of domain knowledge by allowing practitioners to implement customized probabilistic programs. A learning-driven search algorithm is developed on top of the probabilistic language abstraction, which delivers competitive performance with state-of-the-art frameworks. In the future, we will explore and modularize declarative API for various hardware environments. Therefore, we will open-source our framework and hope it could enable broader collaboration between the machine learning deployment engineers and intelligent machine learning algorithms for tensor programs.

A.4 Search Space Construction and Customization for New Hardware
Our previous experience suggests that adapting to new hardware would require months of effort and deep expertise in the compilation stack, across code generation, transformations, and search across the legacy codebase.
Take TensorCore GPUs as an example. Unfortunately, Ansor in TVM does not support TensorCore, and it is highly non-trivial to extend its support to cover TensorCore without a dramatic revamp to its IR; Similarly, AutoTVM does not deliver comparable TensorCore performance without a dramatic re-design of its schedule tree, even though it supports a limited set of tensor intrinsics (e.g. WMMA).
In both cases, with the assumption of upgrading the core design, we anticipate that it will take months of engineering effort, and more to first get familiar with the codebase.
MetaSchedule enables doing that without such coupled complexity and no prerequisite experience on code generation or existing transformations, where each transformation is modeled as an independent primitive in the probabilistic language. Incorporating TensorCore, as an independent program transformation, could be engineered within 2 days by a grad student, and then composed into the existing system without having to re-design or affect any existing functionality.

A.5 Tuning Time
MetaSchedule makes an orthogonal contribution as it is a probabilistic language for composable search space construction rather than speeding up tuning. To provide a fair comparison of tuning speed, we reproduced Ansor's search space in our MetaSchedule probabilistic language in Table 1 TVM

A.6 End-to-End Integration Workflow with Deep Learning Framework
From frontend frameworks, for example, TensorFlow, PyTorch, or JAX, the tensor program to be optimized is generated from their computational graph. The generation process is generic to the shapes/ranks of tensors. The MetaSchedule program is generated from the tensor program obtained, based on which the search algorithm is performed. Sample an integer given the probability distribution sample-perfect-tile Sample the factors to perfect tile a specific loop sample-compute-location Sample a compute-at location of the given block