Adaptive Estimation and Learning under Temporal Distribution Shift

Experiments on real dataset

Please see the response to Reviewer c9dT under the same comment title.

How to do optimal wavelet selection.

This is a well-known model selection problem in wavelet denoising, not a drawback unique to our work. We acknowledge this in Section 7. Practitioners typically analyze data trends—if they follow a piecewise polynomial of degree $k$, a wavelet of order $k+1$ is chosen, as it effectively models such structures. However, higher-order wavelets introduce numerical instabilities and variance, though the latter can be mitigated with more data. Selecting a wavelet basis is akin to choosing a kernel in Gaussian Processes—application-specific and guided by practical intuition. A more comprehensive exposition can be obtained from [5].

Theoretical justification for using higher order wavelets

Please see the comment to Reviewer 4jCD under the same comment title.

Computational efficiency in comparison to prior work (Mazzetto and Upfal , 2023)

Applying the algorithm from (Mazzetto and Upfal , 2023) requires $O(\log n)$ calls to an ERM oracle. While our method requires exactly one call as detailed in Section 4. This is the reason for computational efficiency.

How to use Algorithm 1 for binary classification?

Please refer to the statement of Theorem 9. We compute refined the loss estimate of a model at most recent timestamp. For this, loss sequence $\ell(f(x_n),y_n),\ldots,\ell(f(x_1),y_1)$ is given as input to Algorithm 1. ERM is performed using the returned loss estimate.

Algorithm uses fixed soft-thresholding for wavelet denoising, which may not always be optimal

The guarantee in Theorem 2 is obtained via the universal soft-thresholding based denoising. The bound is known to be minimax optimal in light of (Mazzetto and Upfal , 2023).

The algorithm assumes that the ground truth sequence has sparse wavelet coefficients but it's hard to verify this assumption

We would like to point out that such an assumption is not made. Instead the bounds naturally adapts to the sparsity of the wavelet coefficient and hence the degree of non-stationarity of the groundtruth without any prior knowledge.

The temporal shift setting the algorithm considers is limited and not applicable to online learning...

We clarify that our setup is offline: given all data $y_n, \dots, y_1$, the goal is to estimate the ground truth at $t=1$. For online settings, Algorithm 1 can be used iteratively to update estimates over time. The multi-resolution nature of wavelets enables efficient updates with $O(\log n)$ complexity per round. Unlike typical online methods with cumulative regret guarantees, our approach provides stronger per-round estimation error guarantees.

Study omits comparisons with transformers and RNN architectures

Our method operates in a data-scarce regime. For example, to estimate the groundtruth at timestamp 50 we only have just 50 data points. Advanced models like transformers and RNNs excel with abundant data. Hence, we compare against methods suited for low-data settings that are also known to have robust theoretical guarantees on estimation error.

Results on adaptive tuning schemes for thresholding

As per reviewer’s suggestion, we conducted experiments on other well studied thresholding schemes namely SUREShrink [1] and Energy based thresholding [3] heuristic. The MSE results are reported below for Doppler signal and different wavelet types. Name in brackets indicates the thresholding scheme. Similar observations hold true for Random ground truth signal.

This experiment shows that some thresholding schemes may outperform universal soft thresholding empirically, but no single method is best overall. While SUREShrink is known to achieve minimax MSE optimality in non-parametric regression, its extension to point-wise bounds remains unclear. We will include these insights in the manuscript.

References

[1] https://www.jstor.org/stable/2291512

[3] https://digital-library.theiet.org/doi/full/10.1049/iet-smt.2016.0168

[5] https://www.sciencedirect.com/book/9780123743701/a-wavelet-tour-of-signal-processing

We thank the reviewer for their comments. Please see the responses below.

Experiments on real dataset

As an application of our proposed methods, we conduct a model selection experiment using real-world data. We evaluate our method on data from the Dubai Land Department (https://www.dubaipulse.gov.ae/data/dld-transactions/dld_transactions-open), following the setup identical to that of [2]. The dataset includes apartment sales from January 2008 to December 2023 (192 months). Each month is treated as a time period, where the goal is to predict final prices based on apartment features. Data is randomly split into 20% test, with two train-validation splits: (a) 79%-1% and (b) 75%-5%.

For each month $t$, we train Random Forest and XGBoost models using a window of past data where we consider window sizes $w \in [1,4,16,64,256]$, yielding 10 models per month. Validation MSEs from past and current months are used to refine the current month’s estimate of MSE via Algorithm 1 or ARW from [2]. The refined validation scores are used to select the best model for final MSE evaluation on test data. We report the average MSE of this model selection scheme over 192 months and 5 independent runs, comparing our method to ARW [2]. In the table below, HAAR and DB8 are versions of Algorithm 1 with the corresponding wavelet basis given as input.

Case 1: 79%-1% Train-Validation split for each month

Case 2: 75%-5% Train-Validation split for each month

We see that wavelet based methods shine especially when the validation data is scarce. This allows us to include more data for training while still allowing to obtain high quality estimates for validation scores. Such a property can be especially helpful in data-scarce regimes. Unlike the synthetic data experiments, here we find that Haar wavelets perform better than DB8. This can be attributed to the following facts: i) the noise in the observations depart from iid sub-gaussian assumption and; ii) the high degree of non-stationarity in the pricing data as indicated by Fig.6(b) in [2] makes the underlying trends to have a low degree piecewise polynomial structure. This makes the groundtruth irregular (or less smooth) which can be suitably handled by lower order Haar wavelets which are also less smooth and abrupt (see Fig. 2 in Appendix).

Lemma 4: please cite if it came from somewhere else or state that this is a previously unknown result

We have mentioned in the proof that Lemma 4 that it is adapted from (Achille and Soatto, 2018). For clarity, we will add the phrase “(adapted from Achille and Soatto, 2018)” in the Lemma statement itself.

Def 6: I do not understand the last line, could you explain that more clearly? With the T1=Ts if T = Tr, etc

Thanks for pointing this out! There is a typo and the sentence should be “Here T represents training data T≡T_r (or testing data T≡T_s) and T1 ≡ T_s (or T1 ≡ T_r).” We will revise the manuscript accordingly. We will also add the following explanation into the manuscript. “Specifically, under the Training Distribution Shift Scenario, the training data T≡T_r consists of (1) samples from the testing distribution T_1≡T_s, and (2) samples from a dissimilar distribution T_2. When, more and more dissimilarly-distributed samples from T_2 are added into the training distribution, this leads to higher dissimilarity ratio $\beta$. Under the Testing Distribution Shift Scenario, the testing data T≡T_s consists of (1) samples from the training distribution T_1≡T_r, and (2) samples from a dissimilar distribution T_2. When, more and more dissimilarly-distributed samples from T_2 are added into the testing distribution, this leads to higher dissimilarity ratio $\beta$.

Typos

Thanks for pointing out the typos. We will fix them. $H$ in Algorithm 1 should be $W$ the matrix for wavelet transform. The example below the definition of the TV class must be sum of squared differences as the reviewer pointed out. This is a typo that will be corrected.

References

[1] Adapting to Unknown Smoothness via Wavelet Shrinkage, David L. Donoho and Iain M. Johnstone, Journal of the American Statistical Association, 1995

[2] Model Assessment and Selection under Temporal Distribution Shift, Elise Han, Chengpiao Huang, and Kaizheng Wang, ICML 2024

We thank the reviewer for their comments. Please see the responses below.

Why DB8 can outperform Haar wavelets in the synthetic experiments

The main reason why higher order wavelets can outperform Haar is because, $k+1$-th order wavelets provide an ideal basis for sparsely compressing the information in a gorundtruth signal that has a piecewise polynomial trend with degree upto $k$. A sparse set of wavelet coefficients of the groundtruth signal leads to a low value for the bound in Lemma 1. We attempt to capture this intuition empirically as demonstrated in Fig.4 in Appendix and pointed the readers to that figure in the main text itself (L165, col 1).

Theoretical justification for using higher order wavelets

Wavelet basis of order $k+1$ is known to be an ideal basis for sparsely representing piecewise polynomials of degree (upto) $k$, or more generally groundtruths with low $k$-th order total variation. $k$-th order total variation measures the the variation incurred in the $k$-th order (discrete) derivatives of a groundtruth. It is known in the literature that low $k$-th order TV for a groundtruth is equivalent to sparse (measured in the sense of L1 norm) wavelet coefficients in a $k+1$-th order wavelet basis. For a precise quantification see Theorem 1 in [4].

Extension to the multivariate setting

We acknowledge that we study the estimation in a univariate setting. Extension of our point-wise estimation proposal to multivariate time series that takes into account inter-task correlations is an interesting direction to explore as a future work.

Experiments on real dataset

Please see the response to Reviewer c9dT under the same comment title.

Intuition behind achieving Lemma 1 and Theorem 2 without assuming prior bounds on $|\theta_i|$ or $|\theta_{i+1} - \theta_i|$

We first address informally why no prior bound on $|\theta_i|$ is required. The wavelet soft-thresholding estimator applies a soft-threshold at the level of $\tau := \sigma \sqrt{\log n}$. Suppose $\alpha_i$ is a wavelet coefficient of the groundtruth data. Consider the case when $|\alpha_i| > \tau$. Then the shrinkage caused by soft-thresholding only introduces a bias of at-most $\tau$. Similarly when $|\alpha_i| \le \tau$, then the bias is capped at $|\alpha_i|$ itself. This intuition can be extended to noisy wavelet coefficients upto constants using concentration arguments.

Eventhough we do not assume any prior knowledge on $|\theta_{i+1} - \theta_i|$, note that the bound in Theorem 2 naturally gets worse when there is lot of intra-sequence total variation. We remark that in the bound $\bar theta_{1:t} = (\theta_1 + \ldots + \theta_t)/t$. Hence higher intra-sequence total variation leads to larger values for the bound. The highlight here is that the bound (which is minimax optimal as per the results of Mazzetto and Upfal 2023) in Theorem 2 is obtained with no such prior knowledge on the intra-sequence variation. The reason for attaining such an adaptive bound is primarily algebraic as highlighted in the proof of Theorem 2.

References

[4] Minimax Estimation via Wavelet Shrinkage, David L. Donoho and Iain M. Johnstone, Annals of Statistics, 1998

We thank the reviewer for correctly recognizing and appreciating our contributions.

Rates in Theorem 11

The displayed rates are indeed the optimal rates for squared and absolute losses.

Paragraph below corollary 8

We will fix the typos and rewrite it for better clarity.

We thank the reviewer for their comments. Please see the responses inline.

Experiments on real dataset

Please see the response to Reviewer c9dT under the same comment title.

Experiments are based on ablation

Ablation on noise level was conducted primarily to test the validity of the algorithm for various signal-to-noise levels.

Existence of Linear Representation of Ground Truth

Wavelet basis are universal basis for $\mathbb{R}^n$. This means that they are of full rank and any ground-truth sequence in $\mathbb{R}^n$ can be written as a linear combination of the basis vectors. So the existence of such linear combination is not an assumption, but a structural property that always holds true.

Extending to other noise distributions

In this work we primarily focused on subgaussian noise, as the main motivation for our methods is to get high quality estimates of the groundtruth without any prior knowledge on how it evolves. This is common in many studies from non-paramteric regression (see for eg [6] for a comprehensive exposition). That said, extending to other noise distributions that are heavy tailed is a possible future direction of research.

Missing Constants in Lemma 1

We believe that we have defined all constants pertaining to Lemma 1 in its statement.

$H$ in Algorithm 1

Thanks for pointing this typo, instead of $H$, it should be the Wavelet Transform Matrix $W$ that is given as input to the algorithm.

Meaning of $\bar \theta_{1:t}$

It is the average of last $t$ ground truth values. $\bar \theta_{1:t} := (\theta_1+\ldots+\theta_t)/t$. We will update the manuscript to reflect this.

Assumptions for estimation

The only assumption we make is the fact that the noise in the observations are iid sub-gaussian as stated in Lemma 1. No further assumptions are made.

Connection to random process

Random processes such as Gaussian processes impose prior structural assumptions (which are often hard to verify) on how the ground-truth evolves over time. However, no such structural assumptions are placed while deriving our bounds in Lemma 1 and Theorem 2. Instead, the bound naturally adapts to the degree of smoothness in the ground-truth. We will mention this in the context of random process.

However, viewing our problem setup through the lens of GPs holds the potential to unlock new properties that are simultaneously connected to estimation quality and uncertainty quantification. This is indeed a good direction to explore in the future.

The proof is adapted from (Achille and Soatto, 2018)

In the paper, we have acknowledged in Appendix C that we adapt the proof from (Achille and Soatto, 2018). Note that (Achille and Soatto, 2018) did not make a distinction between the training distribution $D_{Tr}$ and the predictive distribution $p_{\theta}( y | \mathbf{x} )$. As the result, it is important to include our full proof so that avid readers are able to understand and recreate all of the mathematic derivations. The proof will also serve as a reference for researchers who want to improve our performance bounds in Theorem 7.

Infinite number of training data and perfect ML model training in Section 3

Our goal with Section 3 was just to act as a compelling motivation for Section 4. In particular, to formally show how distribution shift between training and test can degrade model performance. Though we deal with population level losses, extension to finite samples can be realized with standard concentration inequality arguments. Our focus was just to show formally that if the test distribution is a mixture of training and some other contaminant distribution, then the model performance degrades as the contamination fraction increases.

This sets the stage that performance degradation is inevitable even in an ideal case for perfect model training and it is important to develop algorithms whose performance degrades gracefully under distribution shift in practical use-cases. Note that to develop our algorithms in Sections 4, we already considered realistic scenarios with limited training data and imperfect machine learning training processes.

References

[6] Adaptive Piecewise Polynomial Estimation via Trend Filtering, Ryan Tibshirani, Annals of Statistics 2014