Time Series Data Cleaning: From Anomaly Detection to Anomaly Repairing (Technical Report)
Aoqian Zhang Shaoxu Song Jianmin Wang Philip S. Yu
cs.DB




4


Time Series Data Cleaning: From Anomaly Detection to Anomaly Repairing (Technical Report)
    
Aoqian Zhang
 

       School of Software, 
Tsinghua University

       zaq13@mails.tsinghua.edu.cn
Shaoxu Song
 

       School of Software, 
Tsinghua University

       sxsong@tsinghua.edu.cn
Jianmin Wang
 

       School of Software, 
Tsinghua University

       jimwang@tsinghua.edu.cn

Philip S. Yu
 

       University of Illinois at Chicago

       Institute for Data Science, Tsinghua University

       psyu@cs.uic.edu

    Final Version, December 30, 2023
=================================================================================================================================================================================================================================================================================================================================================================================================================================================




Errors are prevalent in time series data, such as GPS trajectories or sensor readings. 
Existing methods focus more on anomaly detection but not on repairing the detected anomalies. 
By simply filtering out the dirty data via anomaly detection, applications could still be unreliable over the incomplete time series. 
Instead of simply discarding anomalies, we propose to (iteratively) repair them in time series data, 
by creatively bonding the beauty of temporal nature in anomaly detection with 
the widely considered minimum change principle in data repairing.
Our major contributions include: 
(1) a novel framework of iterative minimum repairing (IMR) over time series data, 
(2) explicit analysis on convergence of the proposed iterative minimum repairing, 
and 
(3) efficient estimation of parameters in each iteration. 
Remarkably, with incremental computation, we reduce the complexity of parameter estimation from O(n) to O(1).
Experiments on real datasets demonstrate the superiority of our proposal compared to the state-of-the-art approaches.
In particular, we show that (the proposed) repairing indeed improves the time series classification application.





§ INTRODUCTION



Time series data are often found with dirty or imprecise values, such as GPS trajectories, sensor reading sequences <cit.>, or even stock prices <cit.>.
For example,  
the price of SALVEPAR (SY) is misused as the price of SYBASE (SY), both of which share the same notation (SY) in some sources. 
It is different from the interesting anomaly that actually happens in real life, e.g., 
the temperatures sudden change from 20C to 10C in one day when cold air rushes in.
To distinguish such cases, we propose to employ some labeled truth of dirty observations.
(See more detailed examples on dirty data and their labeled truth in Example <ref>.)



 §.§ Motivation on Anomaly Repairing



Applications, such as pattern mining <cit.> 
or classification <cit.>, 
built upon the dirty time series data are obviously not reliable. 
Anomaly detection over time series  is often applied to filter out the dirty data
(see <cit.> for a comprehensive and structured overview of anomaly detection techniques).
That is, the detected anomaly data points are simply discarded as useless noises. 
Unfortunately, with a large number of consecutive data points eliminated,
the applications could be barely performed over the rather incomplete time series. 

Recent study <cit.> shows that repairing dirty values could improve clustering over spatial data.
For time series data, we argue that repairing the anomaly can also improve the applications such as time series classification <cit.>. 
A repair close to the truth helps greatly the applications.





 §.§ Potential Methods for Repairing



A straightforward idea is to directly interpret the predication values in anomaly detection, 
e.g., by AR <cit.> or ARX <cit.>, 
as repairs (see details in Section <ref>). 
A data point is considered as anomaly if its (truth) predication significantly differs from (noisy) observation. 
Unfortunately, noisy/erroneous data are often close to the truth in practice, 
under the intuition that human or systems always try to minimize their mistakes, 
e.g., misspellings (John Smith vs. Jhon Smith), 
typos (555-8145 vs. 555-8195) 
as illustrated in <cit.>,
rounding off (76,821,000 vs. 76M) or unit error (76M vs. 76B) 
as shown in <cit.>.
Owing to such disagreement, the repairing performance of directly applying anomaly detection techniques is poor, as illustrated in both Example <ref> below and experiments in Section <ref>.


On the other hand, constraint-based repairing SCREEN <cit.>, strictly following the minimum change principle in data repairing <cit.>, 
heavily relies on a proper constraint of speeds on value changes. 
The repairing is performed based on two consecutive points, i.e., considering only one historical point, and thus 
does not sense the temporal nature of errors. 
As shown in the following Example <ref>, 
the speed constraint-based SCREEN is effective in repairing  spike errors, 
but can hardly handle a sequence of consecutive dirty points.

In short, the anomaly detection method does not expect the minimized mistakes in practice, 
whereas the constraint-based repairing is not effective in addressing the temporal nature of errors.




 §.§ Intuition of Our Proposal



Since completely automatic data repairing might not work well in  repairing time series data (such as SCREEN <cit.> observed in our experiments in Figure <ref>),
enlightened by the idea of utilizing master data (a single repository of high-quality data) in data repairing <cit.>, 
we employ some labeled truth of dirty observations to advance the repair.
The truth can be obtained either by manual labeling or automatically by more reliable sources.  
For instance, accurate locations are manually marked in the map by user check-in activities (and utilized to repair the imprecise GPS readings). 
Moreover, periodical automatic labeling may take place in certain scenarios, e.g., precise equipments report accurate air quality data  (as labeled values) in a relatively long sensing period, 
while crowd and participatory sensing generates unreliable observations in a constant manner <cit.>.

Being aware of both the error nature in anomaly detection and the minimum change principle in data repairing, 
we propose iterative minimum repairing (IMR). 
The philosophy behind the proposed iterative minimum repairing is that the high confidence repairs in the former iterations could help the latter repairing.
Specifically, IMR minimally changes one point a time to obtain the most confident repair only, referring to the minimum change principle in data repairing that human or systems always try to minimize their mistakes. 
The high confidence repairs, 
together with the labeled truth of error points, 
are utilized to learn and enhance the temporal nature of errors in anomaly detection, and thus generate more accurate repair candidates in the latter iterations. 






Figure <ref> presents an example segment of sensor readings, denoted by black line. 
Suppose that sensor errors occur in the period from time point 370 to 385, 
where the observations  are shifted from the truth, 
e.g., owing to granularity mismatch or unit error.
To repair the errors, the truth of several observations are labeled, including time points 
{370, 371, 372, 379, 387}.

Existing speed constraint-based cleaning (SCREEN) <cit.> could not effectively repair such continuous errors in a period 
(which is indeed also observed in Example 1 in <cit.>). 
The reason is that speed constraints, restricting the amount of value changing relative to time difference, can detect sharp deviations such as from time point 369 to 370, but not continuous errors, e.g., in 383 and 384.  
The exponentially weighted moving average (EWMA) <cit.> algorithm also hard to find a proper way to clean the trace. These two methods have similar repair trace.

By considering the predication values in anomaly detection as repairs (see more details in Section <ref>), the result of ARX based repairing is also reported. 
ARX, considering the errors between truth and observations, shows better repair results than EWMA and SCREEN methods. 
 
Finally, our proposed IMR approach, with both error predication and minimum change considerations, obtains repairs  closest to the truth.


The iterative minimum repairing leads to new challenges: 
(1) whether the repairing process converges; and 
(2) how to efficiently/incrementally update the parameter of the temporal model over the repaired data after each iteration.
Both issues in anomaly repairing are not considered in the  anomaly detection studies.



  
*Contributions


Our major contributions in this paper are summarized as follows.

(1) We formalize the anomaly repairing problem, given a time series with some points having labeled truth, in Section <ref>. 
The adaption of existing anomaly detection techniques (such as AR and ARX) is introduced for anomaly repairing.

(2) We devise an iterative minimum-change-aware repairing algorithm IMR, in Section <ref>. 
Remarkably, we illustrate that the ARX-based approach (in Section <ref>) is indeed a special case of IMR with static parameter (Proposition <ref>).

(3) We study the convergence of IMR in various scenarios, in Section <ref>. 
In particular, the convergence is explicitly analyzed for the special case of IMR(1) with order 𝑝=1, 
which is sufficient to achieve high repair accuracy in practice (as shown in the experiments in Section <ref>).
We prove that under certain inputs, the converged repair result could be directly calculated without iterative computing (Proposition <ref>).  

(4) We design efficient pruning and incremental computation for parameter estimation in each repair iteration, in Section <ref>. 
Rather than performing parameter estimation over all the n points, 
matrices for parameter estimation could be pruned by simply removing rows with value 0 (Proposition <ref>). 
It is also remarkable that the incremental computation among different repair iterations (Proposition <ref>) could further reduce the complexity of parameter estimation from 
O(n) to O(1).

(5) Experiments on real datasets with both real and synthetic errors, in Section <ref>, demonstrate that 
our IMR method shows significantly better repair performance than the state-of-the-art approaches, 
including the aforesaid anomaly detection and constraint-based repairing. 
Table <ref> 
lists the notations  frequently used  in this paper.






§ PRELIMINARIES


This section first introduces
the problem of anomaly repairing.  
We then adapt the existing anomaly detection models for anomaly repairing, 
i.e., AR without considering labeled data and ARX supporting labeled data.

The major issues of this simple adaption are:
(1) Applying predications with significant difference to the observation as repairs 
contradicts to the minimum change principle in data repairing <cit.>, 
as discussed in the introduction. 
(2) A static parameter (ϕ in Equations <ref> and <ref> below) needs to be preset, e.g., estimated from the dirty data during the initialization. 




 §.§ Problem Statement




Consider a time series of n  observations, 𝑥 = 𝑥[1],…,𝑥[n],  
where each 𝑥[i] is the value of the i-th data point. 
For brevity, we write 𝑥[i] as 𝑥_i.

Let 𝑦 denote the labeled/repaired sequence of 𝑥. 
Each 𝑦_i is either the labeled truth or the repaired value of 𝑥_i.

Given a time series 𝑥 and a 
partially labeled subset 𝑦 of 𝑥,
the repairing problem is to determine  the repairs 𝑦_i of 𝑥_i that are not labeled in 𝑦.


[Observation 𝑥,  partially labeled 𝑦, and fully repaired 𝑦]

Consider 𝑥 = {6, 10, 9.6, 8.3, 7.7, 5.4, 5.6, 5.9, 6.3,  6.8, 7.5, 8.5} 
with twelve data points of observations
in Figure <ref>, 
where shift (up) errors occur on four points x_2,…,x_5. 
Suppose that five points are labeled with truth, i.e., the partially labeled 𝑦.
By repairing (using the methods presented below), we propose to obtain a fully repaired 𝑦, e.g., 
𝑦 = {6, 5.6, 5.4, 5.2, 5.4, 5.4, 5.6, 5.9, 6.3, 6.8, 7.5, 8.5} as shown in Figure <ref>. 
In the repaired 𝑦, 𝑥_4 and 𝑥_5 are changed from 8.3 and 7.7 to 5.2 and 5.39, respectively. 
The labeled 𝑦_2 and 𝑦_3 will not be modified in the repair result.

Figure <ref> also presents another repair y' by the approach of connecting the dots with the labeled values, 
i.e., linear interpolation <cit.>. 
As shown in Figure <ref> and also indicated in <cit.>, the major issue of this (smoothing) approach is the serious damage of almost all the (unlabeled) data points, 
such as y'_7… y'_11, 
which are originally correct and should not be modified. 
In contrast, 
our proposed method repairs y_4 and y_5 while leaving y_7… y_11 unchanged. 





 §.§ AR Model



Intuitively, anomaly detection techniques could be adapted to  anomaly repairing. 
For instance, we consider the AR (autoregression) model <cit.> as follows: 

    𝑥'_t    = c + ∑_i=1^pϕ_i𝑥_t-i + ϵ_t

where 
𝑥'_t is the prediction of 𝑥_t, 
𝑝 is the order,  
ϕ_i is the parameter of the model, 
c is a constant defined by 
c = μ(1-∑_i=1^pϕ_i), 
μ is mean value of the process and 
ϵ_t is white noise
(usually Gaussian white noise <cit.>, a normal random variable generated according to the Gaussian distribution with mean μ=0 and variance σ^2; in other words, c=0).

If 𝑥'_t significantly differs from the original observation 𝑥_t, having 
|𝑥'_t-𝑥_t|>τ where τ is a predefined threshold, 
this predication is accepted 𝑥_t=𝑥'_t, a.k.a. a repair. 
The intuition behind is that a farther distance indicates the higher probability of being an outlier.
The threshold τ can be decided by observing the statistical distribution of distances between 𝑥_𝑡' and 𝑥_
𝑡,  using the prediction interval <cit.>.

The AR-based repairing procedure is thus: 
(1) replace 𝑥_t by 𝑦_t if it is labeled,
(2) learn parameter ϕ of AR(𝑝) from 𝑥, and 
(3) fill all unlabeled 𝑦_t by AR(𝑝) over 𝑥, having 

    𝑦_t = 
    𝑥'_t      if 𝑦_t  is unlabeled and  |𝑥'_t-𝑥_t|>τ
    𝑥_t       otherwise





[Example <ref> continued]

Consider again 𝑥 = {6, 10, 9.6, 8.3, 7.7, 5.4, 5.6, 5.9, 6.3, 6.8, 7.5, 8.5}
in Figure <ref>. 
For simplicity, we use AR(1) with order 𝑝 = 1 and c = 0, i.e., 𝑥_t' = ϕ_1𝑥_t-1. 
By ordinary least square <cit.>, we estimate the parameter ϕ from 𝑥, having ϕ_1 = 1.022. 
Let τ = 0.1. 
𝑦_1 is labeled with truth. 
Referring to Equation <ref>, it outputs unchanged 𝑦_1  = 6.
Similarly for 𝑦_2 and 𝑦_3. 
We have 𝑥_4' = ϕ_1𝑥_1 = 1.022 * 5.4 = 5.52. 
Since |5.52 - 8.3| = 2.78 > 0.1, 𝑥_4' is accepted as new 𝑥_4.

Similarly, we have 𝑥_5' = ϕ_1𝑥_4 = 1.022 * 5.52 = 5.64. 
Referring to |5.64 - 7.7| = 2.06 > 0.1, the prediction is accepted.
So on and so forth, we obtain the final repaired result 
𝑦 = {6, 5.6, 5.4, 5.52, 5.64, 5.4, 5.6, 5.72, 5.84, 5.97, 6.10, 8.5}. 
Its RMS error is 0.51 (see Section <ref> for RMS definition).





 §.§ ARX Model



In order to utilize the labeled 𝑦, we consider the ARX model (autoregressive model with exogenous inputs) <cit.> 

    𝑦'_t = 𝑥_t + ∑_i=1^pϕ_i(𝑦_t-i-𝑥_t-i) + ϵ_t

where 𝑦'_t is the possible repair of 𝑥_t,
and others are the same to the aforesaid AR model. 
As shown in Equation <ref>, not only the preceding observations 
𝑥_t-i will affect the determination of 𝑦'_t, 
but also the previously labeled/repaired 𝑦_t-i.

The ARX-based repairing procedure is thus: 
(1) learn parameter ϕ of ARX(𝑝) from 𝑥 and partially labeled 𝑦, and 
(2) fill all unlabeled 𝑦_t by ARX(𝑝), similar to Equation <ref> by replacing 𝑥'_t with 𝑦'_t.


[Example <ref> continued]

Consider  again 𝑥 = {6, 10, 9.6, 8.3, 7.7, 5.4, 5.6, 5.9, 6.3, 6.8, 7.5, 8.5} and the partially labeled 𝑦 over five points, in Figure <ref>.
For simplicity, we also use ARX(1) with order 𝑝 = 1, i.e., 𝑦_t' = 𝑥_t + ϕ_1(𝑦_t-1-𝑥_t-1). 
Similar to AR in Example <ref>, we estimate the parameter ϕ by ordinary least square <cit.>, having ϕ_1 = 0.5. 
Let τ = 0.1.
Again, the labeled 𝑦_3=5.4 is not modified. 
For the fourth point, we have 𝑦_4' = 8.3 + 0.5*(5.4-9.6) = 6.2. Since |6.2 - 8.3| = 2.1 > 0.1, 
we assign 𝑦_4 = 𝑦_4' = 6.2.
Finally, the  repair result by ARX is 𝑦 = {6, 5.6, 5.4, 6.20, 6.65, 5.4, 5.6, 5.9, 6.3, 6.8, 7.5, 8.5} 
with RMS error 0.49, lower than that of AR in Example <ref>.



We consider ARX model since it can capture the difference between the observed errors and labeled truths, 
while other methods such as AR and SCREEN ignore (cannot utilize) such differences.
By modeling such differences between errors and truths rather than original values, 
both ARX and our proposed IMR may not deal with Spike errors (i.e., with # consecutive errors = 1 in Figure <ref>) as good as SCREEN <cit.>.
Nevertheless, with iterative repairing, 
IMR always shows significantly better results in addressing a large number of consecutive errors (see Figure <ref> as well). 
ARX cannot address well such consecutive errors either, since it applies only significant changes which is contradict to the minimum change principle in data repairing as discussed in Section <ref> in the Introduction.





§ REPAIR ALGORITHM


Unlike the existing anomaly detection model that may over-change the data via one-pass repairing (as discussed in the introduction, illustrated in Figure <ref>, and observed in experiments in Section <ref>),
we propose to progressively repair the data, being aware of both the error nature in anomaly detection and the minimum change principle in data repairing, 
so that the high confidence repairs in the former iterations could help the repairing in the latter steps.



 §.§ Iterative Repairing




Let 𝑦^(k) denote the sequence 𝑦 in the k-th iteration, 
where 𝑦^(0) is the partially labeled time series in the input. 
Since 𝑦^(0) is incomplete (partially labeled), 
to initialize, 
we assign 𝑦_t^(0)=𝑥_t if 𝑦_t^(0) is not labeled. 
Recall that the labeled values should not be repaired, i.e., 
𝑦^(k)_t=𝑦^(0)_t if 𝑦^(0)_t is labeled. 
 
Algorithm <ref> presents the iterative minimum repairing procedure, IMR(𝑝),
whose inputs are the observation time series 𝑥 and partially labeled 𝑦^(0). 
It outputs 𝑦^(k) with all the labeled 𝑦^(0)_t unchanged and unlabeled 𝑦^(0)_t repaired. 





The major steps include: 

(S1) Parameter estimation, in Line <ref>, 
learns the parameter of ARX(p) in the k-th iteration, denoted by ϕ^(k), from 𝑥  and the current 𝑦^(k).

(S2) Repair candidate generation, in Line <ref>, 
computes the possible repairs 𝑦̂^(k), according to ARX(𝑝) 
w.r.t. 𝑥, 𝑦^(k) and ϕ^(k).


(S3) Repair evaluation, 
in Line <ref>, 
determines one of the repairs to accept, 
𝑦_t^(k+1) = 𝑦̂_t^(k), 
referring to the minimum change principle in data repairing <cit.>. 


As shown in Line <ref>, the procedure repeats, until the repair converges, 
e.g., having 

    |𝑦_j^(k) - 𝑦_j^(k+1)| ≤τ , 
    j = 1, …,  n.

where τ a threshold of convergence,
or a maximum number of iterations is reached.
Setting max-num-iterations is a remedy to avoid waiting for convergence in practice (see Section <ref> for discussion and evaluation). 

 

[Algorithm overview, Example <ref> continued]

Consider again 𝑥 = {6, 10, 9.6, 8.3, 7.7, 5.4, 5.6, 5.9, 6.3, 6.8,  7.5, 8.5}
in Figure <ref>. 
According to five labeled data points, 
we assign 𝑦^(0) = {6, 5.6, 5.4, 8.3, 7.7, 5.4, 5.6, 5.9, 6.3, 6.8, 7.5, 8.5}, 
where the unlabeled points are initialized by 𝑦_t^(0) = 𝑥_t,  
e.g., 𝑦_4^(0) = 𝑥_4 = 8.3. 

In each iteration, the IMR algorithm 
(1) learns the parameter, e.g., ϕ_1^(0)=0.5 for 𝑝=1;
(2) generates candidates for repairing, such as 𝑦̂^(0)  = {-, -, -, 6.2, 7.7, -, 5.6, 5.9, 6.3, 6.8,  7.5, -}; and
(3) selects one repair to conduct, and form the new sequence, 
say 𝑦^(1) = {6, 5.6, 5.4, 6.2, 7.7, 5.4, 5.6, 5.9, 6.3, 6.8, 7.5, 8.5}.

The procedure repeats until converging.
The final output is 𝑦^(7) = {6, 5.6, 5.4, 5.20, 5.39, 5.4, 5.6, 5.9, 6.3, 6.8, 7.5, 8.5} with RMS error 0.03.
Details on each step are presented in the following examples.




 §.§ Parameter Estimation



The parameter estimation step S1 (in Line <ref> in Algorithm <ref>) 
estimates the parameter ϕ^(k) for ARX(𝑝), 
given 𝑥, 𝑦^(k). 
Existing methods such as 
Ordinary Least Square <cit.>
or 
Yule-Walker Equations <cit.>
can be directly employed. 
For instance, by Ordinary Least Square, 
we have 

    ϕ^(k)   = ((𝑍^(k))'𝑍^(k))^-1(𝑍^(k))'𝑉^(k)

where 

    𝑉^(k) =    [ 𝑦_𝑝 + 1^(k) - 𝑥_𝑝 + 1; 𝑦_𝑝 + 2^(k) - 𝑥_𝑝 + 2;                     ⋮;         𝑦_n^(k) - 𝑥_n ]
    ,     ϕ^(k) = 
    [ ϕ_1^(k); ϕ_2^(k);       ⋮; ϕ_𝑝^(k) ] 
    , 
    𝑍^(k) =    [     𝑦_𝑝^(k) - 𝑥_𝑝 𝑦_𝑝-1^(k) - 𝑥_𝑝-1                 …     𝑦_1^(k) - 𝑥_1; 𝑦_𝑝+1^(k) - 𝑥_𝑝+1     𝑦_𝑝^(k) - 𝑥_𝑝                 …     𝑦_2^(k) - 𝑥_2;                 ⋮                 ⋮                 ⋱                 ⋮; 𝑦_n-1^(k) - 𝑥_n-1 𝑦_n-2^(k) - 𝑥_n-2                 … 𝑦_n-𝑝^(k) - 𝑥_n-𝑝 ].




[Parameter estimation on 𝑦^(0), Example <ref> continued]
Consider 𝑥 = {6, 10, 9.6, 8.3, 7.7, 5.4, 5.6, 5.9, 6.3, 6.8,  7.5, 8.5} and 𝑦^(0) = {6, 5.6, 5.4, 8.3, 7.7, 5.4, 5.6, 5.9, 6.3, 6.8,  7.5, 8.5}. 
Given order 𝑝=1, we have 𝑉^(0) = {-4.4,-4.2,0, 0,0, 0,0,0,0,0,0}' 
with 11 rows and 1 column,
and 𝑍^(0) = {0,-4.4,  -4.2,0,0,0,0,0,0,0,0}'
with 11 rows and 1 column.
Referring to Equation <ref>,
the parameter is estimated by

    ϕ_1^(0) = (-4.4)*(-4.2)/(-4.4)^2+(-4.2)^2 = 0.5.



Owing to the iterative repairing,  online incremental parameter estimation is necessary, which is not studied in the previous studies (see our approach in Section <ref>).



 §.§ Candidate Generation



The repair candidate generation step S2 (in Line <ref> in Algorithm <ref>)
employs ARX(𝑝)  to infer the candidate repair 
𝑦̂^(k) = ϕ^(k)·(𝑦^(k) - 𝑥) + 𝑥,
referring to the estimated parameter ϕ^(k).  
More specifically,
for each point t, 
𝑦̂_t^(k) is given by  

    𝑦̂_t^(k) = ∑_i=1^pϕ_i^(k)(𝑦_t-i^(k)-𝑥_t-i) + 𝑥_t

according to 𝑦_t-1^(k), …, 𝑦_t-p^(k). 
We note that only candidates with |𝑦̂_t^(k) - 𝑦_t^(k)| > τ need to be considered 
referring to the convergence condition in Equation <ref>.
 

[Repair candidate 𝑦̂^(0), Example <ref> continued]
Consider the parameter ϕ_1^(0) = 0.5 estimated in Example <ref>. 
Let threshold τ = 0.1. 
Referring to Equation <ref>, 
we have 
𝑦̂_4^(0) = 0.5 * (5.4 - 9.6) + 8.3 = 6.2
with |𝑦̂_4^(0) - 𝑦_4^(0)| = |6.2 - 8.3| = 2.1 > 0.1, 
and 
𝑦̂_5^(0) = 0.5 * (8.3 - 8.3) + 7.7 = 7.7 
with |7.7 - 7.7| = 0 < 0.1.
The repair candidates are 𝑦̂^(0)  = {+, +, +, 6.2, -, +, -, -, -, -, -+} 
where 
`+' corresponds to the labeled points and 
`-' denotes no candidates. 
That is, we need to consider only one candidate 𝑦̂_4^(0) for repairing.





 §.§ Repair Evaluation



The repair evaluation step S3  
(in Line <ref> in Algorithm <ref>)
selects one repair to accept, 
i.e., assigning 𝑦_t^(k+1) = 𝑦̂_t^(k)
the aforesaid generated repair candidate.
Following the minimum change principle in data repairing <cit.>, 
the repair that minimally differs from its original input is preferred with higher confidence.
The repaired result in each iteration is: 

    𝑦_t^(k+1) = 
      𝑦̂_t^(k)     if  t = _i |𝑦̂_i^(k) - 𝑥_i| 
    𝑦_t^(k)     otherwise  .


Remarkably, only one data point with the minimum change (most confident) is repaired in each iteration, which is more efficient than the NP-hard problem of minimizing the overall changes w.r.t. integrity constraints <cit.>.

[Minimum repair 𝑦̂_t^(1), Example <ref> continued]

Since there is only one repair candidate obtained in Example <ref>, 
i.e.,  𝑦̂_4^(0) = 6.2, it is the minimum repair (among all candidates).
The sequence after the first iteration becomes 𝑦^(1) = {6, 5.6, 5.4, 6.2, 7.7, 5.4, 5.6, 5.9, 6.3, 6.8, 7.5, 8.5}.


Note that the minimum change principle <cit.> in data repairing is based on the intuition that human or systems always try to minimize their mistakes. 
However, it is not guaranteed that the minimum change repair always corresponds to the true value. 
Therefore, similar to other minimum change-based data repairing studies <cit.>, the accuracy of the final results is unlikely to have theoretical guarantees, since there is no constraint on how far the errors may diverge from the truth.
For this reason, we can only evaluate the correctness of the proposed repair by comparing to the ground truth in experiments, similar to other data repairing studies <cit.> as well. 
Nevertheless, 
we can show that the efficient pruning and incremental computation are safe (Propositions <ref> and <ref>), 
i.e., the accuracy of the final results with efficient computing is theoretically guaranteed to be the same as the results of original IMR without pruning and incremental computation.




§ CONVERGENCE ANALYSIS


In this section, we analyze the convergence of iterative repairing, i.e., 
lim_k→+∞∑_i=1^n (
𝑦_i^(k+1) - 𝑦_i^(k))
 = 0,
which is essential to the termination of Algorithm <ref>. 
While the general convergence problem is still open, we study the convergence of the proposed method in certain special cases in this section for two reasons: 
(1) we illustrate that the ARX-based approach is indeed a special case of the proposed IMR with static parameter (Proposition <ref>) in Section <ref>; 
(2) we identify another special case with convergence guarantee in Section <ref>, which enables online repairing over streaming data without iteration 
(also see Sections <ref> and <ref> for more details and experiments).




 §.§ Static Parameter




We study this special case in order to illustrate the relationship between our proposed IMR and the existing ARX.
Let us first analyze the convergence of IMR (Proposition <ref>) 
and then illustrate their equivalence in certain case (Proposition <ref>). 

Rather than dynamically updating the parameter ϕ^(k) in each iteration, 
in Line <ref> in Algorithm <ref>, 
a special case is to specify a static parameter, ϕ^(k)=ϕ^(0), for all the iterations.





With a static parameter ϕ^(k) = ϕ, ∀ k, 
the repair result converges, 
i.e., 

    lim_k→+∞∑_i=1^n (
    𝑦_i^(k+1) - 𝑦_i^(k))
     = 0.




Let t be the point repaired in iteration k. 
We show that 𝑦_t^(k+i)=𝑦_t^(k+1) converges (will not change in the following iterations k+i), 
if the 𝑝 preceding points of t are converged (no longer change). 
The reason is that 𝑦̂_̂t̂^(k+i)=𝑦_t^(k+1),
referring to Equation <ref> of candidate generation 
with static parameter ϕ.
That is, no repair candidate will be generated.

If one of the 𝑝 preceding points of t changes (not converged), say t', 
the aforesaid derivation applies similarly to t'. 
By recursively applying the derivations, we can finally find a point t^* whose 𝑝 preceding points are either labeled or converged, i.e., a new point t^* converges. 
Such t^* always exists, i.e., 𝑝+1, given that the first 𝑝 points in the sequence will not be modified by the repairing 
(referring to the repair candidate generation in Equation <ref>) 
and thus already converged under the static parameter ϕ.


Finally, all the data points gradually converge. 






  §.§.§ Special Case of IMR(1)


We show in the following that
for 𝑝=1,  
the ARX(𝑝)-based repairing in Section <ref> 
is a special case of our proposed IMR(𝑝) with static parameter ϕ^(0).
This equivalence demonstrates the rationale of our proposal. 





For IMR(1) with static parameter ϕ^(k) = ϕ, ∀ k, 
Algorithm <ref> 
is equivalent to 
ARX(1)-based repairing.



Suppose that the first point is labeled. 
Consider data point 2. 
Referring to Proposition <ref>, 
𝑦_2^(k_2) becomes converged in some iteration k_2, 
denoted by 𝑦_2=ϕ_1(𝑦_1-𝑥_1)+𝑥_2, 
where 𝑦_1=𝑦_1^(0) is labeled. 

Similarly, given the converged 𝑦_2, 
𝑦_3^(k_3) will converge in some iteration k_3, 
denoted by 𝑦_3=ϕ_1(𝑦_2-𝑥_2)+𝑥_3. 

By recursively obtaining the converged  𝑦_t which is not labeled, 
it is exactly the procedure of ARX(1)-based repairing. 





 §.§ Converged Parameter





We now consider the dynamically updated parameter ϕ^(k) in each iteration, 
in Line <ref> in Algorithm <ref>. 
As shown in the following Proposition <ref>, if the dynamic parameter converges, 
the repair converges as well. 
This converged parameter case is interesting, since 
the corresponding converged repair result could be directly calculated without iterative computing in certain cases
as illustrated below. 



If the parameter converges, 
lim_k→+∞ϕ^(k) = ϕ, 
then the repair also converges

    lim_k→+∞∑_i=1^n (
    𝑦_i^(k+1) - 𝑦_i^(k))
     = 0.




When the parameter ϕ converges, it becomes static. 
Similar proof for Proposition <ref> applies.





  §.§.§ Special Case of IMR(1)



Again, we consider the special case of IMR(1) with order 𝑝=1. 
To show how the repair results could be directly computed without iterations, we first illustrate that any 𝑦_t^(k) generated during Algorithm <ref> can be represented as follows, 
a.k.a. provenance of 𝑦_t^(k). 



For IMR(1), we can represent each 𝑦_t^(k) by 

    𝑦_t^(k) = ϕ_1^(k_s)ϕ_1^(k_s-1)…ϕ_1^(k_1)(𝑦_t-s^(0)-𝑥_t-s) + 𝑥_t,

where 0<k_1<…<k_s-1<k_s<k denote iteration numbers, 
𝑦_t-s^(0) is labeled truth, 
and 
time points t-s+1, …, t are not labeled. 



Referring to Equation <ref>, 𝑦_t^(k) in k-th iteration is either unchanged (equal to 𝑦_t^(k-1))
or determined by 𝑦̂_t^(k-1)=ϕ_1^(k-1)(𝑦_t-1^(k-1)-𝑥_t-1) + 𝑥_t referring to Equation <ref> with 𝑝=1.

Suppose that 𝑦_t^(k) is unchanged since iteration k_s.
We have 

    𝑦_t^(k) = 𝑦_t^(k-1) = … = 𝑦_t^(k_s+1) 
    = ϕ_1^(k_s)(𝑦_t-1^(k_s)-𝑥_t-1) + 𝑥_t 
    .

Similar provenance applies to 𝑦_t-1^(k_s), …, 𝑦_t-1^(k_1), 

    𝑦_t-1^(k_s) 
    =    … = 𝑦_t-1^(k_s-1+1) 
    = ϕ_1^(k_s-1)(𝑦_t-2^(k_s-1)-𝑥_t-2) + 𝑥_t-1 
    , 
       … ,
    𝑦_t-s+1^(k_2) 
    =    … = 𝑦_t-s+1^(k_1+1) 
    = ϕ_1^(k_1)(𝑦_t-s^(k_1)-𝑥_t-s) + 𝑥_t-s+1 
    , 
    𝑦_t-s^(k_1) 
    =    … = 𝑦_t-s^(0)
    .

Combining these s+1 derivations, the conclusion is proved. 



We denote the labeled points in 𝑦^(0) by multiple (say m) segments. 
Let s(j) and e(j) denote the start and end point of the j-th labeled segment, j=1,…,m.
For instance, there are 3 segments of labeled data points in Figure <ref>, 
having s(1)=1, e(1)=2, 
s(2)=4, e(2)=5, 
s(3)=9, e(3)=9.

Let 𝑧_i = 𝑦_i^(0) - 𝑥_i for all labeled points i.
(We set 𝑧_e(0)=0 over undefined segment 0.)





For IMR(1), if the parameter converges, having 
lim_k→+∞ϕ_1^(k) = ϕ_1,
then the converged repair result can be directly given by

    lim_k→+∞𝑦_i^(k) = 𝑦_i,
 
where

    𝑦_i =  
    𝑦_i^(0)   if  i ∈ [s(j),e(j)] 
    ϕ_1^i-e(j)(𝑦_e(j)^(0) - 𝑥_e(j)) + 𝑥_i
       if  i ∈ (e(j), s(j+1))

and the converged parameter ϕ_1 is a solution to 

    ϕ_1 = ∑_j=1^m(
    ϕ_1^s(j)-1-e(j-1)𝑧_e(j-1)𝑧_s(j)
    + 
    ∑_i=s(j)^e(j)-1𝑧_i𝑧_i+1)
    /∑_j=1^m(
    (ϕ_1^s(j)-1-e(j-1)𝑧_e(j-1))^2
    + 
    ∑_i=s(j)^e(j)-1(𝑧_i)^2 
    )
    .


 
Given lim_k→+∞ϕ_1^(k) = ϕ_1,
according to Proposition <ref>, the result 𝑦^(k) will converge, having lim_k→+∞𝑦_i^(k) = 𝑦_i.

Referring to Lemma <ref>,
the converged value has its provenance from the nearest preceding labeled point, i.e., 𝑦_i = ϕ_1^i-e(j)(𝑦_e(j)^(0) - 𝑥_e(j)) + 𝑥_i, for e(j) < i < s(j+1).

Let 𝑧_t = 𝑦_t - 𝑥_t for all unlabeled points t.
Referring to Equation <ref> of parameter estimation, we have

    ϕ_1    = 𝑧_1𝑧_2
    +…+
    𝑧_n-1𝑧_n/
    (𝑧_1)^2
    +…+
    (𝑧_n-1)^2
    .


Similar to the aforesaid 𝑦_i, we have  
𝑦_i+1= ϕ_1^i-e(j)+1(𝑦_e(j)^(0) - 𝑥_e(j)) + 𝑥_i+1, for e(j) < i+1 < s(j+1).
It follows 
𝑧_i+1=ϕ_1𝑧_i.

By applying 𝑧_i+1=ϕ_1𝑧_i for all e(j) ≤ i < s(j+1)-1, j ∈ [1, m], we have 

    ϕ_1    = 𝑧_1𝑧_2
    +…+
    𝑧_n-1𝑧_n/
    (𝑧_1)^2
    +…+
    (𝑧_n-1)^2
    
       = ∑_j=1^m(
    ∑_i=e(j-1)^s(j)-1𝑧_i𝑧_i+1 
    + 
    ∑_i=s(j)^e(j)-1𝑧_i𝑧_i+1)
    /∑_j=1^m(
    ∑_i=e(j-1)^s(j)-1(𝑧_i)^2
    + 
    ∑_i=s(j)^e(j)-1(𝑧_i)^2
    )
    
       = ∑_j=1^m(
    ∑_i=e(j-1)^s(j)-2ϕ_1(𝑧_i)^2 
    + 
    𝑧_s(j)-1𝑧_s(j)
    +
    ∑_i=s(j)^e(j)-1𝑧_i𝑧_i+1)
    /∑_j=1^m(
    ∑_i=e(j-1)^s(j)-1(𝑧_i)^2
    + 
    ∑_i=s(j)^e(j)-1(𝑧_i)^2
    )
    
       = ∑_j=1^m(
    𝑧_s(j)-1𝑧_s(j)
    + 
    ∑_i=s(j)^e(j)-1𝑧_i𝑧_i+1)
    /∑_j=1^m(
    (𝑧_s(j)-1)^2 
    + ∑_i=s(j)^e(j)-1(𝑧_i)^2
    )


By applying Equation <ref> again, i.e., 

    𝑧_s(j)-1=ϕ_1^s(j)-1-e(j-1)𝑧_e(j-1),
 
the conclusion in Equation <ref> is proved.







  §.§.§ IMR(1) with One Labeled Segment



We consider the case that only one segment with length ℓ is labeled at the beginning of 𝑦^(0), 
i.e., 𝑦^(0)_1, 𝑦^(0)_2,  …,  𝑦^(0)_ℓ are labeled. 
In this special case, the converged parameter and repair result can be directly calculated  without iterating, 
and most importantly 
it enables efficient online computation, 
by interpreting all the historical data as one segment labeled (see Section <ref> for details and evaluation). 
Remarkably, no threshold needs to be set in this case.

The idea is: 
(1) We first show in Lemma <ref> that under certain inputs, the estimated parameter in each iteration is indeed bounded;
(2) Proposition <ref> then illustrates that the bounded parameter leads to converged parameter;
(3) Finally, analogous to Proposition <ref>, given the converged parameter, 
we directly calculate the converged repair without iterative computing in Proposition <ref>.



For IMR(1) with first ℓ data points labeled in 𝑦^(0).
If the input satisfies
|∑_t=1^ℓ-1𝑧_t^(0)𝑧_t+1^(0)|
< ∑_t=1^ℓ-1𝑧_t^(0)𝑧_t^(0), i.e., 

    |∑_t=1^ℓ-1(𝑦_t^(0)-𝑥_t)(𝑦_t+1^(0)-𝑥_t+1)|
    < ∑_t=1^ℓ-1(𝑦_t^(0)-𝑥_t)^2,

then we have |ϕ_1^(k)| < 1 in the iterations k, 0 ≤ k ≤ n-ℓ.




According to Equation <ref>, 
we have 

    ϕ_1^(0) = ∑_t=1^ℓ-1𝑧_t^(0)𝑧_t+1^(0)/
    (𝑧_ℓ^(0))^2+∑_t=1^ℓ-1(𝑧_t^(0))^2
    .
 
Referring to the given condition, it follows

    |ϕ_1^(0)| ≤
    |∑_t=1^ℓ-1𝑧_t^(0)𝑧_t+1^(0)|
    /∑_t=1^ℓ-1(𝑧_t^(0))^2
    
    <1.


We prove |ϕ_1^(k)| < 1, 1 ≤ k ≤ n-ℓ in three cases.




  
*Case 1: ϕ_1^(0) > 0

We will prove 0<ϕ_1^(k)< 1, 1 ≤ k ≤ n-ℓ by induction. 

Basis:
For k = 1, 
referring to the minimal change principle and the candidate generation in Section <ref>, 
we have _i ∈ [1,n]|𝑧̂_i^(0)|=ℓ+1 and  
𝑧_ℓ+1^(1) = 𝑧̂_ℓ+1^(0) = 𝑧_ℓ^(0)ϕ_1^(0). Since only one point will be changed in each iteration, we have
𝑧_i^(1) = 𝑧_i^(0), i ≠ℓ+1.
Referring to Equation <ref>, 
it follows

    ϕ_1^(1)   = 𝑧_ℓ+1^(1)𝑧_ℓ^(1) + ∑_t=1^ℓ-1𝑧_t^(0)𝑧_t+1^(0)/
    (𝑧_ℓ+1^(1))^2 + ∑_t=1^ℓ(𝑧_t^(0))^2


The following derivation shows ϕ_1^(1) > 0.

    ϕ_1^(1)-ϕ_1^(0)   = 
    𝑧_ℓ+1^(1)𝑧_ℓ^(1) + ∑_t=1^ℓ-1𝑧_t^(0)𝑧_t+1^(0)/
    (𝑧_ℓ+1^(1))^2 + ∑_t=1^ℓ(𝑧_t^(0))^2
     - ∑_t=1^ℓ-1𝑧_t^(0)𝑧_t+1^(0)/∑_t=1^ℓ(𝑧_t^(0))^2
    
       = 𝑧_ℓ+1^(1)(𝑧_ℓ^(1)(∑_t=1^ℓ(𝑧_t^(0))^2)
    -𝑧_ℓ+1^(1)(∑_t=1^ℓ-1𝑧_t^(0)𝑧_t+1^(0)))
    /
    (∑_t=1^ℓ(𝑧_t^(0))^2)((𝑧_ℓ+1^(1))^2 + ∑_t=1^ℓ(𝑧_t^(0))^2)
    
       = ϕ_1^(0)(1-(ϕ_1^(0))^2)(𝑧_ℓ^(0))^2
    /
    (𝑧_ℓ+1^(1))^2 + ∑_t=1^ℓ(𝑧_t^(0))^2
     > 0


And it follows ϕ_1^(1) < 1 below.

    1-ϕ_1^(1)   = 
    (∑_t=1^ℓ(𝑧_t^(0))^2-∑_t=1^ℓ-1𝑧_t^(0)𝑧_t+1^(0)) 
    +𝑧_ℓ+1^(1)(𝑧_ℓ+1^(1)-𝑧_ℓ^(1))
    /
    (𝑧_ℓ+1^(1))^2 + ∑_t=1^ℓ(𝑧_t^(0))^2
    
       = 
    (∑_t=1^ℓ(𝑧_t^(0))^2-∑_t=1^ℓ-1𝑧_t^(0)𝑧_t+1^(0)) 
    +(𝑧_ℓ^(0))^2(ϕ_1^(0)(ϕ_1^(0)-1))
    /
    (𝑧_ℓ+1^(1))^2 + ∑_t=1^ℓ(𝑧_t^(0))^2
    
       > 
    (𝑧_ℓ^(0))^2(1+(ϕ_1^(0))^2 - ϕ_1^(0))
    /
    (𝑧_ℓ+1^(1))^2 + ∑_t=1^ℓ(𝑧_t^(0))^2
     > 0


Combining the aforesaid two derivations, we have 0 < ϕ_1^(0) < ϕ_1^(1) < 1.

Inductive step:
We will show that if the conclusion holds 0 < ϕ_1^(k) < 1, 
for k = i-1, 1 ≤ i ≤ n-ℓ, 
then also 0 < ϕ_1^(i) < 1 holds for k=i.



Following the IMR, we have
_j ∈ [1,n]|𝑧̂_j^(i-1)|
=ℓ+i and  
𝑧_ℓ+i^(i) = 𝑧̂_ℓ+i^(i-1) = 𝑧_ℓ+i-1^(i-1)ϕ_1^(i-1).

The following shows ϕ_1^(i) > 0.

    ϕ_1^(i)-ϕ_1^(i-1)
       = 𝑧_ℓ+i^(i)𝑧_ℓ+i-1^(i) + ∑_t=1^ℓ+i-2𝑧_t^(i-1)𝑧_t+1^(i-1)/
    (𝑧_ℓ+i^(i))^2 + ∑_t=1^ℓ+i-1(𝑧_t^(i-1))^2
     - ∑_t=1^ℓ+i-2𝑧_t^(i-1)𝑧_t+1^(i-1)/∑_t=1^ℓ+i-1(𝑧_t^(i-1))^2
    
       = 𝑧_ℓ+i^(i)
    (𝑧_ℓ+i-1^(i)(∑_t=1^ℓ+i-1(𝑧_t^(i-1))^2)
    -𝑧_ℓ+i^(i)(∑_t=1^ℓ+i-2𝑧_t^(i-1)𝑧_t+1^(i-1)))
    /
    (∑_t=1^ℓ+i-1(𝑧_t^(i-1))^2)((𝑧_ℓ+i^(i))^2 + ∑_t=1^ℓ+i-1(𝑧_t^(i-1))^2)
    
       = ϕ_1^(i-1)(1-(ϕ_1^(i-1))^2)(𝑧_ℓ+i-1^(i-1))^2
    /
    (𝑧_ℓ+i^(i))^2 + ∑_t=1^ℓ+i-1(𝑧_t^(i-1))^2
     > 0


It follows ϕ_1^(i) < 1 as illustrated below.

    (1-ϕ_1^(i))( 
    (𝑧_ℓ+i^(i))^2 + ∑_t=1^ℓ+i-1(𝑧_t^(i-1))^2
    ) 
       = (∑_t=1^ℓ(𝑧_t^(0))^2-∑_t=1^ℓ-1𝑧_t^(0)𝑧_t+1^(0)) 
       +(𝑧_ℓ^(0))^2 ×∑_s=1^i∏_j=0^s-2(ϕ_1^(j))^2ϕ_1^(s-1)(ϕ_1^(s-1)-1) 
       > (𝑧_ℓ^(0))^2 ×
    (1 + (ϕ_1^(i-1)-1)(∑_s=1^i(ϕ_1^(i-1))^2s-1)) 
       = (𝑧_ℓ^(0))^2 ×
    (1+ (ϕ_1^(i-1)-1)ϕ_1^(i-1)(1-(ϕ_1^(i-1))^2i)
    /
    1-(ϕ_1^(i-1))^2
    ) 
       = (𝑧_ℓ^(0))^2 ×
    (1 + (ϕ_1^(i-1)-1)ϕ_1^(i-1)-(ϕ_1^(i-1))^2i+1/
    (1-ϕ_1^(i-1))(1+ϕ_1^(i-1))
    ) 
       = (𝑧_ℓ^(0))^2 ×
    1+(ϕ_1^(i-1))^2i+1/
    1+ϕ_1^(i-1) > 0


Thereby 0 < ϕ_1^(k) < 1 holds for k=i.




  
*Case 2: ϕ_1^(0) < 0
Similar proof steps apply.



  
*Case 3: ϕ_1^(0) = 0
We can show that ϕ_1^(k) = 0, 1≤ k ≤ n-ℓ.


The following conclusion illustrates that with a bounded parameter in each iteration,
the parameter converges. 



For IMR(1) with first ℓ data points labeled in 𝑦^(0), if |ϕ_1^(k)| < 1 in the iterations k, 0 ≤ k ≤ n-ℓ,
then the parameter converges, 
i.e., 

    lim_k→+∞ϕ_1^(k) = ϕ_1,




Let 𝑧_t^(k) = 𝑦_t^(k) - 𝑥_t and 
𝑧̂_t^(k) = 𝑦̂_t^(k) - 𝑥_t.
Suppose that point t is selected to repair in the k-th iteration, 
having 

    𝑧_t^(k+1) =𝑧̂_t^(k)= ϕ_1^(k)𝑧_t-1^(k)

where t  = _i ∈ [1,n]|𝑧̂_i^(k)| 
referring to the minimum change.
Since one and only one point will be changed in each iteration, we have
𝑧_i^(k+1) = 𝑧_i^(k), i ≠ t.

For the (k + 1)-th iteration, any point i has

    𝑧̂_i^(k+1)= ϕ_1^(k+1)𝑧_i-1^(k+1).

It follows 

    |𝑧̂_t+1^(k+1)| 
       = |ϕ_1^(k+1)𝑧_t^(k+1)| 
    = |ϕ_1^(k+1)|
    |ϕ_1^(k)|
    |𝑧_t-1^(k)| 
       < |ϕ_1^(k+1)|
    |𝑧_t-1^(k)|
    = |𝑧̂_t^(k+1)|       (since |ϕ_1^(k)| < 1) 
       < |ϕ_1^(k+1)|
    |𝑧_i^(k+1)|
    = |𝑧̂_i+1^(k+1)|

where i ≠ t. 
That is, (t+1)  = _i ∈ [1,n]|𝑧̂_i^(k+1)|,
will be the minimum repair in the (k+1)-th iteration.

Similar repairing applies in m = n - ℓ iterations,
having 𝑧_n^(m) = ϕ_1^(1)…ϕ_1^(m-1)𝑧_ℓ^(0), 
𝑧_n-1^(m) = ϕ_1^(1)…ϕ_1^(m-2)𝑧_ℓ^(0), 
according to Lemma <ref>.
For the estimation of 𝑧̂_i^(m), 
we have n = _i ∈ [1,n]|𝑧̂_i^(m)| 
since |𝑧_n-1^(m)| = min_i ∈ [1,n-1]|𝑧_i^(m)|.
That is, point n will always be  repaired in the following iterations.

According to Equation <ref>, we have

    |ϕ_1^(m+1) - ϕ_1^(m)|
       = |𝑧_1^(m+1)𝑧_2^(m+1)
    +…+
    𝑧_n-1^(m+1)𝑧_n^(m+1)/
    (𝑧_1^(m+1))^2
    +…+
    (𝑧_n-1^(m+1))^2
    
            - 𝑧_1^(m)𝑧_2^(m)
    +…+
    𝑧_n-1^(m)𝑧_n^(m)/
    (𝑧_1^(m))^2
    +…+
    (𝑧_n-1^(m))^2
    | 
        = |𝑧_n-1^(m)
    (𝑧_n^(m+1)-𝑧_n^(m))
    /
    (𝑧_1^(m))^2
    +…+
    (𝑧_n-1^(m))^2
    | 
        = |
    (𝑧_n-1^(m))^2
    /
    (𝑧_1^(m))^2
    +…+
    (𝑧_n-1^(m))^2
    |·
    |ϕ_1^(m)-ϕ_1^(m-1)| 
        < |ϕ_1^(m)-ϕ_1^(m-1)|

where 𝑧_i^(m) = 𝑧_i^(m+1), i < n. 
It is proved that the parameter converges.


It is worth noting that the condition |ϕ_1^(k)| < 1 in Proposition <ref> for the parameter to converge 
could be commonly observed in real data. 
First, referring to <cit.>, most time series in practice are stationary, 
which is guaranteed to have |ϕ_1^(k)| < 1 for 𝑝=1. 
Moreover, for non-stationary cases, a typical processing way is to transform it to stationary via differencing <cit.>.



For IMR(1) with first ℓ data points labeled in 𝑦^(0), 
if the parameter converges, having 
lim_k→+∞ϕ_1^(k) = ϕ_1,
then the converged repair result is

    lim_k→+∞𝑦_t^(k) = 𝑦_t = 
    𝑦_i^(0)     if  i ∈ [1,ℓ] 
    ϕ_1^i-ℓ(𝑦_ℓ^(0) - 𝑥_ℓ) + 𝑥_i      if  i > ℓ

where the converge parameter can be directly calculated by 

    ϕ_1 = 
    (𝑦_1^(0)-𝑥_1)(𝑦_2^(0)-𝑥_2)
    +…+
    (𝑦_ℓ-1^(0)-𝑥_ℓ-1)(𝑦_ℓ^(0)-𝑥_ℓ)
    /
    (𝑦_1^(0)-𝑥_1)^2
    +…+
    (𝑦_ℓ-1^(0)-𝑥_ℓ-1)^2
    .




It is a special case of Proposition <ref> with m = 1.
Similar to the proof of Proposition <ref>, 
we have

    ϕ_1    = 𝑧_1𝑧_2
    +…+
    𝑧_n-1𝑧_n/
    (𝑧_1)^2
    +…+
    (𝑧_n-1)^2
    
       = 𝑧_1𝑧_2
    +…+
    𝑧_ℓ-1𝑧_ℓ
    + ϕ_1(𝑧_ℓ)^2
    + … + ϕ_1^2n-2ℓ-1(𝑧_ℓ)^2
    /
    (𝑧_1)^2
    +…+
    (𝑧_ℓ)^2
    + ϕ_1^2(𝑧_ℓ)^2
    +…+
    ϕ_1^2n-2ℓ-2(𝑧_ℓ)^2
    
       = 𝑧_1𝑧_2
    +…+
    𝑧_ℓ-1𝑧_ℓ/
    (𝑧_1)^2
    +…+
    (𝑧_ℓ-1)^2
     .



It is not surprising that the converged parameter ϕ_1 
in Equation <ref>
in Proposition <ref>
is exactly the solution of 
Equation <ref>  
in Proposition <ref>
with m=1
for the special case of one labeled segment.  




§ EFFICIENT PARAMETER ESTIMATION


Among the three major steps in Algorithm <ref>, 
while the repair candidate generation and evaluation (in Sections <ref> and <ref>) are inevitable for the minimum repair,
we show in this section that the costly parameter ϕ^(k) estimation (in Section <ref>) is optimizable 
in our iterative repairing scenario. 
First, we identify that the matrices 𝑍^(k),𝑉^(k) for parameter estimation could be pruned by simply removing rows with value 0, in Section <ref>. 
Moreover, incremental computation could be designed such that the complexity of parameter estimation is reduced from O(n) to O(1), in Section <ref>.





 §.§ Matrix Pruning





  
*Intuition
Recall that when estimating the parameter ϕ^(k) in Equation <ref>, 
we need to consider two large matrices 𝑍^(k) and 𝑉^(k) with sizes (n-p)× p and (n-p)× 1, respectively. 
The value 𝑦_i^(k) - 𝑥_i in 𝑍^(k) and 𝑉^(k) denotes the difference between the labeled/repaired value 𝑦_i^(k) and the input value 𝑥_i of point i.
In practice, the labeled data is often limited, while the repaired data should not be significantly changed referring to the minimum change principle of repairing. 
That is, most values in 𝑍^(k) and 𝑉^(k) equal to 0. 
We show (in Proposition <ref> below) that the sparse matrices could be pruned by removing rows with value 0. 

Let 𝑧^(k)_i = 𝑦^(k)_i - 𝑥_i for simplicity.  
We rewrite 𝑉^(k),𝑍^(k) in Equation <ref> for parameter estimation, 

    𝑉^(k) =    [ 𝑧_𝑝 + 1^(k); 𝑧_𝑝 + 2^(k);           ⋮;     𝑧_n^(k) ]
    ,   𝑍^(k) =  
    [   𝑧_𝑝^(k)         …   𝑧_1^(k); 𝑧_𝑝+1^(k)         …   𝑧_2^(k);         ⋮         ⋱         ⋮; 𝑧_n-1^(k)         … 𝑧_n-𝑝^(k) ]
    .



The following conclusion states that 
the same parameter ϕ^(k) could still be computed by Equation <ref>, after removing the rows in 𝑍^(k), whose values equal to 0, and the corresponding rows in 𝑉^(k).



For any row in 𝑍^(k), denoted by 𝑍_r^(k), 
if the entire row are all with value 0, 
i.e., 𝑧_r+𝑝-1^(k)=𝑧_r+𝑝-2^(k)=…=𝑧_r^(k)=0,
then it is safe to remove the row 𝑍_r^(k) and the corresponding row 
𝑉_r^(k) = [ 𝑧_𝑝+r^(k) ] 
from matrices 𝑍^(k) and 𝑉^(k),  respectively, 
which still compute that the same ϕ^(k).



We rewrite Equation <ref> as
ϕ^(k) = (𝐴^(k))^-1𝐵^(k), 
where 
𝐴^(k) = (𝑍^(k))'𝑍^(k) 
and 
𝐵^(k) = (𝑍^(k))'𝑉^(k)
as illustrated in Equations 
<ref>
and 
<ref>,
respectively. 
To prove the same parameter ϕ^(k), 
we show that all the values in 
𝐴^(k) 
and 
𝐵^(k) 
are unchanged after removing the r-th row with value 0.

(1) For 
a^(k)_ii, 1≤ i≤ p in Equation <ref>, 
we have 

a^(k)_ii
= ∑_l=𝑝+1-i^n-i𝑧_l^(k)𝑧_l^(k)
=∑_l≠ r+𝑝-i𝑧_l^(k)𝑧_l^(k)
+ 𝑧_r+𝑝-i^(k)𝑧_r+𝑝-i^(k)
= ∑_l≠ r+𝑝-i𝑧_l^(k)𝑧_l^(k) + 0.
That is, a^(k)_ii will not be affected without considering 
𝑧_r+𝑝-i^(k) in the r-th row  𝑍_r^(k).

(2) For 
a^(k)_ij, a^(k)_ji ,
1≤ i≤ p, 1≤ j≤ p in Equation <ref>, 
we have
a^(k)_ij = a^(k)_ji 
=∑_l=𝑝+1-i^n-i𝑧_l^(k)𝑧_l-j+i^(k)
=∑_l≠ r+𝑝-i𝑧_l^(k)𝑧_l-j+i^(k)
+ 𝑧_r+𝑝-i^(k)𝑧_r+𝑝-j^(k)
=∑_l≠ r+𝑝-i𝑧_l^(k)𝑧_l-j+i^(k) + 0.
Again, a^(k)_ij, a^(k)_ji will not be affected without considering 
𝑧_r+𝑝-i^(k), 𝑧_r+𝑝-j^(k) in the r-th row  𝑍_r^(k).

(3) For
b^(k)_i  ,
1≤ i≤ p in Equation <ref>,  
we have 
b^(k)_i 
= ∑_l=𝑝+1^n𝑧_l^(k)𝑧_l-i^(k)
= ∑_l≠ r+𝑝𝑧_l^(k)𝑧_l-i^(k)
+ 𝑧_r+𝑝^(k)𝑧_r+𝑝-i^(k)
= ∑_l≠ r+𝑝𝑧_l^(k)𝑧_l-i^(k) + 0.
Finally, b^(k)_i will not be affected without considering 
𝑧_r+𝑝^(k), 𝑧_r+𝑝-i^(k) in the r-th row  𝑍_r^(k) as well.





[Parameter estimation with matrix pruning, Example <ref> continued]
Consider again 𝑥 = {6, 10, 9.6, 8.3,  7.7,  
5.4, 5.6, 5.9, 6.3, 6.8, 7.5, 8.5} and 𝑦^(0) = {6, 5.6, 5.4, 8.3,  7.7, 5.4,  5.6, 5.9, 6.3, 6.8, 7.5, 8.5} 
in Example <ref>. 
Given order 𝑝=1, we have 𝑉^(0) = {-4.4,-4.2,0,0,0, 0,0,0,0,0,0}' 
with size 11×1,
and 𝑍^(0) = {0,-4.4,-4.2,0,0,0, 0,0,0,0,0}'
with size 11×1.
Referring to Proposition <ref>, all the rows except the fourth and fifth rows in 𝑍^(0) and the corresponding rows in 𝑉^(0) can be removed.
After matrix pruning, 
we have 𝑉^(0) = {-4.2,0}' 
with size 2×1, and 
𝑍^(0) = {-4.4,-4.2}' 
with size 2×1.
Referring to Equation <ref>, the parameter is estimated by

    ϕ_1^(0) = (-4.4)*(-4.2)/(-4.4)^2+(-4.2)^2 = 0.5.

The computed parameter is the same as in Example <ref>.




 §.§ Incremental Computation





  
*Intuition
In each iteration in Algorithm <ref>, 
the parameter ϕ^(k) is estimated by Equation <ref> 
w.r.t. 𝑍^(k) and 𝑉^(k)
over all the n points. 
However, referring to the minimum change principle in Section <ref>, 
there is only one point, say r, which is changed in each iteration, 
having 𝑦_r^(k)≠𝑦_r^(k-1). 
That is, most values in 𝑍^(k) and 𝑉^(k) are the same to 𝑍^(k-1) and 𝑉^(k-1). 
We show (in Proposition <ref> below) that 
ϕ^(k) can be incrementally computed by considering only the changed values 
instead of the entire 𝑍^(k) and 𝑉^(k).
The time complexity of parameter estimation in each iteration is thus reduced from 
linear
time to constant time. 

 



  §.§.§ Recursive Formula


To enable the incremental computation, we rewrite Equation <ref> for parameter estimation as follows, 

    ϕ^(k)   = (𝐴^(k))^-1𝐵^(k)
    𝐴^(k)   = (𝑍^(k))'𝑍^(k) = 
     [ a_11^(k) a_12^(k)        ⋯ a_1p^(k); a_21^(k) a_22^(k)        ⋯ a_2p^(k);        ⋮        ⋮        ⋱        ⋮; a_p1^(k) a_p2^(k)        ⋯ a_pp^(k) ]
    𝐵^(k)   = (𝑍^(k))'𝑉^(k)= 
     [ b_1^(k); b_2^(k);       ⋮; b_p^(k) ]

where 

    a^(k)_ii   = ∑_l=𝑝+1-i^n-i𝑧_l^(k)𝑧_l^(k) ,
           1≤ i≤ p 
    
    a^(k)_ij   = a^(k)_ji =∑_l=𝑝+1-i^n-i𝑧_l^(k)𝑧_l-j+i^(k),
           1≤ i≤ p, 1≤ j≤ p , j>i 
    
    b^(k)_i   = ∑_l=𝑝+1^n𝑧_l^(k)𝑧_l-i^(k) ,
           1≤ i≤ p


The following conclusion illustrates that 
𝐴^(k) and 𝐵^(k), 
with sizes p× p and p× 1, respectively, could be recursively computed from 
𝐴^(k-1) and 𝐵^(k-1). 
The pre-defined p in the AR(p)/ARX(p) model is a fixed value in the algorithm and 
has p ≪ n. 
The computational cost of Equation <ref> over 𝐴^(k) and 𝐵^(k) with sizes on p 
will be significantly lower than 
Equation <ref> 
w.r.t. 𝑍^(k) and 𝑉^(k)
over all the n points.



Let r be the changed point in the current repair iteration, having 
𝑦_r^(k)≠𝑦_r^(k-1) or equivalently 
𝑧_r^(k)≠𝑧_r^(k-1). 
𝐴^(k), 𝐵^(k) could be recursively computed from 
𝐴^(k-1), 𝐵^(k-1). 

That is, for 1≤ i ≤ p, we have 

    a_ii^(k)   = 
    a_ii^(k-1) + 
       
    0 
         if  r < 𝑝+1-i  or  r > n-i 
    𝑧_r^(k)𝑧_r^(k) -𝑧_r^(k-1)𝑧_r^(k-1)     if 𝑝+1-i ≤ r ≤ n-i

For 1≤ i ≤ p, 1≤ j ≤ p, i<j, we have 

    a_ij^(k)   = a_ji^(k)
    =a_ij^(k-1) + (𝑧_r^(k)-𝑧_r^(k-1))×
       
    0 
         if  r < 𝑝+1-j  or  r>n-i 
    𝑧_r+j-i^(k-1)     if  p+1-j ≤ r < p+1-i 
    𝑧_r-j+i^(k-1)     if  n-j < r ≤ n-i 
    
    (𝑧_r+j-i^(k-1)+𝑧_r-j+i^(k-1))
         if  p+1-i ≤ r ≤ n-j


For 1≤ i ≤ p, we have 

    b_i^(k)   
    =b_i^(k-1) + (𝑧_r^(k)-𝑧_r^(k-1))×
       
    0 
         if  r<𝑝+1-i 
    𝑧_r+i^(k-1)     if  p+1-i ≤ r < p+1 
    𝑧_r-i^(k-1)     if  r>n-i 
    
    (𝑧_r+i^(k-1)+𝑧_r-i^(k-1))
         if  p+1 ≤ r ≤ n-i






We show correctness of Equations <ref>-<ref> as follows.


(1) For a_ii^(k) in Equation <ref>, 
it is the summation of 𝑧_l^(k) values for data points l starting from p+1-i to n-i.
Recall that there is only one repaired point r in each iteration 
with value changed from 𝑧_r^(k-1) to 𝑧_r^(k). 
If r is in the range of summation, i.e., 
𝑝+1-i ≤ r ≤ n-i as illustrated in Figure <ref>(a), 
we need an update on a_ii^(k-1) for a_ii^(k) by adding 
𝑧_r^(k)𝑧_r^(k) and removing 𝑧_r^(k-1)𝑧_r^(k-1). 
Since other points are not modified, the correctness of Equation <ref> is illustrated. 

(2) For a_ij^(k) in Equation <ref>, 
the changed point r may affect two terms,  
𝑧_r^(k)𝑧_r-j+i^(k)
and 
𝑧_r+j-i^(k)𝑧_r^(k), 
in the summation with l=r and r+j-i, respectively. 
A term 𝑧_l^(k)𝑧_l-j+i^(k) is represented by `×' in Figure <ref>(b).
There are four cases:


i) If both terms are not in the summation range, i.e., 
r<𝑝+1-j or r>n-i, 
a_ij^(k) will not change.


ii) For 𝑝+1-i≤ r < 𝑝+1, 
the term 𝑧_r+j-i^(k)𝑧_r^(k) will be under consideration in the summation. 
We need an update on a_ij^(k-1) for a_ij^(k) by adding 
𝑧_r+j-i^(k)𝑧_r^(k) and removing 
𝑧_r+j-i^(k-1)𝑧_r^(k-1). 
Since point r+j-i is not changed, we have 
𝑧_r+j-i^(k)=𝑧_r+j-i^(k-1).
The correctness of the second case in Equation <ref> is illustrated. 


iii) Similarly, for n-j< r ≤ n-i, 
the term 𝑧_r^(k)𝑧_r-j+i^(k) is under consideration. 
An update 
(𝑧_r^(k)-𝑧_r^(k-1))×𝑧_r-j+i^(k-1)
needs to be applied on a_ij^(k-1) for a_ij^(k).


iv) If both terms are in the range, i.e., 
p+1-i ≤ r ≤ n-j, 
we introduce an update on a_ij^(k-1) for a_ij^(k) by adding 
(𝑧_r^(k)-𝑧_r^(k-1))×(𝑧_r+j-i^(k-1)+𝑧_r-j+i^(k-1)).

(3) For b_i^(k) in Equation <ref>, 
similar to a_ij^(k), 
the changed point r may affect two terms,  
𝑧_r^(k)𝑧_r-i^(k)
and 
𝑧_r+i^(k)𝑧_r^(k), 
in the summation with l=r and r+i, respectively. 
There are four cases to consider as well, as illustrated in Figure <ref>(c).
Similar proofs in (2) for each case apply.   







  §.§.§ Recursive Algorithm


Algorithm <ref> shows the procedure of incrementally estimating the parameter ϕ^(k).
For k=0 in the first iteration, computing 𝐴^(0), 𝐵^(0) by Equations <ref> and <ref> w.r.t.  matrices 𝑍^(0),𝑉^(0), however, is inevitable. 
Nevertheless, the efficient pruning of rows with value 0 in Proposition <ref> can be applied 
as presented in Line <ref> in Algorithm <ref>. 



For the following iterations k>0, the recursive computing of 𝐴^(k), 𝐵^(k) from
 𝐴^(k-1), 𝐵^(k-1) performs. 
As presented in Proposition <ref>, 
all the p^2+p values in 𝐴^(k), 𝐵^(k) can be recursively updated in constant time. 
Consequently, 
the complexity of parameter estimation is reduced from O(n) (referring to Equations <ref>-<ref>) to 
O(1) in Equations <ref>-<ref>.



[Parameter estimation using incremental computation, Example <ref> continued]

Consider again 𝑥 = {6, 10, 
9.6, 8.3, 7.7, 5.4, 5.6, 5.9, 6.3, 6.8, 7.5, 8.5} and 𝑦^(0) = {6, 5.6, 5.4, 8.3, 7.7, 5.4, 5.6, 5.9, 6.3, 6.8, 7.5, 8.5} 
in Example <ref>.
We have 𝑉^(0) = {-4.4,-4.2,0,0,0, 0,0,0,0,0,0}' 
and 𝑍^(0) = {0,-4.4,-4.2,0,0,0, 0,0,0,0,0}'.
Given 𝑝=1, the matrices 𝐴^(0), 𝐵^(0) have only one element, 
with a_11^(0)=(-4.4)^2+(-4.2)^2=37 and b_1^(0)=(-4.4)*(-4.2)=18.48
initialized by Line <ref> in Algorithm <ref>.

According to Examples <ref> and <ref>, 
the repaired point is 𝑦_4^(1)=6.2.
We have 𝑧_4^(1)=6.2-8.3=-2.1, while 𝑧_4^(0)=0.
Line <ref> in Algorithm <ref> incrementally computes 
𝐴^(1), 𝐵^(1) from 
𝐴^(0), 𝐵^(0), 
i.e.,  a_11^(1)=a_11^(0)+(-2.1)^2=41.41
by using the incremental update in Equation <ref>, and
b_1^(1)=b_1^(0)+(-2.1-0)*(-4.2+0)=27.3 referring to Equation <ref>.
Finally, the parameter ϕ_1^(1) is computed according to 
𝐴^(1), 𝐵^(1)
by using Equation <ref>, 
i.e., ϕ_1^(1)=27.3/41.41=0.66.









§ EXPERIMENT


In this section, we experimentally compare our proposed methods IMR with the state-of-the-art approaches,
including the anomaly detection methods using 
(1) AR <cit.>, (2) ARX <cit.>, (3) ARIMA <cit.>,
(4) Tsay <cit.> as models, 
(5) the smoothing-based method EWMA <cit.>, and 
(6) the constraint-based approach SCREEN <cit.>.




  
*GPS data with real errors 

In the GPS dataset,
real errors are naturally embedded and the corresponding ground truths are manually labeled. 
It collects GPS readings by a person carrying a smartphone and walking around at campus. Since we know exactly the path of walking, a number of 186 dirty points out of 742 points in trajectory are  manually identified. 
True locations of dirty points are also manually labeled, as ground truth. 
(See major results in Section <ref>.)



  
*ILD data with synthetic  errors 

The Intel Lab Data (<http://db.csail.mit.edu/labdata/labdata.html>,
ILD) 
includes a number of measurements taken from 54 sensors for every 31 seconds in about 38 days. 
Taking 31 seconds as one epoch and omitting the missing data, a dataset of 4912 points is obtained in sensor 1 from Feb 29th to Mar 1st. 
We synthetically inject errors into the data, by shifting the values for an amount of 3 with variance 0.1 under Gaussian distribution (see some examples in Figure <ref>). 
Such “shifting” errors are very common in practice, for example the sensor is stuck for a short while, or unit error in collection in a period. 
(See major results in Section <ref>.)









  
*Criteria

RMS error <cit.> is employed to evaluate the repair.
Let 𝑥^truth be the ground truth of clean sequence, 𝑥^dirty be the observation sequence with faults embedded, and 𝑥^repair be the repaired sequence.  
The RMS error <cit.> is given by: 

    Δ(𝑥^truth,𝑥^repair)=
    √(1/n∑_i=1^n(𝑥^truth_i-𝑥^repair_i)^2).

The measure evaluates the distance between the ground truth and its repair result. 
Low RMS error is preferred. 











 §.§ Experiments on Real Errors



The experiments on real errors over GPS data consider various algorithm settings, 
including 
(1) order 𝑝, 
(2) convergence threshold τ, 
(3) max-num-iterations,
and 
(4) labeling rate. 
Similar results are also observed in ILD and omitted. 






  §.§.§ Example Results




Figure <ref> illustrates an example part (in latitude, after transformed) of the GPS dataset, including the collected observations with errors, the labeled truth, and the repair results by different methods. 
Owing to various influences such as buildings,  GPS readings may deviate from the truth.  
For instance, the data between time points 400 and 500 are collected from a place near a high building, where significant errors are observed. 
As shown, the proposed IMR shows a repair closest to the truth, compared to other methods.


 



























  §.§.§ Varying Order 𝑝







Figure <ref> presents the results on varying order p, for 
AR(𝑝), ARX(𝑝) and IMR(𝑝). 
First, as shown in Figure <ref>(a), AR-based method shows better performance with the increase of order p, where more historical values take effect in the predication of each value.
An excessively large p, however, does not show further improvement, 
since this simple model may not be able to capture the complicated semantics in a large window.  
Similar results are also observed in ARX for the same reason. 
Remarkably, owing to the iterative strategy with minimum repair in each iteration, our proposed IMR method already achieves low RMS error of repairing even with p=1. 
The results verify the necessity of analyzing the special case of IMR(1) with p=1 in Section <ref>.



It is not surprising that the iterative IMR needs higher time costs in Figure <ref>(b) than other existing methods with only one pass through the data. 
In addition to the original IMR in Algorithm <ref>, 
we also present the results of IMR with matrix pruning (IMR-MP) in Section <ref> 
and incremental computation (IMR-IC) in Algorithm <ref> for efficient parameter estimation. 
IMR, IMR-MP and IMR-IC show exactly the same accuracy results in Figure <ref>(a).
Both efficient estimation methods improve the time costs in Figure <ref>(b). 
In particular, IMR-IC for incremental parameter estimation with constant time significantly reduces time costs. 

In summary, as illustrated in Figures <ref> and <ref> over the GPS and ILD datasets, respectively, our proposed IMR has no clear preference of order 𝑝 in repairing accuracy, 
while larger order 𝑝 leads to higher time cost.



  §.§.§ Varying Convergence Threshold τ



Figure <ref> reports the results by varying the threshold τ (with p=3). 
By setting a small τ, IMR needs more iterations to converge.
The corresponding time costs in Figure <ref>(b) are higher. 
Better repairing performance is achieved by IMR with a small τ, as shown in Figure <ref>(a).
However, by further reducing the threshold τ, e.g., from 0.1 to 0.01, the repair accuracy could hardly be further improved, while the corresponding iterations and time costs significantly increase. 
On the other hand, by increasing the threshold τ, the time costs reduce. 
Indeed,
the threshold τ provides a trade-off between repair accuracy and time costs for IMR. 

In summary, a lower threshold indeed leads to better results (lower RMS error) and needs more iterations (higher time costs). 
See Figure <ref> on ILD for more clear impact of the threshold.





  §.§.§ Specifying Maximum Number of Iterations



In Section <ref>, we analyze several special cases, 
where repairing is guaranteed to converge in theory under some conditions.  
For general cases where such conditions are not met, 
(although all the experiments converge in Section <ref> 
under various settings with/without theoretical convergence guarantee),  
one may specify the maximum number of iterations, as a remedy in practice to avoid waiting for convergence.
That is, Algorithm <ref> terminates when the iteration number reaches max-num-iterations, even if the convergence condition in Line <ref> is not met. 

Figures 
<ref> evaluates various settings of maximum number of iterations 
(average time costs are reported by repeating each test 10 times).
As illustrated, a moderately large number of iterations already achieve good repair results, i.e., close to the (right-most) converged results.



  §.§.§ Varying Labeling Rate


Figure <ref> illustrates the results on various labeling rates. 
A labeling rate 0.1 denotes that 10% data points are labeled with truth in the dataset. 
It is not surprising that the higher the labeling rate is,  the better the repair performance of IMR and ARX will be, which utilize the labeled truth, 
as illustrated in Figure <ref>(a). 
An interesting result is that with the increase of labeling rate, 
the corresponding time costs in Figure <ref>(b), first increase and then drop. 
The reason is that for a small labeling rate (say 0.06) with points barely modified, 
the iterative repair can quickly converge, while leaving most dirty data unchanged. 
The corresponding RMS error is high in this case as shown in Figure <ref>(a). 
With more data labeled in the input, more dirty points will be identified and repaired by the algorithm, leading to higher computation costs. 
When the labeling rate is large, such as 0.25, a great number of dirty points may be labeled. Thereby, the iterative repair could  converge quickly again. 

As anomaly detection methods, the results of Tsay and ARIMA are generally similar to those of AR and ARX. 
SCREEN and EWMA methods are not affected by  order p, threshold τ and  labeling rate. 
It is not surprising that SCREEN performs weakly, which verifies our motivation and analysis in the Introduction. 
Similarly, since EWMA does not utilize the labeled truth, its performance is weak.



















 §.§ Experiments on Synthetic Errors




The experiments on synthetic errors over ILD data focus on varying the errors.
Again, similar results are also observed in the other dataset GPS and thus omitted. 








  §.§.§ Evaluation on Various Errors





We consider Shift and Innovational errors <cit.>, as the example illustrated in Figure <ref>, where Spike errors are considered as a special case with # consecutive errors = 1. 
The corresponding accuracy results are reported in  Figure <ref>. 
In general, similar results are observed over Innovational and Shift errors. 
That is, while our proposed IMR may not deal with Spike errors (i.e., with # consecutive errors = 1 in Figure <ref>) as good as SCREEN <cit.>, 
IMR always shows significantly better results (lower RMS measure) in addressing a large number of consecutive errors, on both Shift and Innovational error patterns. 
The results demonstrate again that 
our proposal works well in repairing consecutive errors.

As illustrated in Figure <ref>(a), not only the proposed IMR but also ARX with a small threshold τ suffers from an overcorrection when the time series shifts back to non-anomalous data. 
The number attached to each method, e.g., IMR-0.05, denotes τ=0.05 for IMR.
As shown, there is a trade-off in both IMR and ARX: 
a smaller threshold τ shows better results in dealing with consecutive errors, but leads to overcorrection when the shift ends.  
Nevertheless, as presented in Figures <ref> and <ref>, 
a smaller threshold τ generally has better overall accuracy.
It is also worth noting that existing methods, EWMA smoothing and SCREEN, cannot handle well the repairing either when the time series shifts back to non-anomalous data, as the example illustrated in Figure <ref>.













  §.§.§ Evaluation on Online Computing




As long as online labeling is available (discussed in Section <ref>), 
the proposed IMR is applicable. 
Remarkably, by interpreting all the historical data as one labeled segment, 
the direct calculation of repairs without iterating in Section <ref> can be applied for efficient online computing. 

Figure <ref> presents the results of online repairing. 
Since errors are randomly introduced in the dataset, to obtain reliable results, we repeat 10 times for each test with random error introducing and report the average.
Troughput of IMR is stable and comparable to others.
IMR again shows the best repair.
Improvement of IMR compared to ARX is not as significant as in other experiments. 
It is not surprising, referring to the similar Equation <ref> for ARX and Equation <ref> for IMR(1).
The advantage of IMR is that no threshold parameter is required for IMR(1) with one labeled segment in Section <ref>, 
while ARX needs to set threshold τ.





















  §.§.§ Experiments on Temperature



To evaluate over cyclic time series, we employ another data set on temperature in years (<http://data.cma.cn>,
with cyclic patterns as illustrated in Figure <ref>). 
Similar to ILD, we inject synthetic errors to the temperature data. 
Figure <ref> presents the results over different error types.
Again, the results are generally similar to those on ILD in Figure <ref>.
That is, the proposed IMR shows significantly better results when dealing with a large number of consecutive errors. 
SCREEN shows no better results in this dataset under Spike errors (i.e., # consecutive errors = 1), since the clean data also contain a large number of Spikes as illustrated in Figure <ref> and cannot be distinguished from errors.




§ RELATED WORK


The idea of performing repair in multiple iterations has also been studied <cit.>, where the past repairs could help in recommend more accurate repairs in the future.
The continuous data cleaning approach <cit.> however is not directly applicable in our problem, 
since it employs FD constraints which is not available over time series data. 




 §.§ Anomaly Detection over Temporal Data




AR and ARX indeed have been widely used for anomaly detection in various areas such as economics and social surveys <cit.>.
We consider ARX <cit.> in this study, since this approach can utilize the labeled truth and could be adapted to cooperate with the minimum change principle in data repairing
(as shown in Section <ref>). 

Hellerstein <cit.> surveys methods 
on cleaning errors in quantitative attributes of large databases, 
where time series data are also discussed as a special scenario. 
For instance, 
Tsay <cit.> presents unified methods for detecting and handling outliers and structure changes in a univariate timeseries. Iterative procedures consist of specification-estimation-detecting-removing cycles to handle one-by-one the most significant disturbance.

ARIMA <cit.> is a  general parametric family of time series, consisting of autoregressive process and moving-average process. It can also incorporate a wide range of nonstationary series.





 §.§ Smoothing-based Cleaning



Smoothing techniques are often employed to eliminate noisy data.
For example, the simple moving average (SMA) <cit.> smooths time series data by computing the unweighted mean of the last k points.
Instead of weighting equally, 
the exponentially weighted moving average (EWMA) <cit.> assigns exponentially decreasing weights over time.
As indicated in <cit.>, 
also illustrated in Figure <ref> in Example <ref> 
and observed in Figure <ref> (EWMA) in the experiments,
the smoothing methods may seriously alter the original correct data, and thus have low repair accuracy.
In contrast, our minimum change-based approach, applying only high confidence repairs, could preserve most the original values with a better repair accuracy. 



 §.§ Constraint-based Cleaning



Constraint-based repairing is widely considered in cleaning dirty data, such that the repaired data satisfies some given constraints and the repair modification is minimized <cit.>.
The constraints are often defined over multiple attributes <cit.>, which are not available in a univariate time series considered in this study. 
To clean sequential data, 
existing study <cit.> employs a class of speed constraints declaring that the speeds of value changes should be bounded. 
The repairing is thus to modify the sequence towards the satisfaction of such speed constraints.
This constraint-based repairing falls short in two aspects: 
(1) it cannot handle a sequence of continuous errors, and 
(2) the labeled truth is not utilized. 
Rather than hard constraints, the cleaning is further extended to maximizing the likelihood w.r.t. speeds <cit.>.
Finally, the constraints are also utilized to clean the timestamps <cit.> and the qualitative event data <cit.>, which are different and not applicable to the quantitive time series data studied in this paper. 





§ CONCLUSION



In this paper, we study the problem of repairing dirty time series data, 
given the labeled truth of some data points. 
(1)
While existing anomaly detection techniques could be adapted to repairing, we argue that 
significant deviation (between observation and predication) based anomaly detection
is inconsistent with the minimum change principle in data repairing.  
Our experiments over real datasets illustrate such inconformity of applying anomaly detection in anomaly repairing. 
(2)
We thereby propose an iterative minimum repairing (IMR) algorithm. 
By creatively performing one minimum repair in each iteration of error predication, 
the algorithm bonds the beauty of capturing temporal nature in anomaly detection 
with the minimum change in data repairing. 
Again, the experiments demonstrate the superiority of our proposal.
Remarkably, in contrast to anomaly detection approaches AR and ARX, our proposed IMR is not sensitive to the setting of order 𝑝, i.e., a small 𝑝 is sufficient to achieve high repair accuracy with low time costs. 
(3)
The convergence of IMR is explicitly analyzed. 
In particular, we show that the converged repair result could be directly calculated without iterative computing in certain cases, 
which enables efficient online repairing over streaming data. 
It is worth noting that unlike the existing ARX, no threshold needs to be specified for IMR in online computing.

(4) Finally, we design efficient pruning and incremental computation, 
which reduce the complexity of parameter estimation from linear time to constant time. 
Experiments illustrate the significant improvement on time performance by pruning and incremental computation. 




§ ACKNOWLEDGMENT


This work is supported in part by
National Key Research Program of China under Grant 2016YFB1001101;
China NSFC under Grants 61572272, 61325008, 61370055, 61672313 and 61202008;
Tsinghua University Initiative Scientific Research Program.
Shaoxu Song is a corresponding author.


abbrv






    < g r a p h i c s >

