Task-oriented Time Series Imputation Evaluation via Generalized Representers

Authors: zhixian  wang
linxiao  yang
liang  sun
qingsong  wen
yi  wang
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 06:38:15

Abstract: Time series analysis is widely used in many fields such as power energy, economics, and transportation, including different tasks such as forecasting, anomaly detection, classification, etc. Missing values are widely observed in these tasks, and often leading to unpredictable negative effects on existing methods, hindering their further application. In response to this situation, existing time series imputation methods mainly focus on restoring sequences based on their data characteristics, while ignoring the performance of the restored sequences in downstream tasks. Considering different requirements of downstream tasks (e.g., forecasting), this paper proposes an efficient downstream task-oriented time series imputation evaluation approach. By combining time series imputation with neural network models used for downstream tasks, the gain of different imputation strategies on downstream tasks is estimated without retraining, and the most favorable imputation value for downstream tasks is given by combining different imputation strategies according to the estimated gain.

Main paper: 
 introduction: 
 Time series analysis plays a crucial role in many real-world applications, such as energy, finance, healthcare, and other fields [1, 2, 3] . For example, forecasting load series forms the basis for further decision-making in power dispatch in the power grid system, thereby generating a significant amount of economic benefits [4, 5, 6] . However, collecting time series data, especially high-quality ones, is challenging. Due to the instability of the external environment, sensor failures, and even ethical and legal privacy issues, missing values are prevalent in time series data [7] . For instance, in the BDG2 load series dataset [8] , widely used in building energy analysis, the ratio of complete time series data is less than 10%.To handle missing values in time series data, numerous methods have been proposed for time series imputation in the literature. Based on the features of the imputation methods, these approaches can be divided into statistical and machine learning methods, such as ARIMA and KNN [9, 10] , as well as deep learning-based methods [11, 12, 13, 14] . Both types of methods generally use reconstruction errors of missing values to guide learning and perform evaluation. Recently, some researchers have turned their attention to evaluation strategies based on downstream task performance [15] . However, in most cases, downstream tasks are classification tasks [13] , while forecasting tasks, as another important branch of time series-related tasks, have not been fully considered. The main challenge for time series forecasting is that the time series serves as both input and label (output) for the model during training, whereas in classification tasks, it only serves as input for the model.In supervised learning, training labels influence the calculation of the loss function, which in turn affects the optimization of model parameters and, ultimately, the performance of the model on the test set. [16] indicates that noise in input data (missing data can be considered a type of noise) often has a limited impact on forecasting results. In contrast, label noise can significantly affect the model and, consequently, the final test results from the beginning to the end of the time series. Therefore, when evaluating the impact of different time imputation methods on downstream forecasting tasks, it is essential to focus on assessing the quality of training labels constructed through various imputation methods.To evaluate how the quality of the imputation labels affects downstream forecasting tasks, it is important to clarify that an excellent imputation strategy does not necessarily mean that the imputed value at each time step is superior to any other method. [15] provides a benchmark for various methods in time series imputation tasks. Although SAITS [13] , as one of the latest SOTA methods, has achieved remarkable results, there are still methods that surpass SAITS in some cases. This demonstrates that time series imputation is a complex task, making it difficult to find a universal method capable of handling all situations, let alone considering the performance of downstream forecasting tasks. A more realistic scenario is that while one method may perform better overall, it may not outperform other methods locally. In time series, this means that one method may excel in some time steps, while others do better in different time steps. However, to examine the impact of each time step on forecasting, retraining the forecasting model multiple times is necessary, which is impractical due to time and computational costs. Therefore, an efficient estimation method is needed to examine the impact of each time step with different imputation methods. Additionally, since finding a universal method is difficult, it is natural to shift focus toward combining the advantages of current methods to obtain a better imputation strategy. Consequently, determining how to combine different strategies becomes a challenge.Based on the above situation, we have developed a task-oriented time series imputation evaluation strategy. Specifically, we summarize our contributions into the following points.1. We propose a strategy that evaluates each time series imputation method by estimating the impact of missing (imputed) labels at each time step on downstream tasks without requiring multiple retraining, which significantly reduces time and computational consumption. To the best of our knowledge, we are the first to consider the impact of missing values in time series as labels on downstream forecasting tasks.2. We introduce a simple and effective similarity calculation method based on the characteristics of long time series to estimate the impact of imputed values more quickly, striking a balance between performance and computational cost.3. We develop a time series imputation framework guided by maximizing the gains of downstream tasks, enabling the combination of advantages from different time series imputation strategies to achieve a better one. This results in improved performance in the downstream forecasting task.Time Series Imputation Time series imputation can be primarily classified into two categories: traditional techniques and neural network-based techniques. Traditional methods replace missing values with statistics, such as the mean value or the last observed value [17] . They include simple statistical models like ARIMA [18] , ARFIMA, SARIMA [19] , and machine learning techniques such as KNNI [20] , TIDER [21] , MICE [22] , BayOTIDE [23] . In recent years, deep learning imputation methods have demonstrated remarkable capabilities in capturing intricate temporal relationships and complex variation patterns inherent in time series data. These methods employ deep learning models like Transformers [24, 13] , generative neural networks such as VAEs [12, 25] , GANs [26, 27] , and diffusion models [28] to capture complex dynamic relationships within time series data. Although different methods exhibit various advantages, no universal method currently outperforms others in all scenarios and datasets. This observation inspires us to consider combining existing advanced methods in this work to achieve better time series imputation strategies.Sample-based Explaination Sample-based explainable methods can be divided into two categories [29] . One is based on retraining, which evaluates the importance of corresponding data by comparing the impact of removing data points on the model and even the final test results [30, 31, 32] .Among them, the introduction of shapley value by [30] naturally ensures the fairness of data attribution. The other type is based on gradient methods, which directly estimate the influence of data points without the need for retraining. This type of method can be subdivided into three main categories, which are based on representative theories, Influence Function, and training loss trajectory. [33] is a representative of the method of the first type, whose core idea is to fix the other layers and only focus on the last layer of the neural network, so that the influence of each sample can be explicitly calculated. On the other hand, the Influence Function [34] is based on the assumption of convergence and uses the Hessian matrix to estimate the influence of samples. The last type of method takes into account the entire training process of the neural network, continuously tracking the impact of samples on each parameter update [35] . In addition, [36] summarizes the gradient-based method and unifies them as generalized representers. 
problem statement: 
 Consider a multivariate time series dataset represented by {(X i , y i )} n i=1 , incorporating n samples. In this dataset, X i ∈ R D×L1 corresponds to a feature matrix containing D distinctive features over L 1 temporal intervals, whereas y i ∈ R L2 signifies the target time series, which spans L 2 temporal intervals. It is crucial to recognize that y i may include several missing entries, a common complication within real-world datasets. For example, in the context of electrical load forecasting, X i encompasses daily weather-related time series data, comprising variables such as temperature, humidity, and wind speed, while y i represents the electrical load time series of a given day, possibly containing missing entries due to issues in data collection or transmission. Addressing missing values in {y i } through imputation is a fundamental preprocessing step for machine learning tasks involving this data, underscoring the necessity to assess the effectiveness of various imputation methods. Consider {y (1) i } and {y (2) i } as two time series resulting from the imputation of {y} via two different methods. The goal is to ascertain whether the imputation performed on {y (2) i } is superior to that on {y (1) i }. Moreover, we seek to evaluate the quality of imputation at each temporal interval, determining if the imputation of the l-th interval in y (2) i is more accurate than that in y (1) i . Conventionally, the quality of imputation is quantified by measuring the discrepancy between the imputed values and the actual data, favoring methods that minimize this deviation. In this study, however, we propose to assess imputation quality based on the performance of subsequent tasks.One step further, we evaluate the quality of imputation on a timestep basis, examining if the imputation for the l-th interval in y(2) i exhibits improved efficacy over y (1) i , thereby offering a more nuanced and comprehensive evaluation of imputation methodologies.Let us define the loss function for the downstream task as L(f (X, θ), y), where f (•, θ) denotes the model used in the downstream task parametered by θ. And let {(X v i , y v i )} m i=1 constitute a test dataset that will be used to gauge model performance. We denote y i , respectively. According to our intuition, if y(2) i,l is superior to y (1) i,l , swapping y (1) i,l for y (2)i,l should result in a decrease in the test set's loss. Guided by this rationale, we define the indicator function I(i, l), which discerns whether y(2) i,l is preferable over y (1) i,l as follows:EQUATIONk )EQUATIONEQUATIONHere, y i (2,l) denotes a vector identical to yi , except at the l-th entry, which matches that of yi . Clearly, I(i, l) ≥ 0 implies that the substitution of y (1) i,l with y (2) i,l leads to a decreased test set loss, suggesting that y(2) i,l is superior. Conversely, if I(i, l) < 0, it suggests that y (1) i,l is preferable to y (2) i,l .Despite the effectiveness of the definition provided by Equation (1), computing I(i, l) for every missing value in the dataset is impractical due to the extensive model retraining required, which can be prohibitive in terms of time. To overcome this challenge, in the next section, we put forth an efficient methodology for estimating I(i, l) without retraining the model.To compute I(i, l) efficiently, we propose a retrain-free method in this subsection. As both y (1) and y (2) are imputation of y, then we assume that y(1) i,l is close to y (2)i,l , with which we approximate I(i, l) using the first order Tolyer expansion asEQUATIONEquation ( 2) provides an approximation for computing I(i, l), where ∂f (X v k ,θ) ∂y i,l measures how the training target y i,l affect the prediction of the test data, and∂L(f (X v k ,θ),y v k ) ∂f (X v k ,θ)computes how the changing of the prediction of X v k affect the final loss. Note that the symbolic expression∂f (X v k ,θ) ∂y i,l can be conceptually broken down into ∂f (X v k ,θ) ∂θ• ∂θ ∂y i,l , elucidating the role of the label y i,l in shaping the model parameters θ throughout the training process. This, in turn, has repercussions on the model's prediction when evaluated on unseen data from the test set, i.e. f (X v k , θ). By focusing on the derivative∂f (X v k ,θ) ∂y i,l, our goal is to assess the extent to which changes in label values y i,l influence the model's predictions on the test set, thereby affecting overall model efficacy.When it comes back to the estimation, dispite∂L(f (X v k ,θ),y v k ) ∂f (X v k ,θ)and yi,l -yi,l are easy to compute, estimating ∂f (X v k ,θ) ∂y i,l is difficult. The difficulty comes from two aspects. Firstly, for the complex f (•, θ), the final parameter is not only affected by the training data, some other factors, such as the structure of the network and learning rate during the learning process. Secondly, all of the n training samples affect the parameters of the model, leading to the mixture of the effect of data points on the final model. Thus isolating the effect of a single data point is difficult.To overcome these two difficulties, we propose to approximate ∂f (X v k ,θ) ∂y i,l using a white-box model, where how each training datapoint affects the final prediction is clear from the design of the model. To this end, we propose to approximate ∂f (X v k ,θ) ∂y i,l using a kernel machine, i.e. α T i,l K (X i , X v k ), where K (X i , •) is a kenerl between the training sample X i and test samples measuring the similarity between the X i and X v k , and α is a learnerable hyperparameter. It can be proven that the indicator function based on this definition satisfies many desirable properties (please see Appendix for details) to construct an axiomatic attribution. Formally, the coefficient α i,l ∈ R L2 can be computed by solving the following optimization problem:EQUATIONTo solve this problem,EQUATIONi,l -y(2) i,lcan act as a substitute of ∂f (Xj ,θ) ∂y i,l since there is no ground truth. However, the problem is still not practical to solve because we can not obtain f (X j , θ 2 ) without retraining the model. Even though it isn't necessary to traverse all i, j, and l to obtain the complete data, it still goes against our original intention to calculate I(i, l) efficiently. Furthermore, it is also difficult for us to determine how much data is sufficient to ensure the accuracy of the solution. Fortunately, with the help of Remark 1, we can bypass the process of finding enoughf (Xj ,θ1)-f (Xj ,θ2) y (1) i,l -y (2) i,land directly approximate ∂f (Xj ,θ) ∂y i,l.Remark 1. Given two infinitely differentiable functions f (x) and g(x) in a bounded domain D ∈ R n , ||f (x) -g(x)|| is always less than ϵ. For any given δ and ϵ 2 , there exists an ϵ such that, in the domain D, the measure of the region I that satisfying || ∂f (x) ∂x -∂g(x) ∂x || > δ is not greater than ϵ 2 , i.e, m(I) ≤ ϵ 2 .Remark 1 gives us the intuition that we can avoid retraining our downstream models by first-order approximation. However, some issues still need to be clarified. First, we consider neural networks used in the downstream task as infinitely differentiable functions since in practical applications, it is unlikely for computed floating-point numbers to precisely equal non-differentiable points. Second, Remark 1 limits the definition domain to a bounded region D. Time series data is usually bounded (for example, the renewable generation sequence cannot be greater than the installed capacity), making this assumption reasonable. Finally, Remark 1 can be rephrased as the better the approximation of the original function, the better the approximation of its derivative, that is, we can have a m(D)-ϵ2 m(D) probability of fitting the derivative well. Therefore, the optimization problem (3) can be transformed into an easier one. (Note that we indicate the existence of ϵ that meets the conditions in this remark, but make no restrictions on ϵ, while it is often difficult for us to make the approximation error of the original function sufficiently small. If the ϵ in the remark is infinitely close to 0, practical applications will encounter difficulties. However, in gradient descent-based neural network training, such situations often do not hinder our practical applications. Due to the space limit, the full theoretical discussion is provided in the Appendix.)EQUATIONαi,l = ∂ α′ i ∂y i,l .Now the problem is converted to solving the problem (4). Intuitively, we solve it by projecting it onto the RKHS subspace spanned by the kernels,EQUATIONwhere K l is the kernel gram matrix defined as K l,ij = K (X i , X j ) l . Considering the first-order optimality condition, α′ i,l = -1 n ∂L( fK (Xi),f (Xi,θ))∂ fK (Xi) l [36, 37] . Recalling our goal of estimating the relationship between α′ and y i,l in (5), their relationship is still unclear. This is because the objective of the optimization problem is to construct an approximation model without considering the role of label value y i,l . To clarify the effect of y i,l , we introduce an approximation by trigonometric inequality that L(f K (X i ) , f (X i , θ)) ≤ L (y i , f (X i , θ)) + L (f K (X i ) , y i ).The second term on the right side is the loss function corresponding to the downstream model, and in the case of training convergence, this should be close to a constant. Therefore, the optimization problem (6) can be rewritten asEQUATIONFormally, the solution are α′i = -1 n ∂L( fK (Xi),y i ) ∂ fK (Xi)and αi,l = -1 n ∂ 2 L( fK (Xi),y i ) ∂ fK (Xi)∂y i,l. Furthermore, since f K (•) is used to approximate f (•, θ), we use f (•, θ) to replace f K (•) to simplify the calculation. With the above preparation, the I(i, l) can be represented as follow with NTK kernel [38] ,EQUATION)In the previous section, to gauge the effects of substituting yi,l with yi,l on the downstream task, we utilized the Neural Tangent Kernel to assess the similarity between the model outputs for inputs X i and X v k . Given that the model's output length is L 2 , the computational complexity of calculating I(i, l) for all time steps in y scales as O(mL 2 P ), where P denotes the total number of parameters in the model f (•, θ), i.e., |θ|. In numerous time series applications, such as forecasting, L 2 can be substantially large (e.g., 128), rendering the evaluation process for all imputations timeconsuming. To mitigate this challenge, we propose a method to compress the size of ∂f (Xi,θ) ∂θ . Our approach is inspired by the observation that in time series forecasting, the output of f (•, θ) typically exhibits smooth variability across different l values. Therefore, we posit that the model output f (X i , θ) resides in a low-dimensional space spanned by a limited number of smooth basis functions. In mathematical terms, f (X i , θ) ≈ A † Af (X i , θ) where A ∈ R r×L2 consists of rows each representing a predefined smooth vector, A † ∈ R L2×r is the pseudo-inverse of A, and r, which is significantly smaller than L 2 , denotes the number of basis functions employed to approximate f (X i , θ). Consequently, we can approximate ∂f (Xi,θ) ∂θ as:EQUATIONThis approximation allows for the compression of the model output, thereby reducing the number of gradients that require computation. Through this simplification, the computational complexity for calculating I(i, l) decreases to O(mrP ), substantially less than the original complexity.In our experiments, we further simplify by assuming A to be a block diagonal matrix, defined as blkdiag(1 1 , 1 2 , . . . , 1 c ), where1 1 = • • • = 1 c-1are vectors of length ⌊L 2 /c⌋ with all elements equal to 1, and 1 c ∈ R L2-(c-1)⌊L2/c⌋ is a vector with all elements equal to 1.We have introduced a method for computing the indicator I(i, l), which assesses if replacing yi,l results in a reduced loss for the downstream task. Given two sets of imputation results, {yi } n i=1 and {yi } n i=1 , derived from distinct imputation techniques, we can evaluate I(i, l) across all samples and time steps, and identify that y(2) i outperforms y(1) i at time step l if I(i, l) greater than zero and vice versa if the value is lesser. In contrast to conventional evaluation strategies, our proposed method does not necessitate the availability of the ground truth y, thereby enhancing its practical utility in myriad real-world scenarios where actual values remain unobtainable. This feature renders our approach significantly more adaptable to situations where empirical truths are elusive.Given that our proposed method can evaluate the quality of two imputations at the time step level, a natural extension is to combine these two sets of imputations to derive an improved result. Specifically, we can generate a new set of imputations y ′ i for the i-th sample, where the l-th entry is yi,l if I(i, l) > 0, and yi,l otherwise. Based on the definition of I(i, l), we anticipate that a model trained using y ′ i will incur a lower loss compared to one trained with yi , thereby yielding better imputation results. It is important to note, however, that the calculation of I(i, l) is predicated on the scenario where only the lth timestep in the ith sample from {y (1) } n i=1 is substituted with y(2) i,l . This consideration omits the potential interactions among samples. Consequently, in practical implementations, we opt to substitute only the timesteps that rank within the top c% of I(i, l) values. Our experiments, detailed in 3.3.2, confirm the efficacy of our proposed task-oriented imputation ensemble method. Our proposed ensemble process is summarized in Algorithm 1.Algorithm 1: Task-oriented Imputation Emsemble.Data: Training data {(X i , y (1) i )} n and {(X i , y (2) i )} n ; Estimated gain on Validation data {I(i, l)} i=1:n,l=1:L2 ; Downstream task model f (•, θ); Replacing percentage c% 1 Calculate threshold = P 1-c ({I(i, l)|I(i, l) > 0}) 2 for i = 1 to n, l = 1 to L 2 do 3 if I(i, l) > threshold then 4 Let y ′ i,l = y (2) i,l 5 else 6 Let y ′ i,l = y (1) i,l 7 end 8 end 9 Train f (•, θ) on {(X i , y ′ i )} n Result: f (•, θ) 3 ExperimentQuestion: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes] Justification: We cover some limitations of this work in our conclusion section (the part left for future work). Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
datasets and experiment setup: 
 To validate our method, we conduct experiments on six datasets: the GEF load forecasting competition dataset with the corresponding temperature [39] , the UCI dataset (electricity load and air quality) [40] , the Traffic dataset containing road occupancy rates2 , and two transformer datasets, ETTH1 and ETTH2 [41] . Note that we use the hourly resolution version of the UCI electricity dataset from [42] in the main experiment. In our main experiment, we set the downstream task as univariate time series forecasting, with both input sequence and prediction lengths set to 24. In addition to the GEF dataset, we implement our method on the 'OT' time series in ETH1 and ETH2, the mean value of the road occupancy rate in Traffic, the temperature in the UCI air quality dataset, and the total electricity consumption of 321 users in the UCI electricity dataset. It is important to note that there are no original missing values in these datasets. To simulate the missing values situation, we randomly set masks with lengths in [2, 4, 6, 12, 24, 48, 96, 120] , and replace the original values with the average value at the corresponding positions as the baseline. For the missing rate setting, if the missing rate is too low, the difference between different imputation methods may be small, while if the missing rate is too high, the even best imputation method will also be difficult to obtain reasonable results. Based on the above considerations, we mainly consider 40% missing rates as our main experimental setup. Meanwhile, we will provide experimental results under other missing rate settings in [30%, 50%, 60%] in the Appendix.To verify the performance of our strategy in evaluating different time imputation methods, we introduce multiple advanced time imputation methods. Firstly, as mentioned in the last subsection, we use the basic mean value imputation as the baseline. Secondly, we consider several time series imputation methods based on deep neural networks. They are GPVAE [25] and USGAN [27] based on generative neural network, mRNN [43] and BRITS [11] based on RNN structure, SAITS [13] , and ImputeFormer [44] based on attention mechanism. All the implementations of the above models are with the help of the toolkit called PyPOTS [45] . In addition, we also implement the network structure of a spatiotemporal graphs-assisted method called SPIN [46] for time series imputation. All the details of the neural network setting will be described in the Appendix. In this section, we examine the estimated gains of imputation for each time step. We divide the GEF dataset into a training set and a test set, where the training set includes load data from 2011 and the test set includes load data from 2012. We replace 40% of the load data in the training set with the average value at the corresponding position as the baseline, apply the linear interpolation method to replace the corresponding baseline, and use them as training labels separately. In the training set, we have obtained a total of 8760 training samples. Therefore, we replace the labels of each sample one by one to construct new 8760 sets of samples, retrain 8760 forecasting models, test their performance on the test set, and compare the performance of the new model with the model trained on the mean value-based samples. Note that although our estimation is done at each time step of each sample, the time consumption of replacing each time step and retraining is too high. Therefore, we replace all the time steps in each sample and sum up the benefit estimates for each time step as the benefit estimate for the entire sample after replacement. We adopt the 3-layer MLP structure used in [35] and extend it to be used for time series. In addition, we have also added a simple and widely used forecasting model, DLinear [47] , as our forecasting backbone model. As shown in Figure 1 , we use two metrics to compare the performance of different methods, where Corr represents the correlation between estimated gain and actual gain, and accuracy represents the percentage of the same sign between actual gain and estimated gain. In the MLP model, the Influence Function achieves a good correlation. However, in terms of accuracy, the accuracy of the Influence Function rapidly decreases as the percentage of estimated data to total data increases. On the contrary, our method exhibits good characteristics at different percentages. In addition, the segmented acceleration method generally shows that the more segments there are, the better the performance. Although it performs poorly in small percentages, it exhibits good performance in a wide range of data. In the DLiner model, both our method and the segmented acceleration method present similar situations to those in the MLP model. However, the Influence Function exhibits poor performance, with estimation accuracy and correlation significantly lagging behind our method.In addition, Table 1 reports the total time required for different methods to run on the GTX4090, indicating that retraining requires a significant amount of time while our method achieves a good balance between performance and time consumption. It is worth mentioning that when we focus not only on the small portion of time steps with the greatest impact but also examine the impact of relatively large time steps on the forecasting results, our segmented acceleration method achieves very good performance while reducing the required time.Similar to the previous section, we use mean value as a baseline to estimate the benefits on the validation set (not test set) obtained by imputation at each time step and replace the 10% time step with the highest benefits to train the forecasting model. In addition, the missing values in the time series damage the original characteristic of the time series. Although imputation values can repair it to some extent, they may still have adverse effects on downstream tasks. Therefore, we also introduce the Influence Function as a comparison, using it to estimate the 10% points with the worst impact on the forecasting results after imputation and remove them, which is often the case in the application of Influence Function [34] .Table 2 reports the comparison of MSE for downstream forecasting tasks (here we use the DLinear as the backbone model), and we will analyze from two aspects. 2 demonstrate that combining different imputations can enhance the performance of downstream tasks on all datasets. The improvements in the GEF dataset, Electricity dataset, and Traffic dataset are significant, while the enhancement in others is relatively less noticeable. The primary reason is that, in those datasets, the performance of other imputation methods considerably lags behind the original one, making the replacement of the original label with a newly imputed label seem less impactful. However, after incorporating our estimation, there is still a slight improvement in forecasting performance. On the other hand, when considering combining two imputation methods with similar original performance, incorporating our method will bring significant gains (as seen in the GEF datasets). It's worth noting that the mean value used as a baseline outperforms other imputation methods in most cases, primarily because the training data input is also based on mean value imputation, facilitating unified and convenient comparisons. In practical applications, we can also apply other advanced imputation methods to the input data and modify the labels based on the estimated benefits.II. Comparison between Gain Estimation and Influence Function. Part I, II, and III Table 2 indicate that discarding a certain number of samples according to the Influence Function can indeed improve the performance of forecasting; however, such improvement is not universal. In the AIR dataset, discarding some data can negatively impact the performance of most methods. This may be due to the small amount of data contained in the AIR dataset, resulting in a greater adverse effect when discarding data. The operation of discarding data can only consider one imputation method, while our method can combine any two imputation methods to achieve better results. Furthermore, the results of repeated experiments show that the strategy of modifying values at specific time steps can make performance more stable, as its variance is significantly smaller than that of discarding data.In addition to the univariate input results displayed in Table 2 , we also include the results of multivariate inputs, which are common in practical applications. For instance, when predicting power loads, temperature is a crucial external variable. A large amount of research has focused on studying the relationship between load and temperature [48, 49] . In this experimental setting, unlike multivariate forecasting, temperature plays an auxiliary role in load forecasting while there is no need to forecast temperature itself. Consequently, we conduct experiments on the GEF dataset, inputting temperature as an external variable into the model to forecast loads. Table 3 presents the performance of our model under multivariate input, which is consistent with the univariate input scenario, and incorporating our method proves to be beneficial. Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes]Justification: We have provide the public accessible link to both the code and dataset we used for experiment.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
conclusion and future work: 
 In this work, we propose to evaluate the imputation values at each time step for the impact on downstream forecasting tasks. On the one hand, our method can accurately estimate the gain of each imputation value without retraining. On the other hand, our method can also combine different time series imputation strategies based on the estimation of gain to obtain better imputation for downstream tasks. To ensure the applicability of this method in practical scenarios, we also provide an accelerated calculation method. In the future, we will focus on further downstream tasks, such as optimization tasks based on prediction values, and build an end-to-end evaluation strategy. 
acknowledgement: 
 The A Why We Need Task-oriented ImputationHere we use a toy example to illustrate that, in some cases, we can not directly use the accuracy of imputation instead of downstream tasks to evaluate the imputation method. We want to point out that better imputation accuracy does not always mean better forecasting performance, and we simulate a dataset based on the GEF dataset to illustrate this viewpoint, experimenting with a predicted length of 24. Suppose that we only observed the value at the time step nk (k≥0) and nk+1 (k≥1), just for the convenience of linear interpolation. In the first case(represented by I), we set n = 4, fill the missing value with linear interpolation, and uniformly add Gaussian noise N (0.05,0.3). In the second case(represented by II), we set n = 6 and only do the linear interpolation (shown in Figure 2 ). We put two data sets into MLP and calculated the forecasting error as shown in the following Table 4 .  
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: We cover all these information in our experiment section.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] ,Justification: Every experiment is done 3 times.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).• The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] ,Justification: We cover these information in our appendix sections.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes]Justification: The societal impacts of time series forecasting are provided in the Introduction section.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes]Justification: We have properly cited all the models used for the experiment and ensured the license and terms of use are explicitly mentioned and properly respected.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.• If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA] Justification: This paper does not release new assets.. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] Justification: This paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA] Justification: This paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.• We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: In [36] , there are several axioms desirable for a fair attribution. In this section, we modify them to a version suitable for our task and demonstrate that our method satisfies such properties. Here we use I(i, l, X k ) to represent the impact of the perturbation of the l-th step of the i-th sample on sample X k .Definition 1 (Efficiency Axiom). For any model f (•, θ), and test point X v k , an indication function I(•, •, •) satisfies the efficiency axiom iff:EQUATIONEQUATIONThis is a counterpart of the efficiency axioms in Shapley values [50] , and our method is naturally satisfying.Definition 2 (Self-Explanation Axiom). An indication function I(•, •, •) satisfies the self-explanation axiom iff there exists any training point X i having no effect on itself, i.e. I(i, l, X i ) = 0, the training point should not impact any other points, i.e. I(i, l, X) = 0 for all X.Similar to the dummy axiom in the Shapley values, our method naturally satisfies.Definition 3 (Symmetric Zero Axiom). An indication function I(•, •, •) satisfies the symmetric zero axiom iff any two training points X i , X j such that if I(i, •, X i ) ̸ = 0 and I(j, •, X j ) ̸ = 0, thenEQUATIONThis situation holds if and only if K(X i , X j ) = 0, therefore, our method satisfies.Definition 4 (Symmetric Cycle Axiom). An indication function I(•, •, •) satisfies the symmetric cycle axiom iff for any set of training points X t1 , . . . X t k , with possible duplicates, and X t k+1 = X t1 , it holds that:EQUATIONThe original definition of I(•, •, •) does not satisfy the above properties, as our goal is to estimate the impact of label perturbations on the loss. However, when removing∂L(f (X k ,θ),y k ) ∂f (X k ,θ)from the definition, our method satisfies this property.•, •) satisfies the continuity axiom iff it is continuous wrt the test data point X, for any fixed training point X i : EQUATIONEQUATIONA sufficient condition for an attribution A(•) to satisfy the irreducibility axiom is forEQUATIONWhen selecting the NTK kernel, this property naturally satisfies.Remark 1. Given two infinitely differentiable functions f (x) and g(x) in a bounded domainD ∈ R n , ||f (x) -g(x)|| is always less than ϵ. For any given δ and ϵ 2 , there exists an ϵ such that, in the domain D, the measure of the region I that satisfying || ∂f (x) ∂x -∂g(x) ∂x || > δ is not greater than ϵ 2 , i.e, m(I) ≤ ϵ 2 .Correctness. Firstly, we can relax the restrictions on the function by requiring that each dimension of the function on R n be continuous. Then, we can simplify the problem into a one-dimensional case on R. Secondly, we can let h(x) = f (x) -g(x) (note that |h(x)| < ϵ), and then our problem can be transformed into proving that for any given δ and ϵ 2 , there exists an ϵ such that, in the interval [a, b], the measure of the region satisfying |h ′ (x)| > δ is less than ϵ 2 .Let the domain I represent the set of all x that satisfy |h ′ (x)| > δ. We first need to prove that I can be rewritten as the union of several disjoint intervals I i that satisfy ∀I i , ∀x ∈ I i , h(x) > δ or h(x) < -δ.Since h ′ (x) itself is a continuous function, the division is obvious. The problem is that we need to prove the number of I i is countable.Assuming I i is uncountable, we perform n bisection on [a, b] to obtain numbers of smaller intervals like [a + ( 12 ) n , a + ( 1 2 ) n+1 ]. If I i is not countable, then no matter how large n is, we can always find two points x and y in the interval that satisfy h ′ (x) > δ, h ′ (y) > δ (here we mainly consider greater cases, the smaller case is the same), then|x -y|δ ≤ |h(x) -h(y)| ≤ |h(x)| + |h(y)| ≤ 2ϵ,and this equation does not hold as ϵ → 0, which makes contradiction. Therefore, we can rewrite I as ∞ i=1 I i . To simplify, assume that I can be written in the form of a finite number of unionsI = N i=1 I i . If the remark does not hold, we have max i m(I i ) = k, then N k ≥ ϵ 2 . Obviously, |δk| ≤ 2ϵ, then k ≤ | 2ϵ δ |. Finally, we have N | 2ϵ δ | ≥ N k ≥ ϵ 2 ,which makes contradiction as ϵ → 0. If the number of I i is infinite, we can always find N that is big enough to satisfy N k ≥ ϵ 2 and the rest is the same.Discussion on Application. As shown in the analysis above, for some given δ and ϵ 2 , we only require that |ϵ| < |x -y|δ and |ϵ| < ϵ2δ 2N . Even though the real boundary is highly related to the specific scenarios and is difficult to tell in real applications, we can have a look at some widely used examples. Here, we mainly discuss two kinds of optimizers that are widely used in neural network training and they are SGD [51] and Adam [52] .For the loss function L to be set as MSE. Note that we use X and X test to represent the input data in the training set and test set, separately. y is the training label and y i,l is the value of a time step in it. We use t to represent the training epoch and there will be T epochs total. Our goal is to use g(X test , y i,l ) to approximate f (X test , θ T , y i,l ) first and then use For SGD (as well as its variants SGD with momentum), θ T = θ 0 -T t=1 ηh t (y i,l ). Therefore, the∂ 2 θ T ∂y 2 i,l= 0, which means that the gradient of the function we want to approximate is constant. In this case, our approximation will be pretty good.For Adam, on the one hand, [35] has claimed that the first-order approximation in the SGD situation remains valid, as long as a small step-size η is used in the update. On the other hand, let θ T = θ 0 -T t=1 ηv t (h t (y i,l )), where v t represent the terms in Adam. In this situation, ∂θ T ∂y i,l will be an algebraic function with only a finite number of monotonic intervals. Therefore, for any given δ and ϵ, the ϵ 2 will not be really high since the y i,l that makes || ∂f (Xtest,θ T ,y i,l ) ∂y i,l -∂g(Xtest,y i,l ) ∂y i,l|| > δ can only appear near the local extreme point, whose number is finite.Table 5 summarizes the dataset partitioning we used. Except for the GEF data, the rest are multivariate datasets. We forecast the 'OT' sequences in ETH1 and ETH2, as well as the combined electricity consumption of all users in ELECTRITY, the average occupancy rate of all roads in TRAFFIC, and the temperature series in the AIR dataset. We include two models in our experiment. The first one is a 3-layer MLP in which the input size and output size are both 24 while the hidden size is 128. In addition, we mainly apply the simple and high-performance DLiner with default setting in [47] as our forecasting model backbone. In addition, to adapt to situations where the input and output dimensions are different, we constructed an output layer at the end of the DLiner model, mapping the output of multiple variables to the output of a specified dimension. We use the torch.SGD optimizer [53] to optimize the parameters of the model, where the learning rate is set to 0.1. The maximum epochs for each training are 300, and the patience is set to 10.We have introduced a total of five imputation methods for comparison, and all experiments were based on pyPOTS [45] except for SPIN. The hyperparameters for each method are set as shown in Table 6 . In addition, referring to [13] , we set the total training epoch to 100 and the patience to 10, while other hyperparameters are the default setting. In section 2.2, we gives the estimation of I(i, l) as follows,EQUATION)However, depending on the solution to the optimization problem (7) , we may have different forms of estimation for I(i, l). Referring to [36, 35] , when we no longer only consider the downstream model parameters θ at the moment of training convergence but also consider the entire training process, we can obtain another form of solution to the optimization problem that αi,l = -t∈[T ]:i∈B (t)η (t) |B (t) | ∂ 2 L(f (t-1) (Xi,θ),yi)∂f (t-1) (Xi,θ)∂y i,l , where t, B (t) , and η (t) here represent the t th epoch, the batch size, and the learning rate, respectively. And the I(i, l) in this situation will beEQUATION)Compared to the original estimation, repeated calculations will bring a significant computational burden. However, based on our acceleration method, the time required for such multiple calculations is still controlled within a reasonable range. In application, we calculate for each parameter updates in each epoch.In section 3.3.1, we compared the performance of our method with the Influence Function that we modified. Below, we will describe how to modify the Influence Function to fit our task.For a training point (X, y), define y l,δEQUATIONEQUATIONTherefore, I(l,X test ) = -∇ θ L (f (X test , θ), y test ) ⊤ H -1 θ ∇ y l ∇ θ L(f (X, θ), y) δ.Note that we applied the Conjugate gradients mentioned in [34] to accelerate its computation and compare it with our methods.We use 1 NVIDIA GTX 4090 GPU with 24GB of memory for all our experiments.Our estimation may not be 100% accurate compared to the actual situation, so it is possible to introduce bias in the evaluation among different imputation strategies, which may further have adverse effects on downstream tasks.In our practice, we mainly examine the benefits of modifying each time step on downstream tasks. Therefore, we mainly focus on whether the gain estimation is positive or negative without providing precise values. Based on this idea, we provide methods for accelerating calculations in Section 2.3. Here, we present a comparison between the accelerated estimate and the original estimate. Note that we conduct this experiment on three datasets and they are GEF, ELECTRICITY, and a generated time series, denoted by Brown, based on the following Python code. from fbm import FBM f = FBM(n=2281, hurst=0.75, length=1, method='daviesharte')The result is shown in the Figure 3 . Note that here we also use correlation and accuracy mentioned in the main text. 0.1666(0.0007) 0.0522(0.0000) 0.1796(0.0000) 0.0972(0.0006) 0.4182(0.0022) 0.1490(0.0001) Mean+BRITS 0.1704(0.0000) 0.0522(0.0000) 0.1795(0.0000) 0.1078(0.0000) 0.4332(0.0000) 0.1507(0.0000) Mean+MRNN 0.1707(0.0000) 0.0523(0.0000) 0.1795(0.0000) -0.4333(0.0000) 0.1508(0.0000) Mean+GPVAE 0.1708(0.0000) 0.0522(0.0000) 0.1795(0.0000) 0.1069(0.0001) 0.4308(0.0004) 0.1507(0.0000) Mean+USGAN 0.1704(0.0001) 0.0522(0.0000) 0.1795(0.0000) 0.1076(0.0002) 0.4251(0.0002) 0.1506(0.0000) Mean+SPIN 0.1693(0.0013) 0.0523(0.0001) 0.1800(0.0003) 0.1047(0.0001) 0.4302(0.0010) 0.1502(0.0005) Mean+ImputeFormer 0.1666(0.0002) 0.0522(0.0000) 0.1794(0.0000) 0.0991(0.0002) 0.4203(0.0014) 0.1498(0.0000) III.With Seg-2 Gain estimation Mean+SAITS 0.1686(0.0005) 0.0522(0.0000) 0.1799(0.0001) 0.1003(0.0005) 0.4212(0.0017) 0.1491(0.0001) Mean+BRITS 0.1724(0.0000) 0.0522(0.0000) 0.1795(0.0000) 0.1105(0.0000) 0.4355(0.0000) 0.1507(0.0000) Mean+MRNN 0.1730(0.0000) 0.0523(0.0000) 0.1795(0.0000) -0.4356(0.0000) 0.1508(0.0000) Mean+GPVAE 0.1730(0.0000) 0.0522(0.0000) 0.1795(0.0000) 0.1097(0.0001) 0.4335(0.0003) 0.1507(0.0000) Mean+USGAN 0.1724(0.0001) 0.0522(0.0000) 0.1795(0.0000) 0.1098(0.0003) 0.4290(0.0001) 0.1506(0.0000) Mean+SPIN 0.1733(0.0008) 0.0523(0.0000) 0.1803(0.0004) 0.1093(0.0003) 0.4343(0.0001) 0.1503(0.0005) Mean+ImputeFormer 0.1688(0.0004) 0.0523(0.0001) 0.1795(0.0000) 0.1021(0.0002) 0.4231(0.0008) 0.1498(0.0000) Table 7 summarizes the performance of our acceleration methods Seg-4 and Seg-2 in downstream forecasting tasks. Overall, performance gradually improves as the number of segmented segments gradually increases. In addition, even if it is only divided into two segments, our method can still bring some gain with minimal additional computational burden.For large-scale data, we applied our method to a 15-minute resolution UCI electricity dataset (with approximately 100000 training points) and we adjusted our experimental setup to input 96 points and output 96 points, and here is the result. In addition to our solution, some kinds of methods, such as robust time series forecasting, to deal with missing (anomaly) values have been proposed these days. Here we combine our method with [16] , which is one of the SOTA of such kind of method to illustrate that this kind of method is not contradictory to our approach but can be combined. Note that the hyperparameters are the same as the original paper in [16] and we replace the dataset with ours. We conduct our method on multivariate forecasting tasks and give a small example of the ELEC-TRICITY dataset. Note that we apply our method to the first three users (columns) in the dataset. NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.Justification: We confirm that our research adheres to the NeurIPS Code of Ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction).• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.https://pems.dot.ca.gov 


Decision: # Summary
This paper proposes a novel approach to evaluating time series imputation methods, focusing on their impact on downstream forecasting tasks rather than just the accuracy of the imputed data. The proposed approach efficiently estimates the effect of different imputation strategies using a similarity-based method, avoiding the need for multiple model retrainings.   Additionally, the authors present a framework that combines the strengths of various imputation techniques to enhance overall task performance, demonstrating significant improvements in downstream forecasting tasks.  

The strengths of the paper include its innovative perspective, strong theoretical foundation, and practical relevance. It addresses a gap in the field by providing a more meaningful assessment of imputation quality in real-world applications.

# Concerns and Author Feedback
However, there are concerns about unclear mathematical explanations, limited scope (focusing only on forecasting), and the computational efficiency of the proposed method. Some reviewers also noted that the paper lacks comparisons with the latest imputation techniques and questioned its generalizability to other tasks and multivariate data.

In their feedback, the authors clarified several points, added comparisons with state-of-the-art methods, took steps toward demonstrating the scalability and efficiency of their approach, and provided limited results for the multivariate setting showing consistent improvements as with univariate. They also provided additional experiments to address concerns about the necessity of focusing on downstream tasks.

# Pros and Cons
- Pros:
  - Innovative approach that shifts the evaluation of imputation methods to their impact on downstream tasks.
  - Strong theoretical foundation and practical relevance.
  - Demonstrated forecast performance improvements and scalability.
- Cons:
  - Mathematical explanations and some experimental setups are unclear.
  - Limited scope, focusing only on forecasting tasks.
  - Concerns about computational efficiency and lack of comparisons with state-of-the-art imputation methods.
  - Questions about the generalizability to other time series tasks and multivariate data.

# Final ratings
After rebuttal and discussions 4/5 reviewers were on the accept side, with one full accept, all others including reject on the borderline, with various different shortcomings listed above being the reason for maintaining borderline ratings. This is a majority on the accept side with at least one feeling more strongly than borderline for acceptance. 

# Final Decision Justification:
The decision to accept the paper is based on its innovative and practically valuable approach, which outweighs the concerns raised. The authors' responses and additional experiments further strengthen the paper's contributions, making it a worthwhile addition to the field.

Reviews: review: paper_summary: paper_summary: This paper focuses on the imputation of missing values in time series. By noticing that different imputation methods might affect the downstream forecasting tasks, this paper proposes the imputation evaluation approach regarding the downstream tasks' performance. 
Then, the authors also developed some methods to improve the efficiency of the evaluation approach and did experiments against some baselines.
strength_weakness: strength_weakness: 1. The evaluation criteria defined in (1) are based on binary comparison and do not look realistic. In real-world applications, the difference in prediction errors between the two methods can be marginal and insignificant, and calculating this binary indicator makes little sense. Meanwhile, in real-world data, the percent of missing values would be very low, and it is doubtful that imputed missing values would lead to significant differences in downstream predictions, especially considering that there are various techniques to regularize models for more robust predictions.  

2. Some mathematical presentation is confusing. For instance, how is Eq.(2) derived?  $f(X_k^v, \theta)$ is irrelevant to $y$, and thus why is there the term of the derivative of  $f(X_k^v, \theta)$ w.r.t. $y$ in Eq.(2)

   According to line 148,  Eq.(3) is for approximating $\frac{\partial f(X_k^v, \theta)}{y_{i, l}}$, while on line 156, it presents the substitute of $\frac{\partial f(X_k^v, \theta)}{y_{i, l}}$ again. This is very confusing. 

3. The experiments have unjustified and confusing setups and results. For instance, the description between lines 259 and 267 is very unclear, e.g., what is the label in a time series prediction task? why specifically replace 40% of data? 
In Fig. 1, what is the percentage w.r.t. on the x-axis?
questions: questions: See above.
limitations: limitations: No potential negative societal impact is found in this work.
paper_summary: paper_summary: This paper presents a novel strategy to evaluate time series imputation methods based on their impact on downstream tasks, without requiring multiple model retrainings. The proposed method leverages a similarity calculation to estimate the effect of imputed values efficiently, balancing performance and computational cost. Furthermore, the authors introduce a framework that combines various imputation strategies to enhance downstream task performance, demonstrating significant improvements in forecasting accuracy.
strength_weakness: strength_weakness: 1. Although the paper introduces a retrain-free method to estimate the impact of missing labels, it does not thoroughly discuss the computational efficiency of this approach compared to traditional methods, particularly for large-scale datasets.

2. The motivation of this paper is not very clear: why considering the impact of missing values as labels in downstream forecasting tasks is a valuable formulation, and why is this the optimal way to combine existing advanced methods
questions: questions: 1. Is the proposed method a post-processing evaluation or an in-processing method during model training?
2. Is the method only comparing one pair of models at a time?
3. What is the approximation error of the indicator function and how does this impact the imputation performance?
limitations: limitations: Yes
paper_summary: paper_summary: The paper proposes a task-oriented time series imputation evaluation approach that assesses the impact of different imputation strategies on downstream forecasting tasks, rather than just the accuracy of the imputed values. The authors introduce a similarity-based method to efficiently estimate the impact of imputed values on downstream task performance, and develop a time series imputation framework that combines the advantages of different imputation strategies. The paper also discusses several axioms that the proposed method satisfies, such as the efficiency, self-explanation, and symmetric zero axioms.
strength_weakness: strength_weakness: 1. The authors state in the abstract that missing values can be observed in many time series analysis tasks, but the task-oriented imputation evaluation method proposed by the authors is only designed for forecasting tasks. In the problem definition, the label of a time series dataset is similarly defined as a time series, however in the classification task, the label is a discrete value. In addition, for some unsupervised learning time series tasks, such as clustering [1], missing values can affect the performance as well. In conclusion, the authors' method is too limited in its scope of application compared to other imputation methods.

[1] Time series cluster kernel for learning similarities between multivariate time series with missing data. PR 2018.

2. The authors' motivation is to provide imputation values that are more beneficial to downstream forecasting tasks, while ignoring the accuracy of the imputation values themselves. However, for the time series imputation task, performing accurate imputations is also one of the non-negligible purposes [2]. The authors ignored the requirement for imputation accuracy when designing their methodology, and did not show experimentally whether their method would reduce imputation accuracy.

[2] Deep learning for multivariate time series imputation: A survey. Arxiv 2402.04059.

3. The authors' ultimate aim is to make more accurate predictions in the presence of missing values, which is highly relevant to the study of robust time series forecasting [3,4]. However the authors don't mention them in the paper and don't compare them as baselines in their experiments. In addition the authors should discuss the advantages and disadvantages of imputation-based and missing value robust-based forecasting methods.

[3] Weakly Guided Adaptation for Robust Time Series Forecasting. VLDB24.

[4] RobustTSF: Towards Theory and Design of Robust Time Series Forecasting with Anomalies. ICLR 24.
questions: questions: Please see weaknesses.
limitations: limitations: Yes
paper_summary: paper_summary: The authors propose a strategy that evaluates the effectiveness of various imputation methods used to fill missing values ​​at different timestamps (time series). The effectiveness of each imputation method is evaluated based on the downstream task gain. Subsequently, rather than filling missing values with a single imputation method, each missing value at a specific timestamp is filled with the most accurate imputation method at that timestamp.
strength_weakness: strength_weakness: - The experimental section 3.3.1 is unclear and very difficult to follow. Authors must refer to the questions to know the reviewer's concerns;

- The proposed strategy may fail when the labels are multivariate time series (MTS). Indeed, in the case of MTS, many imputation methods/models (e.g. deep learning models) consist of filling missing values based on all observed values, i.e. within and across streams. Consequently, combining these imputation methods/models with others that fill in missing values by considering only the observed values of the corresponding univariate time series (e.g. the empirical mean) can lead to a high variance in the imputed values. This high variance will introduce noise into the model calculation scheme and may affect the gain of the downstream task;

- Using the mean as the sole baseline is not sufficient. Authors should adopt at least two baseline imputation methods to perform robust experiments;


- The limitations of the model are not discussed. The authors simply indicate the directions for future work.
questions: questions: - Can you please formally detail how you go from the indicator function $I(i,l)$ (equation 1) to its approximation with the first order Taylor expansion (equation 2)? It is quite confusing because in the equation we have the losses with respect to the learnable parameters $\theta_1$ and $\theta_2$ and in Equation 2 only the loss with respect to the $\theta$ parameter is used. Also, where is the initial term of the first order Taylor expansion? Did you just assume that it is equal to zero?

- Shouldn't NTK be equal to $\sum_{p=1}^{P}\frac{\partial f (\mathbf{X}_i, \theta)}{\partial \theta_p}\frac{\partial f (\mathbf{X}_k^v, \theta)}{\partial \theta_p}$ in Equation 8?

- Can you please clarify for what purpose linear interpolation is used (line 259-261)? Is it to fill in missing values ​​in the labels and thus obtain two labels, one filled in by the empirical mean and the other filled in by this interpolation?

- Have you really retrain 8760 forecasting models (line 263), or is this a typing error?

- Does the actual gain represent the gain obtained using a single imputation method, and does the estimated gain represent the gain that could be obtained using a combination of imputation methods? If not, can you explain these gains?

- Can you please elaborate more about the influence function?

- When evaluating the proposed solution with multivariate time series (MTS), are the labels also MTS, or are they still univariate time series?
limitations: limitations: - Unless the reviewer has missed something, it appears that the proposed method is only applicable when the labels are univariate time series. This makes it less generalizable, since in many real-world forecasting applications the objective is to forecast multivariate time series at consecutive discrete timestamps;

- The reviewer would suggest adding the computation time to obtain the estimation gain for each model in Table 2;

- Please correct the typo on line 139, Taylor expansion;

- The reviewer suggests that as a future work, the authors empirically evaluate the gain obtained by combining imputation methods/models that adopt the same or different approaches to filling missing values. For example, it would be interesting to see how the emprirical mean imputation works with GRU-D[1], which also uses the empirical mean to fill missing values. Another example could be to study the gain obtained when imputation models based on graph neural networks and message passing are used exclusively or combined with other imputation models, such as those based on recurrent neural networks.


[1] Z. Che, S. Purushotham, K. Cho, D. Sontag, Y. Liu, Recurrent neural networks for multivariate time series with missing values, Scientific reports 8 (1) (2018) 6085.
paper_summary: paper_summary: This paper studies the problem of evaluating time series imputation methods in terms of the performance on downstream tasks. It proposes a fine-grained metric and uses RKHS to efficiently estimate the metric. Experiments demonstrate that the proposed method achieves better estimation than the influence function approach.
strength_weakness: strength_weakness: W1: The proposed Section 2 & Appendix A are basically Sections 4.1 & 3 in [28], respectively. Thus, the paper should clearly attribute these contributions to [28]. 
- [28] Tsai et al. Sample based explanations via generalized representers. NeurIPS, 2023.

W2: Since this paper focuses on evaluation, it should provide a comprehensive evaluation. However, this paper mainly uses older imputation methods, and a number of latest imputation methods are missing (e.g., [SPIN,PoGeVon]).
- [SPIN] Marisca et al. Learning to reconstruct missing data from spatiotemporal graphs with sparse observations. NeurIPS, 2022.
- [PoGeVon] Wang et al. Networked time series imputation via position-aware graph enhanced variational autoencoders. KDD, 2023. 

Minor issue: Some sentences are unclear and hard to understand. For example, Lines 103-104 wrote "$X_i^{D\times L_1}$ corresponds to ... $L_1$ temporal intervals," but what they actually mean seems to be "one time range of length $L_1$." It might be better if the authors further polish their writing to improve clarity.
questions: questions: See weaknesses and limitations.
limitations: limitations: L1: The proposed method uses RKHS to approximate the function. However, it seems difficult to check in practice whether the approximation is accurate or not. 

L2: The problem setting does not seem to be well motivated. As what people usually care about is just an aggregated metric (like Table 2), it is unclear why it is worthy to evaluate in such a fine granularity (i.e., I(i,l)) using more computation. 

L3: This paper focuses on missing labels but does not consider missing input data. As the authors have dicussed in Section 1, missing input data can have a big impact on forecasting results if the time of missing data is close to the test data.

Comments: title: Kindly Remind for the Coming Ending of Discussion Period
Responses to the reviewer xbnY
Responses to the reviewer xbnY
Responses to the reviewer xbnY
Responses to the reviewer xbnY
Responses to the reviewer xbnY
Kindly Request for Reviewer's Feedback
Responses to the reviewer xbnY
Response to reviewer 56K9
Response to reviewer z2pu
Response to reviewer wGbj
Responses to the concerns of reviewer xbnY
Kindly Request for Reviewer's Feedback
Responses to the concerns of reviewer xbnY
Counterexample of Lemma 1
Responses to the concerns of reviewer xbnY(II)
Responses to the concerns of reviewer xbnY(I)
Response to author rebuttal
Responses to the concerns of reviewer xbnY
comment: Dear Reviewer aC85,

Thank you very much for your recognition of our response and encouragement of our work. We are very pleased that our response can answer most of your questions. We noticed that you mentioned updating the rating for our work. As the discussion period is coming to an end, could you please take the time to update the rating?

If you have any further questions about our work, we are delighted to continue discussing with you to alleviate your concerns, which will greatly help improve the quality of our work.
Dear reviewer XbnY

Thank you very much for your recognition of our work. We will organize the relevant discussions during this period and incorporate them into our revised version. Thank you again for taking the time to provide us with valuable suggestions.
Thanks for your quick response. 

Regarding your new experiment, it seems that FBM does have higher error than GEF and ELE, validating your hypothesis that the proposed matrix $A$ is more suitable for smooth $f$. Please discuss this limitation in the paper. By the way, are the results absolute or relative errors? If they are absolute errors, please also include relative error results into your paper as true gradients are typically small as well. 

My concerns have been addressed, so I have raised my score. Please revise the paper according to our discussion above.
Dear reviewer xbnY

Thank you for your feedback, and we would like to explain the two issues below.

>Regarding the future data issue

Sorry for the misunderstanding, the ELE dataset here we used is in hourly resolution [1](also cited as [34] in the original paper). Therefore, there are just 24 data points in a day. However, as you mentioned, compared to low-frequency datasets, our method has better application scenarios in high-frequency datasets because it is more conducive to obtaining more data. Therefore, we will add relevant discussions in the revised version.

>Regarding the Sec 2.3

We use the following Python code to generate the series you mentioned.
```python
from fbm import FBM
f = FBM(n=2281, hurst=0.75, length=1, method='daviesharte')
```
Then we use the GEF and ELE datasets to act as the smooth time series(we did not add missing value here). Note that we set the number of training data to 2186(about 3 months), the validation data to 96(3 days), and the predicted length to 24. Then we consider the calculation of the original approach as the true label and calculate the loss between it and the calculation of our approach using the acceleration method in Sec 2.3. The following table shows the error between the accelerated results and the original results under different numbers of segments. The brackets indicate how much the error under this segment increases compared to the error of Seg12.

|  MSE  |         GEF        |         ELE        |       Brownian      |
|:-----:|:------------------:|:------------------:|:-------------------:|
| Seg12 |      7.0735e-9     |      1.3408e-9     |      3.74636e-6     |
|  Seg8 | 7.0737e-9($\uparrow$0.0028%) | 1.3410e-9($\uparrow$0.0014%) | 3.74639e-6($\uparrow$0.0008%) |
|  Seg6 | 7.0739e-9($\uparrow$0.0056%) | 1.3411e-9($\uparrow$0.0022%) | 3.74640e-6($\uparrow$0.0010%) |
|  Seg4 | 7.0741e-9($\uparrow$0.0084%) | 1.3412e-9($\uparrow$0.0029%) | 3.74644e-6($\uparrow$0.0021%) |

It can be seen that as the number of segments increases, the error also increases, but the magnitude is very small, including the generated Brownian motion data. This shows that our acceleration method can speed up the calculation of our method with almost no loss of accuracy.



[1] G. Lai, W.-C. Chang, Y. Yang, and H. Liu, “Modeling long-and short-term temporal patterns with deep neural networks,” in The 41st international ACM SIGIR conference on research & development in information retrieval, 2018, pp. 95–104.
Thanks for your detailed elaboration! I have some follow-up questions. 

Regarding the future data issue: Thanks for your new experiment. To my understanding, the ELE dataset has a relatively high frequency (15 min). Thus, a day gives 96 new observations, so it should be fine to use a day as the validation set. However, it might be time-consuming to obtain so many new observations for low-frequency time series. Hence, your method seems to be more suitable for high-frequency data than for low-frequency data. Could you discuss this limitation in the paper?

Regarding Sec 2.3: Thanks for your explanation. If there is no theoretical guarantee, could you add an ablation study to empirically compare the error of this approximation between non-smooth and smooth time series? For non-smooth time series, you may use the fractional Brownian motion with Hurst parameter > 1/2 (say, 0.75).
Dear reviewer xbnY

Thank you very much for your quick feedback and we will address your concerns below.

>Regarding the "proof"

Thank you for your suggestion. We will continue to adjust the specific text, with the main purpose of explaining the motivation and rationality of our method rather than proof.

>Regarding future data issues.

We acknowledge that the quantity of future data does have an impact on the evaluation. Theoretically, more future data is more likely to provide better evaluation criteria. 

However, this does not mean that our method is difficult to implement when the amount of future data is small. We add a simple experiment based on the ELE dataset. We first adjusted the task to predict the value of the next 4 hours. Then we set the training set to **3 months**, the validation set for evaluation to **1 day**, and then set the test set to the next **3 months**. The rest are the same as the setting in the original paper. The specific experimental results are shown in the table below.
|  MSE  |  Mean |       SAITS      |     BRITS    |     GPVAE    |     USGAN    | ImputeFormer |
|:-----:|:-----:|:----------------:|:------------:|:------------:|:------------:|:------------:|
|   -   | 0.853 |   0.804(0.005)   | 0.854(0.001) | 0.865(0.006) | 0.894(0.017) | 0.839(0.017) |
| +Ours |   -   | **0.784(0.000)** | 0.796(0.000) | 0.797(0.004) | 0.807(0.003) | 0.793(0.005) |

Even in such an extreme environment (the validation set that can be used for evaluation is only **1 day**, and the test set is used to show whether the evaluation is reasonable), our method can still help the forecasting model achieve better results, that is, obtain imputation values ​​that are more conducive to downstream forecasting tasks.





>Regarding the section 2.3

Yes, we concur that $\boldsymbol{A}^{\dagger}\boldsymbol{A}$ is not close to the identity matrix. Indeed, $\boldsymbol{A}^{\dagger}\boldsymbol{A}$ acts as a projection matrix, where $\boldsymbol{A}^{\dagger}\boldsymbol{A}f(x)$ projects the values of $f(x)$ onto the row space of $\boldsymbol{A}$. The fitting error of $\boldsymbol{A}^{\dagger}\boldsymbol{A}f(x)$ relative to $f(x)$ in the $\ell_2$ norm equals the value of $\min_{\boldsymbol{\alpha}}\|f(x)-\boldsymbol{A}^T\boldsymbol{\alpha}\|_2$. To see why this is so, without loss of generality, let us assume that $\boldsymbol{A}$ is a full row-rank matrix. The $\boldsymbol{\alpha}$ that minimizes $\|f(x)-\boldsymbol{A}^T\boldsymbol{\alpha}\|_2$, denoted as $\boldsymbol{\alpha}^{\star}$, is given by $(\boldsymbol{A}\boldsymbol{A}^T)^{-1}\boldsymbol{A}f$. This implies that $\boldsymbol{A}^T\boldsymbol{\alpha}=\boldsymbol{A}^T(\boldsymbol{A}\boldsymbol{A}^T)^{-1}\boldsymbol{A}f=\boldsymbol{A}^{\dagger}\boldsymbol{A}f$. Therefore, we have $\|f(x)-\boldsymbol{A}^T\boldsymbol{\alpha}^{\star}\|_2=\|f(x)-\boldsymbol{A}^{\dagger}\boldsymbol{A}f\|_2$. Based on this insight, we understand that the fitting error of $\boldsymbol{A}^{\dagger}\boldsymbol{A}f(x)$ depends on the design of $\boldsymbol{A}$. 

Given that $f(x)$ represents a time series, we can leverage properties of time series to design $\boldsymbol{A}$, such as sparsity in the frequency domain or smoothness. Here, we primarily consider $f(x)$ as a smooth function that can be approximated by a piece-wise constant function. Thus, designing $\boldsymbol{A}$ in a form as presented in the paper, the fitting error equates to using an optimal piece-wise constant function to approximate $f(x)$. Since $f(x)$ is typically quite smooth, the fitting error tends to be relatively small.
Thanks for your update. My concern about Lemma 1 has been addressed. One minor issue is that the current Appendix B still "looks like" a proof. To avoid misleading the reader, please remove "contradiction" (because there is no contradiction in the counterexample) and "countable" (because you only discussed the case where $N$ is finite) and focus on the construction of the intervals instead. 

I have raised my score from 3 to 4. To help me decide whether to further increase my score, let's continue to discuss about my remaining concerns.

Regarding the future data issue: In your examples for L2, you mentioned that one can observe some future values and use them for evaluation. My concern is that near-future observations might not provide much information for evaluation because most forecasting methods should be very good at predicting near-future values. Consequently, one may need to acquire distant-future values (i.e., a large number of future values) to obtain sufficient information for evaluation. 

Regarding Sec 2.3: You claimed that $A^\dagger Af\approx f$. However, this is not necessarily a good approximation because the pseudo-inverse in general only ensures that $A(A^\dagger A)=A$, that $(A^\dagger A)A^\dagger=A^\dagger$, and that $A^\dagger A$ is Hermitian but does not ensure that $A^\dagger A$ is close to the identity matrix. The difference between $A^\dagger A$ and the identity matrix depends on the specific matrix $A$. For the matrix $A$ given in the paper, do you have any theoretical guarantee about this approximation?
Dear reviewer xbnY

Please note that the above is only a draft that we have revised. If you have any comments, please let us know.
Dear reviewer xbnY

Thank you very much for your reply. We are more than happy to provide you with our ideas for modifying the text and hope that you can provide some valuable suggestions, which will greatly help us continue to improve the quality of our work.

Overall, based on the reasons you mentioned, we will define Lemma1 as a remark, which is to demonstrate the intuition of our work. The specific description will be highlighted below.

>Remark 1. Given two infinitely differentiable functions $f(\boldsymbol{x})$ and $g(\boldsymbol{x})$ in a bounded domain $D \in R^n$. Let $\epsilon = sup_\boldsymbol{x}||f(\boldsymbol{x}) - g(\boldsymbol{x})||$. For any given $\delta$ and $\epsilon_2$, there exists  $\epsilon$ (not necessarily have to be 0 in some cases) such that, in the domain $D$, the measure of the region $I$ that satisfying $||\frac{\partial f(\boldsymbol{x})}{\partial \boldsymbol{x}} - \frac{\partial g(\boldsymbol{x})}{\partial \boldsymbol{x}}|| > \delta$ is not greater than $\epsilon_2$, i.e, $m(I)\leq \epsilon_2$. 

We will keep the derivation process (not proof) in the appendix and add some comments below. 

>As shown in the derivation process above, for some given $\delta$ and $\epsilon_2$, we only require that $|\epsilon|<|x-y| \delta$ and $|\epsilon|<\left|\frac{\epsilon_2 \delta}{2 N}\right|$. Even though the real boundary is highly related to the specific scenarios and is difficult to tell in real applications, we can have a look at some widely used examples.

Then we will give some analysis of widely used cases for MSE loss functions with SGD and ADAM optimizers. Before, we want to analyze here why the counterexample seems to look like a contradiction to our remark. 

The fundamental reason for this is that in this example, (we consider f as 0), we need to consider a function h (x) whose derivative is a periodic function and **the number of monotonic intervals is infinite**, By adjusting the parameters, we can make the period length infinitely small. In this case, even if the derivative is large (a significant difference from the derivative of f), the function can still achieve a situation where it first increases and then decreases, ensuring that the values within the domain are small (i.e., the difference from f is small).

> For the loss function $\mathcal{L}$ to be set as MSE. Note that we use $X$ and $X_{test}$ to represent the input data in the training set and test set, separately. $y$ is the training label and $y_{i,l}$ is the value of a time step in it. We use t to represent the training epoch and there will be T epochs total. Our goal is to use $g(X_{test},y_{i,l})$ to approximate $f(X_{test},\theta_T,y_{i,l})$ first and then use $\frac{\partial g(X_{test},y_{i,l})}{\partial y_{i,l}}$ to approximate $\frac{\partial f(X_{test},\theta_T,y_{i,l})}{\partial y_{i,l}}$.

>During the training process of a forecasting model, we let $h_t(y_{i,l})$ represent $\frac{\partial \mathcal{L}(f(X,\theta_t,y_{i,l}), y)}{\partial \theta_t}$ and $\frac{\partial^2 h_t(y_{i,l})}{\partial y_{i,l}^2}$ will be zero. Recalling that our goal is to approximate the gradient $\frac{\partial f(X_{test},\theta_T,y_{i,l})}{\partial y_{i,l}}$ by approximating  $f(X_{test},\theta_T,y_{i,l})$ first and then take the gradient. However, the function that we really want to approximate is $\frac{\partial \theta_T}{\partial y_{i,l}}$ since $\frac{\partial f(X_{test},\theta_T,y_{i,l})}{\partial y_{i,l}} = \frac{\partial f(X_{test},\theta_T,y_{i,l})}{\partial \theta_T}\frac{\partial \theta_T}{\partial y_{i,l}}$ and $\frac{\partial f(X_{test},\theta_T,y_{i,l})}{\partial\theta_T}$ is a constant for given $X_{test}$.

>For SGD (as well as its variants SGD with momentum), $\theta_T=\theta_0-\sum_{t=1}^T \eta h_t\left(y_{i, l}\right)$. Therefore, the $\frac{\partial^2 \theta_T}{\partial y_{i,l}^2} = 0$, which means that the gradient of the function we want to approximate is constant. In this case, our approximation will be pretty good. 

>For ADAM, on the one hand, [1] has claimed that the first-order approximation in the SGD situation remains valid, as long as a small step-size $\eta$ is used in the update. On the other hand, let $\theta_T=\theta_0-\sum_{t=1}^T \eta v_t(h_t\left(y_{i, l}\right))$, where $v_t$ represent the terms in ADAM. In this situation, $\frac{\partial \theta_T}{\partial y_{i,l}}$ will be an algebraic function with only **a finite number of monotonic intervals**. Therefore, for any given $\delta$ and $\epsilon$, the $\epsilon_2$ will not be really high since the $y_{i,l}$ that makes $||\frac{\partial f(X_{test}, \theta_T, y_{i,l})}{\partial y_{i,l}} - \frac{\partial g(X_{test}, y_{i,l})}{\partial y_{i,l}}|| > \delta$ can only appear near the local extreme point, whose number is finite.

[1] Pruthi G, Liu F, Kale S, et al. Estimating training data influence by tracing gradient descent[J]. Advances in Neural Information Processing Systems, 2020, 33: 19920-19930.
Sorry for the late response. Thanks for the authors' answers.

Most of my questions are addressed. I will update the rating accordingly.
Dear Reviewer aC85,

Since the End of author/reviewer discussions is coming in ONE day, may we know if our response addresses your main concerns? If so, we kindly ask for your reconsideration of the score. Should you have any further advice on the paper and/or our rebuttal, please let us know and we will be more than happy to engage in more discussion and paper improvements.

Thank you so much for devoting time to improving our paper!
Thanks for your understanding. Could you let me know how you plan to modify the wording? I need to make sure that it is a rigorous mathematical statement, or otherwise I cannot raise my score.
Dear reviewer xbnY

Based on your suggestion, we will make some modifications to the relevant wording. Thank you again for taking the time to provide valuable suggestions.
Dear reviewer 56K9

Thank you very much for your feedback. We are also pleased that our response has addressed your concerns. On behalf of all authors, I sincerely thank you for your recognition of our work and valuable suggestions.
Dear reviewer z2pu

Thank you very much for your feedback. Here we would like to make a simple clarification just in case of any missing.

>Regarding the baseline

We have added the comparison with a robust forecasting baseline called RobustTSF. The results show that the relationship between our approach and them is not competitive, but cooperative. That is, we can help them further obtain better forecasting results.

| MSE | Mean | SAITS | BRITS | USGAN | GPVAE | ImputerFormer | SPIN |
| :--- | :---: | :---: | :---: | :---: | :---: | :---: | :---: |
| +RobustTSF | 0.056 | 0.050 | 0.092 | 0.084 | 0.099 | 0.053 | 0.076 |
| +RobustTSF+ours | - | **0.046** | 0.050 | 0.052 | **0.046** | 0.048 | 0.051 |

>Regarding the application

Our main contribution is to propose utilizing forecasting performance to adjust imputation methods. We believe that time series forecasting is an important task, therefore it is necessary to propose specialized methods for it. As we mentioned in the introduction, it is difficult to find a perfect method that is applicable to all scenarios[1]. Combining different excellent methods for different scenarios is a possible strategy to solve this situation. Therefore, thank you very much for your suggestion. In future work, we will also consider designing corresponding strategies for different types of tasks.

Thank you again for the kind feedback.

[1] Wang J, Du W, Cao W, et al. Deep learning for multivariate time series imputation: A survey[J]. arXiv preprint arXiv:2402.04059, 2024.
Dear reviewer wGbj

Thank you very much for your feedback, and we are also very pleased that our response can address most of your concerns. We want to make a small addition regarding the additional time burden just in case of misunderstanding. We list the time for applying our method to obtain the result in Table 2 of the original paper, as shown in the table below.

As the table below shows, our method only requires about 30 seconds of computation time on the ELE dataset. For imputation methods like SAITS (which performs well in our comparison methods), the time required for imputation will exceed 1500 seconds. This means that our additional computational burden will be less than **2%**.

|         |   GEF  |  ETTH1 |  ETTH2 |   ELE  | Traffic | AIR    |
|:-------:|:------:|:------:|:------:|:------:|:-------:|--------|
| Time(s) | 19.908 | 20.233 | 20.908 | 30.170 |  29.132 | 14.701 |

Thank you again for your kind feedback.
Dear authors, thank you for your detailed answers. My concerns were clarified. I keep the score unchanged.
$\epsilon=0$ means that $f=g$ everywhere. If you allow $\epsilon=0$, then Lemma 1 holds trivially with $f=g$, and there would be nothing that needs you to prove.

I agree that your intuition helps you design your method, but intuition is not a mathematical proof. If you just want to emphasize your intuition, then do not state it as a lemma. Otherwise, if you decide to state it as a lemma, I will have to ensure it to be a rigorous mathematical statement.
>Regarding the correctness of Lemma 1

We apologize for any misunderstandings caused. Perhaps we should re-examine Lemma 1 from two perspectives: mathematical proof and practical application.

**From a mathematical proof perspective:**  We confirm that Lemma1 is mathematically correct because we are proving the existence of $\epsilon$ in Lemma1. In fact, in the example you provided, $\epsilon = 0$ can satisfy all the values of $\epsilon_2$ and $\delta$. The existence of $\epsilon$ that meets the conditions has not been challenged. Therefore, Lemma1 does not have any mathematical errors.

**From a practical application perspective:** In fact, we have not claimed that we can use Lemma1 for specific quantitative analysis, and we also acknowledge that $\epsilon$ is difficult to equal to 0 in practical applications. Lemma 1's role here is to provide us with the rationality to simplify the problem (note that rationality does not necessarily mean 100% accuracy, analogical linear approximation can introduce significant errors in some cases, but it is still widely used for simplifying problems). The proof process of Lemma1 tells us that $\epsilon$ does not necessarily need to be constant at 0. Therefore, in practical situations, we can consider using this to simplify the problem.

**We apologize again for making you feel that we are ignoring the issue and we sincerely appreciate the efforts you paid to help us.** However, the issue actually did not challenge the mathematical correctness of Lemma1. If so, you should provide a counterexample that Lemma1 does not hold when $\epsilon = 0$ (at which point $\epsilon$ that meets the condition does not exist). On the contrary, **your issue is challenging us for whether Lemma1 can help us in the real application(for example, in the counterexample you raised)**. Therefore, we have tried to explain to you our rationality in practical application.

>Detailed analysis in the application and modification of Lemma1. 

Here we would like to explain to you in more detail why your example is not seen in most of the cases. In fact, the conclusion cannot be obtained trivially, because it is related to the loss function and optimizer we choose during training. We need to analyze the specific situation. Here we take the widely used MSE loss function and two optimizers SGD and ADAM as examples and we set learning rate as $\eta$.

Remind that our goal is to calculate $\frac{\partial f\left(X_k^v, \theta(y_{i,l})\right)}{\partial y_{i, l}}$, where $\theta(y_{i,l})$ is the model parameters and related to $y_{i,l}$. Let $g_t(y_{i,l}) = \frac{ \partial \mathcal{L_t} (f_\theta(X),y)}{\partial \theta}$. Note that $\frac{\partial^2 g_t(y_{i,l})}{\partial y_{i,l}^2} = 0$ for MSE loss function. For SGD optimizer, $\theta_{t} = \theta_{t-1}-\eta g_t(y_{i,l})$. Therefore, $\theta_T(y_{i,l}) = \theta_0 -\sum_{t=1}^{T} \eta g_t(y_{i,l})$, where T is the training epoches. In this case, $\frac{\partial^2 \theta_T(y_{i,l})}{\partial y_{i,l}^2} = -\sum_{t=1}^{T} \eta\frac{\partial^2 g_t(y_{i,l})}{\partial y_{i,l}^2} = 0$ and we only need to consider the linear case. While with the ADAM optimizer, things will be a little bit complicated. For simplity, let $\theta_T(y_{i,l}) = \theta_0 -\sum_{t=1}^T\eta v(g_t(y_{i,l}))$, where v() represents the terms in ADAM. Then $\frac{\partial^{k+2} \theta(y_{i,l})}{\partial y_{i,l}^{k+2}} = -\sum_{t=1}^T\eta \frac{\partial^k(\frac{\partial^2 v}{\partial g^2}(\frac{\partial g}{\partial y_{i,l}})^2+\frac{\partial v}{\partial g}\frac{\partial^2 g}{\partial y_{i,l}^2})}{\partial y_{i,l}^k}$($k \geq 1$), which is $O(y_{i,l}^{\frac{k+1}{2k+1}})$ and tends to 0 as $k \rightarrow \infty$. The counterexample you provide is not this kind of case, that is, the derivatives of g(x) is not in the form of $O(x^{\frac{k+1}{2k+1}})$.

In the subsequent changes, we will consider incorporating these common situations into Lemma1's constraints and analyzing the changes in errors under different conditions. Thank you again for your help.
Thanks for the responses. Since some concerns remain unsolved, e.g., applications and baselines in experiments, I will keep the score.
I am just saying that your Lemma 1 is mathematically wrong (NOT saying that it is wrong in real-world scenarios because real-world scenarios have way more underlying assumptions than your Lemma 1). Basically, if you require that the second derivative is zero and that the function is infinitely differentiable, then the function is necessarily a quadratic (or linear or constant) function. In this case, Lemma 1 becomes trivial.

As your proof of Lemma 1 does not use any assumption other than those in the statement of your Lemma 1, my counterexample still applies to your Lemma 1. If you still don't see the issue, please give your "relatively large" $\epsilon_2$ and $\delta$ where you believe your Lemma 1 holds, and I will construct the corresponding counterexample.
Dear Reviewer,

Since the End of author/reviewer discussions is coming soon, may we know if our response addresses your main concerns? If so, we kindly ask for your reconsideration of the score. Should you have any further advice on the paper and/or our rebuttal, please let us know and we will be more than happy to engage in more discussion and paper improvements.

Thank you so much for devoting time to improving our paper!
Thank you again for your quick feedback.

>Regarding the counterexample

Thank you very much for providing the counterexample. However, we are sorry to say that the example you provided does not seem to meet most scenarios in real application and we will explain why in the following.

 In section 2.2 of the original paper, we state that our goal is to calculate $\frac{\partial f\left(X_k^v, \theta\right)}{\partial y_{i, l}}$, which is a symbolic expression and can be conceptually broken down into $\frac{\partial f\left(X_k^v, \theta\right)}{\partial \theta} \cdot \frac{\partial \theta}{\partial y_{i, l}}$, elucidating the role of the label $y_{i, l}$ in shaping the model parameters $\theta$ throughout the training process. This, in turn, has repercussions on the model's prediction when evaluated on unseen data from the test set, i.e. $f\left(X_k^v, \theta\right)$. Therefore, our goal is actually to estimate the impact of label $y_{i,l}$'s changes on the model parameter $\theta$. As we mentioned in the **L1 of rebuttal**, for widely used loss functions $\mathcal{L}$ like MSE,  $\frac{\partial^2 \frac{\partial \mathcal{L}\left(f_\theta(\boldsymbol{X_i}), \boldsymbol{y_i}\right)}{\partial \theta}}{\partial y_{i, l}^2}=0$, which indicates that the second derivative of the $y_{i,l}$ should make no effect on the final result. **In the example the reviewer provided, the second derivative of $x$ is not equal to zero, which does not meet most scenarios.**

During the above process, Lemma1 plays the role that allows us to change our goal from directly approximate $\frac{\partial f\left(X_k^v, \theta\right)}{\partial y_{i, l}}$ to approximate $f\left(X_k^v, \theta\right)$ first and then calculate derivatives rather than quantitative analysis. However, in Lemma1, we did not provide corresponding restrictions to obtain a generalized situation, which comes at the cost that in some cases, Lemma1 may require a relatively large $\epsilon_2$ and $\delta$, making our restrictions relatively relaxed. We apologize again for the misunderstanding caused.
Regarding Lemma 1: You cannot let $\delta$ and $\epsilon_2$ be arbitrarily large because $\delta$ and $\epsilon_2$ are **given** (according to your statement of Lemma 1). If $\delta$ and $\epsilon_2$ are allowed to be large, then your Lemma 1 becomes trivial and insufficient to support the design of your method. 

If you still don't see the issue, please check your Lemma 1 with the following counterexample. Let the domain be $D=[0,1]$, let the two functions be $f(x)=0$ and $g(x)=\epsilon\sin\big(\frac{\pi x}{\min\\{\epsilon,1\\}}\big)$, and consider $\delta=0.5$, $\epsilon_2=0.3$. It is easy to see that $f$ and $g$ are infinitely differentiable w.r.t. $x$, and that $\sup_x|f(x)-g(x)|=\epsilon$. However, no matter what $\epsilon$ is, the measure of the region where $|f'(x)-g'(x)|>\delta$ is at least $\frac23\min\\{\epsilon,1\\}\big\lfloor\frac1{\min\\{\epsilon,1\\}}\big\rfloor\ge\frac13>\epsilon_2$. (A similar counterexample can be constructed for other given $\delta$ and $\epsilon_2$ as well.)
>Regarding L2

We do not agree that evaluating **each imputed value** is pointless. We can only consider an imputation method to be completely superior to the comparison method if it outperforms the comparison method at every time step. However, in reality, such a situation is almost impossible to exist. From this point of view, we will never find the real best imputation. 

1. The importance of the fine-grained information

Let's take the first example mentioned as an example and further extreme it. Suppose that imputation method **A** only fills 99.97MW as 100.01MW while imputation method **B** only fills 5MW as 5.05MW. The rest imputation values are the same, so from the perspective of aggregated metrics, **A** will be better than **B** because the error of 0.05MW is greater than 0.04MW. However, **A** will bring much greater economic losses than **B**. In this situation, we clearly can not conclude that method **A** is superior to method **B**. 

As for the example of measuring blood pressure in patients, it is even more convincing. If evaluating at every point is not important, why do doctors still need to design round-the-clock checks? **Obviously, their goal is not to calculate the patient's average blood pressure throughout the day.** Now since **knowing every time step in the measurement is important**, then there is no reason to claim that evaluating each imputed value is pointless.

2. The future data issues and why we need forecasting tasks as supervision

When evaluating imputation, a portion of the future values can be obtained through observation over a while in the future. What truly cannot be obtained is the true label of the missing value, as it occurred in the unchangeable past. Therefore, the existing evaluations can only be used on simulated data, using aggregated metrics, which is truly deviating from real applications. Our approach of using real labels observed in the future (not simulated) and forecasting tasks as supervision is the evaluation method that can be used in real missing scenarios. Furthermore, our method can also combine existing advanced methods to achieve better effects than a single method.
Thank you again for the quick feedback.

>Regarding the Lemma 1

Sorry for the misunderstanding again. What we describe in Lemma1 is that ''$\||f(\boldsymbol{x})-g(\boldsymbol{x})\||$ is always less than $\epsilon$''. For continuous functions defined in bounded region, $\epsilon$ here is equal to $\sup_x\\|f(\boldsymbol{x})-g(\boldsymbol{x})\||$ and must exist (can not $\rightarrow$ $\infty$). This is **not the assumption** but the **definition**. **We did not make any assumption to let $\epsilon \geq \sup_x\\|f(\boldsymbol{x})-g(\boldsymbol{x})\||$ but we use $\epsilon$ to represent $\sup_x\\|f(\boldsymbol{x})-g(\boldsymbol{x})||$**. Therefore, we do not need to prove that $\epsilon \geq \sup _x\||f(\boldsymbol{x})-g(\boldsymbol{x})\||$ but rather examine the range of $\epsilon$ when Lemma1 holds for any given $\delta$ and $\epsilon_2$.

Obviously, when $\epsilon$ can be arbitrarily small, the conclusion naturally holds. However, even though $\epsilon$ is not arbitrarily small, as we mentioned above, we only require that $\epsilon<|x-y| \delta$ and $|\epsilon|<\left|\frac{\epsilon_2 \delta}{2 N}\right|$. As long as the $\delta$ and $\epsilon_2$ we set are sufficiently large (even $\rightarrow$ $\infty$), these two inequalities are bound to hold, and the cost we pay is that our fitting accuracy will decrease. For this case, we have already demonstrated the feasibility of our method through a large number of experiments, because compared to [28], our task does not require particularly high fitting accuracy.

>Regarding novelty

We don't agree with the point that our work can be simply categorized as the combination of the downstream task evaluation + interpretation. In fact, we only used a similar RKHS approximation technique to [28] (RKHS approximation is not equal to interpretation since this technique was first proposed by [29]. However, this does not mean that [28] lacks technological innovation, as they used this technique to solve new problems and properly handle issues that would arise from trivial applications). 

In our case, it is **impossible** to obtain gradients by applying the RKHS approximation technique in a trivial manner, which has been carefully analyzed in our paper. This is also why we need to change the order of gradient calculation and approximation and analyze the rationality of doing so using Lemma 1. In addition, due to the need to consider the length of multi-step outputs in our method, the computational complexity will be high. For such cases, we also emphasize that we have provided the method in section 2.3 to **alleviate the computational burden** to make the approach practical in real applications.

[28] C.-P. Tsai, C.-K. Yeh, and P. Ravikumar, “Sample-based explanations via generalized representers,” Advances in Neural Information Processing Systems, vol. 36, 2024.

[29] B. Schölkopf, R. Herbrich, and A. J. Smola, “A generalized representer theorem,” in International conference on computational learning theory. Springer, 2001, pp. 416–426.
Thanks to the authors for the detailed response! I have gone through the rebuttal and most of my concerns have been addressed. Based on the new results, it is evident that the pairwise model comparison might increase the computational complexity of the overall procedure, given that the authors claim the computational complexity is within an acceptable range. I will keep my rating.
Thanks for your quick response! 

Regarding the proof of Lemma 1: Sorry for the confusion. Let me elaborate my concern more. To correct the proof, you need to prove that the $\\epsilon$ in the contradiction is indeed $\\ge\\sup_x\\|f(x)-g(x)\\|$ as this is your assumption. 

Regarding novelty: I agree that this work provides a new perspective in imputation evaluation. However, my point is that it is a trivial combination of two ideas (downstream task evaluation + interpretation). I believe that this would be an interesting work if it were submitted to the NeurIPS Datasets & Benchmarks track, but I expect more technical novelty in the main track. 

Regarding L1: Sorry for the confusion. What I wanted to emphasize is that "it would be difficult to choose an appropriate kernel for a new dataset in practice." You do not have to reply to this one. As this limitation is not a limitation of your method, my evaluation of your paper was not considering this limitation. Instead, you should discuss this limitation in the paper.

Regarding L2: Sorry for the confusion. The examples you gave are not evaluating imputation methods; your examples are just evaluating each imputed value. In practice, merely evaluating each imputed value is pointless. When people imputation the time series, there is no way to evaluate because future data are not available. Thus, what we really need is to find out the best imputation method to use instead of just how good each imputed value is.
We greatly appreciate the question raised, which is very important for us to improve the quality of our work.
> Regarding W1

1. Proof of Lemma1

We sincerely apologize for the confusion caused and acknowledge that in real-world applications, $\epsilon$ cannot be arbitrarily small. However, the proof of Lemma 1 does not necessitate $\epsilon$ to be arbitrarily small. Instead, it demonstrates that for any given $\delta$ and $\epsilon_2$, there always **exists** $\epsilon$ (not necessarily 0) that satisfies the condition.

We suspect this misunderstanding might stem from our proof process. Let's revisit it. In the proof of Lemma 1, found in the appendix, we use $\epsilon \rightarrow 0$ to demonstrate the invalidity of the inequalities $|x-y| \delta \leq \epsilon$ and $N\left|\frac{2 \epsilon}{\delta}\right| \geq \epsilon_2$, thereby proving Lemma 1.

However, the validity of these two inequalities also depends on $\epsilon_2$ and $\delta$. In other words, we only need to ensure that $\epsilon < |x-y| \delta$ and $|\epsilon| <\left|\frac{\epsilon_2 \delta}{2N}\right|$. It is important to note that $\delta$ and $\epsilon_2$ are set by us. If $\delta$ and $\epsilon_2$ are not extremely small, there is no need for $\epsilon$ to be extremely small.


In practical applications, we don't need to require $\delta$ and $\epsilon_2$ to be sufficiently small since we only need to compare which of the imputation methods is more beneficial for downstream forecasting tasks at each time step. In this case, we only need to know the sign of the first derivative rather than an absolute accurate approximation. Both Table II in the original paper and our newly uploaded Table I in the attached PDF can prove the feasibility of doing so in practical operation.

2.

We acknowledge that our approach shares some similarities with [28]. However, the method in [28] **can not be directly** used to solve our problem. We reiterate our differences, that is, their approach focuses on the importance of training data, while we evaluate the benefits of replacing one imputation value with another in the time series for downstream tasks. 

To our knowledge, we are the first to propose such a viewpoint and provide practical approaches as the work related to the pre-trained model is simply evaluated using some metrics of downstream tasks. On the contrary, we provide a rigorous mathematical definition of our task (see equation 1 in the original paper, which is a bi-level optimization problem) and conduct fine-grained evaluations. 

In addition, **our method can also be applied to combine different series imputation methods to bring benefits to downstream forecasting tasks.** These are properties that the relevant methods do not possess and can not be easily transferred. In addition, we also provided an **accelerated calculation** method in Section 2.3 in the original paper, which plays an important role when the model parameters and the amount of data increase(see table in Regarding L1).

>Regarding L1

We acknowledge that there are general limitations to using RKHS approximation. After our experiments on the original ELE dataset with 15-minute resolution (about 100000 training samples and the predicted length is 96), the NTK kernel achieved a good balance between computational efficiency and accuracy compared to methods like the influence function in our task. (with the **acceleration calculation** in Section 2.3).

|   ELE   |  Mean | SAITS | Mean+SAITS+IF | Mean+SAITS+ours | Mean+SAITS+ours_seg16 | Mean+SAITS+ours_seg8 |
|:-------:|:-----:|:-----:|:-------------:|:---------------:|:---------------------:|:--------------------:|
|   MSE   | 0.249 | 0.307 |     0.248     |      0.238      |         0.236         |         0.237        |
| Time(s) |   -   |   -   |     126.33    |     1053.28     |         264.27        |        118.53        |

In addition, unlike [28], which focuses on explaining model outputs and necessitates high fitting accuracy, our method does not require such precision to deliver benefits, highlighting a key difference between the two approaches.

>Regarding L2

We are delighted to provide some relevant examples. The first example is that the importance is related to the value.  Suppose we have 10 generators, each with a power generation capacity of 10MW, totaling 100MW. Then filling 99.97MW as 100.02MW at the peak and filling 5MW as 5.05MW at the valley will result in the same aggregated error. However, due to the limitations of the generator, a 0.05MW error at the peak will require us to add a new generator to ensure power supply, resulting in significant additional economic expenses. The 0.05MW at the valley has no such concern. The second example considers the importance related to time or external events. When measuring a patient's blood pressure around the clock, the same margin of error in aggregated metrics can have different effects at different times or when patients engage in different physiological activities.
Thanks for your detailed response and new experiments. I still have a few concerns.

Regarding W1: I revisited the proof of Lemma 1. Unfortunately, Lemma 1 does not seem to be valid. The proof of Lemma 1 needs $\\epsilon\\to 0$. However, given two functions $f$ and $g$, the assumption that $\\|f(x)-g(x)\\|\\le\\epsilon$ implies $\\epsilon\\ge\\sup_x\\|f(x)-g(x)\\|$. As long as $f$ and $g$ are not identical, then $\\sup_x\\|f(x)-g(x)\\|>0$. This means that $\\epsilon$ cannot be arbitrarily small, falsifying the proof. Please either correct or remove Lemma 1, as wrong theoretical results should not be published. This is now my biggest concern.

Besides the incorrect Lemma 1, the three points you mentioned seem to be trivial extensions of [28]. Although these ideas are somewhat new in time series imputation evaluation, their perspective novelty is rather limited. For instance, evaluating the performance via downstream tasks is the current mainstream for pretrained language or vision models. Furthermore, none of the three points is of technical novelty as your proposed method is basically [28]. This work seems to be just treating the "interpretation" results as "evaluation" results. 

Regarding L1: Thanks for your detailed explanation. I agree that there have been related empirical evidence, but my point is that it is difficult to evaluate the goodness of approximation. [28] showed that different datasets need different kernels, so it would be difficult to choose an appropriate kernel for a new dataset in practice when we do not have ground truth data. I am not blaming the proposed method; just pointing out a general limitation of RKHS approximation. 

Regarding L2: The examples you mentioned are irrelevent to the evaluation of time series imputation methods. Could you provide an example about how the fine-grained information will be used to evaluate time series imputation methods without using an aggregated metric? 

Regarding W2 & L3: Your responses sound reasonable to me. Thanks for your detailed elaboration. Please discuss about them in the paper.