Empirical Bayes Trend Filtering Through a Variational Inference Framework
Authors: Dongyue Xie
Venue: Accepted by TMLR ()

Abstract: This paper introduces a novel framework for Bayesian trend filtering using an empirical Bayes approach and a variational inference algorithm. Trend filtering is a nonparametric regression technique that has gained popularity for its simple formulation and local adaptability. Bayesian adaptations of trend filtering have been proposed as an alternative method, while they often rely on computationally intensive sampling-based methods for posterior inference. We propose an empirical Bayes trend filtering (EBTF) that leverages shrinkage priors, estimated through an empirical Bayes procedure by maximizing the marginal likelihood. To address the computational challenges posed by large datasets, we implement a variational inference algorithm for posterior computation, ensuring scalability and efficiency. Our framework is flexible, allowing the incorporation of various shrinkage priors, and optimizes the level of smoothness directly from the data. We also discuss alternative formulations of the EBTF model, along with their pros and cons. We demonstrate the performance of our EBTF method through comprehensive simulations and real-world data applications, highlighting its ability to maintain computational efficiency while providing accurate trend estimation.

Published in Transactions on Machine Learning Research (03/2025)
Empirical Bayes Trend Filtering Through a Variational Inference Framework
Dongyue Xie
dynie@uchioogg.ell
Deportment of Statistics University of Chicago
Reviewed on OpenReview: https: openreview. net/ forun? id= AHT
Abstract
This paper introduces a novel framework for Bayesian trend filtering using empirical Bayes approach and a variational inference algorithm. Trend filtering is a nonparametric regression techmique that has gained popularity for its simple formulation and local adaptability, Bayesian adaptations of trend filtering have been proposed as an alternative method, while they often rely on computationally intensive sampling-base based methods for posterior inference. We propose an empirical Bayes trend filtering (EBTF) that leverages shrinkage priors, estimated through an empirical Bayes procedure by maximizing the marginal likelihood. To address the computational challenges posed by large datasets, implement variational inference algorithm for posterior compulation, ensuring scalability and efficiency. Our framework is flexible, allowing the incorporation of various shrinkage priors, and optimizes the level of smoothness directly from the data. We also discuss alternative formulations of the EBTF model, along with their pros and cons. We demonstrate the performance of our EBTF method through comprehensive simulations and real-world data applications, highlighting its ability to maintain computational efficiency while providing accurate trend estimation.
1 Introduction
Nonparametric regression methods have been widely used in many statistical applications such as spatial statistics, time series analysis, and survival analysis (Dabrowska) 1987; Gelfand & Schliep 2016 Moulines ct al. 2007). When the relationship between a predictor and a response variable is nonlincar, nonparametric regression can effectively capture the true underlying relationship by fitting a curve to the predictors. Classical nonparametric methods such as smoothing splines, B-splines, and kernel methods can be recast as penalized linear regressions and they are straightforward to fit (Härdle) 1990). However, they are not adaptive in the sense that they cannot adjust to local changes in the curve.
Trend filtering is relatively new method for nonparametric regression. It penalizes the differences of adjacent signals using an 4 penalty, and the penalty requires parameter tuning via cross-validation. Trend filtering was initially introduced as splines with higher-order total variation regularization (Steidl et al., 2006), without being named trend filtering. Later, trend filtering was independently introduced by Kim et al. (2009) as a modified version of Hodrick-Prescon (H-P) filtering (Hodrick & Prescott 1997), that changes the penalty from la to I norm. A more statistical and theoretical study of trend filtering provided by Tibshirani (2014). which showed that trend filtering achieves the minimax rate over smoothness class defined by bounded total variation, mainly due to its ability to choose basis adaptively from data.
One of the earliest Bayesian adaptations of trend filtering is from Roualdes (2015), and the author borrowed ideas from Bayesian lasso Parr Casella 2008). Other shrinkage priors can be placed on the differences of the signals. Examples are the spike-and- (George & McCulloch 1993). normal-gamma (Brown & Griffin 2210)), generalized double-Pareto Armagan al.] 2013), horseshoe prior (Carvalho al.] 2009), and scale mixture of normal distributions FFuulkner & Minin 2018 A dynamic shrinkage process and the

---

Published in Transactions on Machine Learning Research (03/2025)
corresponding Bayesian trend filtering based dynamic linear model are proposed by Kowal et al. (2019) All of the above Bayesian methods use Gibbs sampler for posterior inference. When the sample size is large, they suffer from high computational cost and low speed.
In this paper, we propose a novel framework for Bayesian trend filtering that is fast, locally adaptive, and accurate. Our method incorporates a shrinkage prior on the differences of the signal and employs an empirical Bayes method estimate the prior by maximizing the marginal likelihood. The posterior distribution is computed using variational inference algorithm. We highlight the advantages of our method as follows: 1. A fast and stable empirical Bayes trend filtering (EBTF), applicable to large-scale datasets; 2. Flexible shrinkage priors that adapt to the best shrinkage operator, not limited to the I penalty; 3. Learns the level of "penalty" from data by optimization; 4. Naturally extends to more complex settings, such as sparse signal, which will be studied in Section 5. We provide a Python implementation of the method, and all the code and analysis are available the package ebtfPy on GitHub (https://gi wareluu...oom//oooeeo/eeeoeee hoa//DDnngeeeeeebbtF2 All the detailed derivations in this paper are in the Appendix D
Denote the flat prior improper prior for 0 0 ~ C(-), and its density function is p(0) over the support of 0. A vector is denoted in bold such B, and when need its elements, the vector is denoted as B = (3; = 1.2,...nn. A diagonal matrix with its diagonal elements is denoted as W diag(w; w; : 1,2, n..
Notation:
2 Empirical Bayes Trend Filtering
In this section, we first give a brief review of the trend filtering problem, then present our models and algorithms. For given integer E N, the kth order trend filtering finds
β arg min BER _ ~ - (K+ B111,
where D(++1) E R(a-A-1)xx is the discrete difference operator of order & + 1. When k = 0, the estimated sequence is piecewise constant, and the difference matrix is
1
0
0 0
0
-1 1
0
0
D =
( 1)xn
(1)
0 0
-1 1
For l 2 1, the difference matrix is defined as (+++1) = (1)) ((),, where D(1) is the (n - k 1) x (n k)) version of equation [ When k = 1, the estimated sequence is piecewise linear; and bwice-differentiable when ==
Existing algorithms for solving the trend filtering problem include primal-dual interior point (Kim al, 2009), path algorithm (Tibshirani & Tay 2011] and alternating direction method of multipliers (ADMM) [Ramdas & Tibshirani 2016). The parameter À controls the smoothness of the estimated sequence and is often selected using cross-validation. Cross-vall validaaiii for trend filtering is typically not random as the folds are often fixed. For a detailed description, see the R function cv. trendfilter (Arnold & Tibshirani 2016
2.1 Model formulation
We formulate trend filtering model as dynamic linear model (DLM) data-generating process. Throughout this paper we will focus on the first-order sequence model, and we shall see that the general order model formulation is straightforward by using the corresponding-orde difference matrix. We consider the following Bayesian variants of first-order trend filtering,
v:/8: ~ N(3,, N(3,,03s?) for i i= 1,..,
B1 ~ C(.),
(2)
(Bi+1 - B3) ~ g()), for j = 1,.... n -1
2

---

Published in Transactions on Machine Learning Research (03/2025)
where C()) denotes the uniform distribution over the entire real line, g(-) is the prior on the difference between two consecutive means, o² is the unknown random error variance and o is the known heterogeneous variance term. The variance oe can be regarded as the inverse weight for each observation, and for the homogeneous sctting sf = 1. The model formulation can be equivalently expressed in a matrix-vector form as
y/B o)),
81 C();
(D8), ~ g(())
where = (vs : = 1,2,..., n). 8 (8, : 1.2.....n) S diag(s;), and D is the first order difference
matrix defined as equation [
2.2 Choice of prior
We choose shrinkage prior 9(-) such that the difference between two neighboring signals are mostly small (close to 0), which would lead to a spatially-structuuee signal. In our model formulation, the choices of shrinkage priors are flexible. In this paper, we focus on mixtures normal distributions including the point-normal and adaptive shrinkage (ash) prior (Stephens) 2017). The prior is represented as
g Σ (,,
(4)
where Σx Th and o² is the random error variance in model equation 2 For point-normal prior (spike and slab prior with normal components). =2 2, and of is often fixed at 0 a very small number, while (T1,0)) are the hyperparameters. For ash prior, K is often large, and all of are known and fixed, and span a large grid (from a small value to large ones), while TH 's are the hyperparameters. In this paper, we make a novel extension of the ash prior such that of are not fixed excepting the first one (the one for the spike component). Both priors have been applied to mean estimation and inference (Castillo Roquain 2020 Willwerscheid & Stephens 2021), matrix factorization (Ning & Ning 2024 Wann & Stephens] (021), sparse
regression (Kim al.) 2024 Ray & Szabó, 2022) and wavelet denoising (Chipman et al., 1997 Xing et al. and they have been shown to have better performance over other choices of priors.
2021
2.3
The variational algorithm
For hyperparameters in the prior, can either fix them before model fitting, or learn them from data. In this paper, we take the empirical Bayes approach for estimating the prior, and the posterior inference is conditional on the estimated prior distribution. While it's possible to use MCMC for sampling from the posterior, it's intractable for large-scale datasets QQuiroz et al.) 2019). We instead propose to use variational inference for posterior computalion. We start this section with a high-level review of empirical Bayes and variational inference.
2.3.1 Review of empirical Bayes and variational inference
An empirical Bayes (EB) approach estimates the prior by maximizing the marginal likelihood p(y:g) = J "py (8)) then the posterior is computed conditional on ğ as p( y, g). EB has been widely applied large-scale multiple testing (Eiron 2004; Stephens
2017 Liu et al. 2024 Xie 2025 shrinkage estimation
[Johnstone Silverman, 2004 Koenker & Mizera 2014 Dey et al.
2018 and other tasks such sparse
regression and model selection Martin et al. 2017 Kim et al.) 2024
Zou et al.
2024
Variational Inference (VI, Blei et al. (2017) turns the posterior inference problem into an optimization
problem by approximating the true posterior with a more tractable distribution. The VI finds
(9( arg min DKL DKL
VEe
3

---

Published in Transactions on Machine Learning Research (03/2025)
where Q is family of approximate densitics, and DE is the Kullback-Leibler (KL) divergence. In practice, maximize the Evidence Lower Bound (EIBO), which is lower bound the poyp(y:))
F. y) log g( DxL (
Ra (logply. ply. 9) g) log (()))
Variational empirical Bayes (VEB) combines variational inference and empirical Bayes in a single optimization problem, expressed as
q (B), 0 = arg max ()).
qEQ.sEG
2.3.2 Variational empirical Bayes for trend filtering
We develop the variational inference algorithm for model equation 2 in this section. For the prior in equation we follow the standard approach for the Gaussian mixture model and introduce the latent variable 2 such that
(az+1 - 8,)1=zk =1 ~ N(0,0k0²) p(z}* = 1) = Tk, for j = 1,2, 11
(5)
For the variational posterior, we consider the following variational distribution class that factorizes over / and
K
(B, = (3)))==()))) = N(8 V
=
II
(6)
k=1
where Ojk Izza (2jk- 1) is the posterior probability that (3,+1 3g) is drawn from the kth mixture component. When V is diagonal matrix, the posterior distribution is fully factorized. But do not make such simplification and assume V is a general covariance matrix. The evidence lower bound for model equation is then
FEBTF E.g log ((g)(BB) + E. 10, (oz)) E, log ((3
For the optimization of the ELBO, we take the VEB approach introduced in the section 2.3. since o², π, and prior variances (a are unknown, we treat them as parameters, which are optimized when maximizing the ELBO. The coordinate ascent algorithm for fitting EBTF model equation 2 has the following updates:
Given qa, the update of the posterior probabilitics ajk
TA -DVD)
K ΣE + o2o))
2. Given az the update for postcrior variance and posterior mcan are 3= Vujo².
Given qa iz, the update for the prior variances of and prior probabilities π are
o E, (n (DVD"))))
Th x Σash
o² o²>, agk
4. Given the rest, let Ω S-1 - 1TTH WD, update o² as
0² = "= 2y 5 13 + TT (r(()))) - 1).
4

---

Published in Transactions on Machine Learning Research (03/2025)
Algorithm VEB algorithm for fitting EBTF equation 2 (outline only)
Data yi, variances s2, for i == 1,,,,.... n. Posterior mean B, posterior precision matrix diagonal d and super-diagonal e, residual variance o², prior weights TH and variances of for k 1.......K. K. repeat
Input:
Init:
1. Update posterior weights aik for i 2, 2,3, ..... and k = 1, K
2. Update posterior precision matrix (its diagonal and super-diagonal elements only);
3. Update 8 by solving a (tridiagonal) banded linear system;
4. Update prior weights, variances, and residual variance.
until converged
The posterior precision matrix V-1 (8-1 + D'WD)/ is tridiagonal matrix, because of the tridiagonal structure of the difference matrix D. This indicates that the sequences (8x) are conditionally independent a posteriori given two adjacent variables. Specifically, for the posterior distribution 9a., the adjacent two neighbors 98_, and g8.t1 are all the information nceded to determine qa.
Although the updates are formulated in matrix multiplication form, the computation cost can be significantly reduced by leveraging the special structure of the difference matrix. The matrix D'WD is tridiagonal and can be caleulated fast, by operations only on wi, yielding its diagonal and super-diagonal elements. An optimized banded system solver (such scipy. scipy.linalg www.lliall....m swwveh_lllehseeee banded) can be used to find B. To find the diagonal of DVDT, the diagonal and super-diagonal elements of V are first obtained by inverting the tridiagonal precision matrix using the recursion algorithu from Usmani 1994). Then DVD D" can be directly calculated using operations only on the diagonal and super-diagonal elements of V. The final algorithm for implementation is summarized in Algorithm [ The algorithm iteratively solves for the maximum of each parameter while keeping all other parameters fixed, ensuring that every update increases the objective function.
Remark 2.1. In variational inference, the ELBO is in general not convex, and the initialization is important for non-convex optimization problems. However, there are fast initialization methods that can provide good starting point for the variational algorithm. We address these initialization issues here. The residual variance o² initialized by applying median absolute deviation (MAD) to the finest level of wavelet coefficients as described in section 4.2 of Donoho Johnstone (1994). For heterogeneous variances where sf are unknown, we may use the running MAD estimator proposed in Gao (1997) the wavelet-based variance estimation in Xing et al. (2021). The precision matrix is initialized to the identity matrix. The posterior mean is initialized to the wavelet denoised mean by applying soft thresholding at a to the Haar wavelet. cocfficients. The wavelet method is chosen because it is fast and the Haar wavelet threshold provides piecewise constant signal estimation. An alternative initialization ion for the posterior mean 3 is cross-validated trend filtering, which generally scalable and provides reasonably good starting value. We compare the two strategics in the simulation study and find that they yield similar results, as indicated by the comparable RMSEs shown in Figure 11 in Appendix A In contrast, if B is initialized at the data points y, the method typically fails to converge, causing the learned curve to remain stuck at y.
3 Alternative formulations of the EB trend filtering
In this section, we present two alternative formulations of the EBTF model, and show the equívalence among all formulations in terms of the objective function ELBO. We further discuss the pros and cons for each formulation.
5

---

Published in Transactions on Machine Learning Research (03/2025)
3.1 Multivariate normal variance prior formulation
The primary model equalion 2 can be formulated in an equivalent way as
y|3 N(B,o2))
D3 ~ NNO.
(8)
W ~ g(-), for j - 2.., -,,
where W is diagonal covariance matrix W diag( w ; 1,
1), and is prior the variances. Specifically, the ÿ corresponding to the prior equation 4 is wi Discrete(of. ..., oi: π), where the discrete distribution is defined as p(wy = o = for k == 1,2, K, with Σx T = 1. This formulation is named the multivariate normal variance (MNV) prior approach, as it introduces a multivariate normal prior on DB, followed by another prior on the variances. We use the VEB framework for prior estimation and posterior computation by maximizing the ELBO
FMN Egg log : g) LEg log q(8, W).
The VEB updates for maximizing the ELBO FNN are the same as those in Section 2.3.2 for fitting the EBTF model oquation if the variational postcrior distribution for model equation is chosen as
(B8, W) = gsw ()) = N(B; V) II (m
(9)
Although the VEB updates remain the same, the VEB algorithm (in Appendix C of the MNV approach has a nice property: alters between a simple update on qa, and general empirical Bayes Gaussian variance (EBGV, Appendix B problem for (ğ, qw). Hence, the MNV model formulation and inference modular. To accommodate different prior distributions problem for these priors, instead of re-deriving the full variational updates. The EBGV solver can then be plugged into the general variational inference iterations.
W. it is sufficient to develop the corresponding EBGV
3.2 Multiple linear regression formulation
The trend filtering problem can be formulated as penalized multiple lincar regression (MLR) problem, as shown in Lemma 2 of Tibshirani (2014) Specifically, let H E Ruxn i be the "inverse" of the first-order difference matrix D, such that D3= Hb = B, where (81, bT)T.
The first element of is 31, which can be regarded as the bascline value, and all the subsequent signals are additions or subtractions it. The vector b captures the piecewise difference among the remaining signals. Thus, model equation 2 can be reformulated as Bayesian sparse multiple linear regression problem:
b ~ N(Hb, o²S),
81 ~ C(-)
(10)
by ~ g(-) for j 1, 2, ..., n 1,
where the prior g(-) is sparsily-induccig and is the same equalion 4
If the variational posterior distribution for model equation 10 is chosen as
n-1
K II
q(b. =) = (b)g (b: 6, v1))
a j
(11)
k=1
then the ELBO for the the multiple linear regression formulation is
FVL E, logr(y))) b) + E. log z) - E., log (b) E. log gz (=(),,
and FMLR is equivalent the ELBO equation for the primary model formulation equation VEB updates for both models the same.
Hence the

---

Published in Transactions on Machine Learning Research (03/2025)
Given the equivalence of the three model formulations equation 2 equation 8 and equation 10 comment their advantages and disadvantages. The multiple linear regression formulation is very general, and there is a large number of methods for Bayesian sparse linear regression. Hence those methods can be readily applicd to Bayesian trend filtering with minimal modifications (though the availability of EB sparse regression methods limited). The multivariate normal variance prior formulation is modular and can easily incorporate different types of priors.
However, from a modeling perspective, the MNV and MLR approaches are specific to the pre-defined trend filtering problem and are not easily generalized to more sophisticated models. On the other hand, the primary dynamic linear model formulation of the trend filtering is more flexible in terms of adding additional model components and adding custom features. In Section 5 we illustrate this perspective by constructing and solving a sparse and spatially-structuree sequencing model. See the related discussion Section 7.2 of Tibshirani 2014
4 Numerical Examples
In this section, we evaluate the performance of our proposed method, EBTF, against other widely used non-parametric regression methods. All experiments are conducted on a Linux system with an 19-10900F processor and 32GB memory. The compared methods are:
1. genlasso-tf: cross-validated zeroth-order trend filtering using the R function ev. www.drreeelee... from the genlasso package (Arnold & Tibshirani 2014
2. wave-hard and wave-Bayes, Haar wavelet denoising (hard thresholding, (Donoho & Johnstone 1994 and Bayesian adaptive shrinkage (Chang et al., 2000)) using the Python function denoise wavelet from the skimage package (Van der Walt et al. 2014
susic-tf:
an empirical Bayes variable selection method extended for trend filtering, using the R function susie trendfiltte ilter from the susieR package (Wang al. 2020). We consider three settings where L 10, 20, 30 in the simulation.
4. BTF:
Bayesian trend filtering with dynamic horseshoe (DHS) prior and normal-inverse-gan (NIG) from the R package dap (Kowal et al. 2019). We use the default setting with 1000 MCMC burn-in iterations and 4000 final MCMC iterations for posterior calculations. We have found that setting = 0 (zeroth-order) generally leads poor fitting and convergence issues, so we set D in the btf function.
prior
5. GP:
Gaussian process regression with the Matern 3/2 kernel using the GaussiaanProeesseereesoo from scikit-learn. When fitting the Gaussian process, we first standardize the data before model fitting and then convert the fitted mean back to the original scale. As noted by one of the reviewers, this approach improves performance and speeds up GP runtime.
4.1 Simulation
We consider six different signal functions: blocks, steps, bumps, Gaussian density (Gauss), linear, and Heavisine. These functions are illustrated in Figure We set the number of samples to be n 1024, and the residual variance to o² 1. The signal-to-noise ratio (SNR) is defined as SNR = Var(8)/o² and is set to 3. Each experiment is repeated 20 times, and we report the averaged root mean squared error (RMSE), mean absolute error (MAF) and the coefficient of determination (RR2) between the estimated and the true signal along with their standard errors. The metrics are defined as
RMSF 111--111 1112, MAF 1 31, R² =1 1
Σ
1= (y/:
11
Σ
(1/:
where for Bayesian methods, A is the posterior mean.

---

Published in Transactions on Machine Learning Research (03/2025)
EBTF
genlasso-tf
Gauss
wave-Bayes
Gauss
wave-hard
susie-tf-10
susie-tf-20
susie-tf-30
GP
Step
BTF-DH5
Step
BTF-NIG
HeaviSine
HeaviSine
Blocks
Blocks
Bumps
Bumps
Linear
Linear
0.00
0.25
0.50
0.75
1.00
1.25
0.0
0.2
0.4
0.6
0.8
1.0
RMSE
R2
Figure 1: Simulation results: RMSE and R² scores of the competing methods across six signal functions. Each bar represents the metric value, while the horizontal lines on the bars indicate the standard error over 20 replications. The left panel presents the Root Mean Squared Error (RMSE), where lower values indicate better performance, while the right panel shows the R² metric, where higher values indicate better fit. Different colors correspond different signal functions, shown in the legend. The results for the Mean Absolute Error (MAE) follow similar trend to RMSE and are provided in Figure g in Appendix A

---

Published in Transactions on Machine Learning Research (03/2025)
Figure 1 presents the RMSE and R scores for all methods across different signal functions. The proposed EBTF mot generally achieves the near-lowest RMSE, and the near-highest R scores. Other strong performers include BTF-DHS,, genlasso-tf, wavelet demoising with hard thresholding, and Gaussian process, which all demonstrate robust performance across various signal functions. Notably, GP performs best on the Gauss and Heavisine signals, benefiting from their smooth nature. The BTF with DHS prior achieves the lowest RMSE a couple of signal functions blocks, bumps and lincar. The DHS prior also consistently outperforms the NIG prior in BTF. EBTF and BTF-DHS have consistently high Rt² scores in all the signal functions, particularly for the bumps function.
The susic-tf-30 method provides the best estimations for blocks and step signals, while EBTF, BTF-DHS also perform well on these piecewise constant functions. Increasing the L parameter in susie generally improves performance, particularly from L = 10 to L = 20, although further increasing it to L = 30 yields only marginal gains. However, susic-tf struggles with the bumps function, showing severe underfitting as the fitted curve often missing most of the bumps.
Figure 2 provides high-level summary of the simulations, illustrating the trade-off between estimation accuracy and runtime. Wavelet-based methods are significantly faster than the rest of methods, as they leverage the efficient pyramid algorithm for wavelet decomposition and reconstruction. EBTF offers consistently strong performance while maintaining a fast runtime; its overall Rº is only slightly smaller than BTF-DHS yet it runs faster than almost all methods except sewelet-based approaches. While susie-tf-20 achieves lower RMSE than susie-tf-10 and wavelet methods, this comes at the cost of substantially higher runtime. Increasing susie's L from 20 to 30 yield slightly performance gains but also results in an increase in computational cost. BTF-DHS achieves the lowest RMSE and highest R² scores overall but at the expense of significantly longer runtime. In contrast, BTF-NIG has a runtime comparable to susie-tf-20 but with lower performance.
4.2
Real data
In this section, we show the applications of EBTF to several real datasets. The first dataset, motorcycle acceleration (Silverman 1985 is classical example used for illustrating the nonparametric regression methods. It provides measurements of head acceleration in simulated motorcycle accident, used to test crash helmets. For comparison, we added fitted curves from the other methods (for susic-tf we used L = 20). The black curve in Figure 3 is the EBTF fit, and clearly il captures the trend of the acceleration over time. On the other hand, both the genlasso-tf and susie-tf exhibit some degree of fitter-.iiing..... The genlasso-tf seems to overshrink the signal around time 60 and time 100. also underestimates the signal around time to 20, as the estimated signal is clearly below all the observations in that time period. The susie-tf seems to underfit the signal in the time period from 70 to 90, as it only produces one big jump there. The GP provides smooth overall fit and captures the acceleration trend, while the wavelet-based method exhibits abrupt jumps and underfits the region around time period 60.
We next evaluate the performance of the methods on datasets with longer sequences and more abrupt changes in trends. The three datasets used in our experiments, sourced from Wu et al. 2021) are described as follows:
This dataset consists of time-series data collected from electricity transformers, including oil temperature, recorded at an hourly frequency between July 2016 and July 2018. The total sequence length is 17420.
ETTh1 dataset:
Illness dataset: This dataset contains 966 weekly records of influenza- like illness (ILI) cases in the United States between 2002 and 2021. captures the ratio of ILI patients to the total number of patients seen.
Weather dataset: This dataset includes CO2 concentration measurements recorded every 10 minutes throughout 2020. The dataset has more than 50000 records.
To visualize the fitted results for the ETTh1 datasot, we select three 500-length segments, as shown in Figure From Figure | we observe that genlasso-tf and susic-tf (L = 500) exhibit noticeable underfitting, while the EBTF, GP, and wavelet-based methods demonstrate better trend-capturing ability, with EBTF and GP
9

---

Published in Transactions on Machine Learning Research (03/2025)
0.98
0.96
0.94
R2
0.92
0.90
0.88
6
-4
-2
0
2
4
log(run time)
(a) R vs. runtime
0.45
0.40
RMSE
0.35
x
0.30
0.25
0.20
-6
-2
log(run time)
BTF-DHS
genlasso-tf
susie-tf-30
BTF-NIG
susie-tf-10 susie-tf-20
wave-Bayes
EBTF
wave-hard
GP
(b) RMSE vs. runtime
Figure 2: Simulation results: plots of log-transformed runtime vs performance metries. The top panel (a) shows the relationship between nog-transformed runtime and R² while the bottom panel (b) displays log-transformed runtime versus RMSE. Each point represents a different method, and the legend indicates the corresponding method labels. The runtime measured seconds and log-transformed for better visualization. Both metrics and runtime values averaged across all signal functions and repetitions for each method. The results for the MAE follow similar trend to RMSE and are provided in Figure 10 in Appendix A
10

---

Published in Transactions on Machine Learning Research (03/2025)
50
Acceleration
data
-50
genlasso-tf
susie-tf
GP
-100
wavelet
EBTF
0
20
40
60
80
100
120
Time Points
Figure 3: Motorcycle acceleration data from Silverman
1985
capturing the most intricate trend details. EBTF and wavelet also appear to capture the trend of outliers, as shown in the top plot of Figure 4 Future work could explore improving robustness.
The fitted curves for the illness and weather datasots are presented Figure 5 For the illness dataset, all methods capture finer seasonal variations, while susie-tf shows some underfitting at peaks, and wavelet exhibits abrupt jumps. The remaining three methods: EBTF, genlasso-tf, and GP, show very similar signal reconstructions. Similarly, susie-tf also underfits regions (such around 700 and 900 minutes) with abrupt changes in the weather data. Increasing L could potentially improve performance but at significant computationns cost. The other methods show very similar trends. Overall, the results highlight EBTF as an effective method for capturing both smooth and abrupt trend changes, while GP and wavelet-based methods provide reasonable approximations.
In addition, we evaluate the methods on holdout test data across the ETTh1, illness, and weather datasets. Specifically, for each dataset, we perform K-fold holdout with structured fold assignment (Arnold & Tibshirani 2014h. In this setup, every Kth data point is assigned to the same fold based its sequential order, ensuring consistency across runs. The first and last points are excluded from all folds and always contribute to building the predictive model. The holdout set error is computed by averaging the predicted values from the two neighboring points of each test sample, which belong to different folds. This approach maintains the temporal structure of the data while providing stable and reproducible results.
We compare the RMSE ratio relative a simple mean predictor (i.e., predicting using the mean of the observed data points) and present the results in Figure Across all datasets and data sizes, EBTF, genlasso-tf, and wavelet-based met hods achieve similar predictive performance, with relatively low RMSE ratios, and EBTF has the smallest RMSE when n = 4096. All methods perform well on the illness data, while on the weather data, GP and susie-tf have similar RMSE ratios that are slightly worse than the other methods when n = 1024.
Figure D also presents the log-scaled runtime comparison across mothods. The wavelet-based method the fastest, while EBTF demonstrates substantial runtime efficiency improvements over genlasso-ff, susie-tf, and GP, making it a competitive choice for large-scale time series analysis. Additionally, EBTF exhibits good scalability, as its runtime increases only marginally when increasing from n 1024 to = 4096, whereas GP shows larger increase in computational cost. Despite susie-tf's underfitting issues, it remains the most computationally expensive method across all datasets. For longer sequences, susie-tf requires a larger I to capture all trends effectively, but this comes at the expense of substantially increased runtime.
Overall, the results suggest that EBTF provides a strong balance between predictive accuracy and computa- tional efficiency, outperforiing GP and susie-tf in both accuracy and runtime while maintaining comparable
11

---

Published in Transactions on Machine Learning Research (03/2025)
27.5
25.0
22.5
Oil Temperature
20.0
17.5
15.0
12.5
10.0
2000
2100
2200
2300
2400
2500
Time Lag (hours)
20
Oil Temperature
16
14
10
10000
10100
10200
10300
10400
10500
Time Lag (hours)
10
Oil Temperature
8
6
15000
15100
15200
15300
15400
15500
Time Lag (hours)
data
genlasso-tf
wavelet
GP
susie_tt tf
EBTF
Figure 4: Fitted results for three selected 500-length segments from the ETTh1 dataset. The genlasso-tf and susie-tf (L = 500) methods underfit the trend, failing to capture the oscillatory wave patterns accurately. contrast, EBTF captures most of the trend details, while GP and wavelet-based methods provide reasonable approximations.
12

---

Published in Transactions on Machine Learning Research (03/2025)
1e6
data
genlasso-tf
wavelet
1.4
susie_tt tf
EBTF
Total Patients
1.0
0.8
0.4
200
400
600
800
1000
Time Lag (weeks)
(a) The Illness dataset represents the number of total patients with ILI recorded weekly. All methods are able to capture finer seasonal variations, while susie-tf shows under-fitting peaks and wavelet shows some abrupt jumps. The rest three methods have similar reconstruction of the signals.
data
genlasso-tf
460
wavelet
CO2 Concentration
450
420
200
400
600
800
1000
Time Lag (minutes)
(b) The Weather dataset records CO2 concentration every 10 minutes. We select the first 1000 data points to visualize. All methods, except susie-tf, effectively reconstruct both smooth and abrupt variations, while GP appears to better capture small fluctuations, particularly around the 600-minute mark.
Figure 5: Comparison of fitted trends for the IlIness and Weather datasets.
13

---

Published in Transactions on Machine Learning Research (03/2025)
performance to wavelet and genlasso-tf. Future work could explore optimizing EBTF further for even greater
robustness and scalability.
Finally we applied the EBTF method to eight more real datasets, and the figures for the fitted signals are
shown in Appendix A The original data were processed by Van den Burg Williams/ (2020), which evaluates several change-point. detection methods. For comparison, we included genlasso-tf fitted signals (blue line)
all plots. Overall, EBTF (black line) provides more visually appealing signal fitting, especially in its ability to capture the data trend without overfitting.
5 Sparse Empirical Bayes Trend Filtering
The primary EBTF model (2) can be viewed general generative prior for a spatially-structured sequence and can be applied in various dynamic model settings. For example, inducing sparsity the signal (Koop & Korobilis 2023 Hamiroz-Hassan
2020 Rockova & McAlinn 2021 A slight modification of the prior lcads
to sparse empirical Bayes trend filtering (sparse EBTF) as lollows:
y|3 ~ N(B,o25),
81 ~ NON, + (1-m0)C(())
(12)
To NO,
+ Σ
for j = 1,,,,.... n1 11,
where of is a pre-chosen small variance value such that N(0, a²ağ) spiky, with Ek=o T = We have added an extra mixture component that induces sparsity on the sequence 8 directly. In particular, each element of the sequence is now mixture of two components: one that promotes sparsity in B, and ceaehhness-indueiing component. The first component is spiky normal distribution at 0 that shrinks & towards O, while the second one is the same as in the original trend filtering model formulation. For the posterior, we again consider the following variational distribution class that factorizes between B and z, as
22 K
- = N (8:33V
2M &
(13)
where askk = qzn (2t = 1) is the posterior probability indicating the mixture distribution. We define ai := for = 2,3,..., , since the prior of B1 has only two mixture components. The detailed development of the VEB update for sparse EBTF is given in Appendix D.4
To illustrate the effect of sparse EBTF estimating sparse and spatially-structured signals, present example where sparse EBTF shrinks the sequence lowards 0 when the underlying signals are Iruly sparse. We generated = 4096 samples from the bump function, in which the signals are mostly at 0 and occasionally jump to large valuos. We fitted a sparse EBTF model to the data and compared the fit with the regular EBTF (without sparsity induction on the sequence). As shown in Figure 7 sparse EBTF is clearly able to shrink the estimated signals towards while estimating the spatially-structured curve. However, without the sigual-sparsity constraint, the estimated signals in sparse areas could be clearly non-zero, especially in areas between two spikes.
6 Extensions and Discussions
In this paper, propose fast and scalable empirical Bayes trend filtering method for nonparametric regression. The method leverages empirical Bayes estimation and variational inference, allowing it to learn the unknown smoothness level from the data. We demonstrated the superior performance of the EBTF method through simulations and real data examples. Our proposed variational posterior family is multivariate for the signal. This approach offers the benefit of fast computation while also maintaining the posterior dependency among all signals. An alternative posterior family to factorize over observations but not over 2, as q(,3,² 2) II, ((,122)) (2)))(z))) However this posterior assumes independence a posteriori and presents more computational challenge we need to track the posterior of over K components.
14

---

Published in Transactions on Machine Learning Research (03/2025)
0.30
4.0
0.25
3.5
3.0
Run Time (log)
RMSE Ratio
0.20
2.5
0.15
2.0
1.5
0.10
1.0
0.05
0.5
0.00
0.0
ETTh1
weather
illness
ETTh1
weather
illness
EBTF
genlasso-tf
wavelet
GP
susie-tf
(a) Holdout evaluation results with 1024. The illness dataset has 966 records.
0.25
5
0.20
Run Time (log)
RMSE Ratio
0.15
3
0.10
2
0.05
0.00
ETThl
weather
ETTh1
weather
EBTF
genlasso-tf
wavelet
GP
susie-tf
(b) Holdout evaluation results with 4096. The illness dataset has 966 records it is not shown in this plot.
Figure 6: Comparison of RMSE ratio and runtime across ETThl, weather, and illness datascts in the holdout test. For each dataset, we select the first n 1024 data points (top figures) and n = 4096 data points (bottom figures) for evaluation, except for the illness dataset, which contains 966 records, and we use all available data. The left panel shows the RMSE ratio relative to a mean predictor, where lower values indicate better predictive performance. The right panel presents the log-scaled runtime comparison. EBIF achieves a balance between predictive accuracy and computational efficiency, demonstrating improvements genlasso-tf, GP, and susie-tf.
15

---

Published in Transactions on Machine Learning Research (03/2025)
80
80
60
60
20
20
-0
20
Figure 7: Illustration of sparse EBTF. The true signal a bump function, as shown in Figure This signal is mostly sparse at O, while occasionally jumps to a large value. The left plot shows the regular EBTF fit, without inducing sparsity on the signals. The right plot shows the sparse EBTF fitted signal. The blue dashed line indicates y 0.
In our approach, we have utilized flexible non-parametric shrinkage prior on the differences. However, there are other Bayesian shrinkage priors that have been proposed recently, such as the widely used global-local shrinkage priors. Our flexible framework can easily incorporate different shrinkage priors, allowing for the selection of the appropriate prior based on the ELBO. The combination of shrinkage priors and variational inference presents an interesting direction for studying uncertainty quantification, particularly in the empirical Bayes setting (Xie & Stephens 2022: Ignatiadis & Wager 2022).
It is straightforward to extend our method to higher-order trend filtering by replacing the difference matrix D with higher-order matrices, and all the results still hold. One difference that when developing software implementations, we need to develop solvers and matrix manipulations for general banded matrices. For example, for k- 1, the posterior precision matrix is pentadiagonal.
As a final note, in real applications, the data may not always be real-valued: for non- non-Gaussian data, count and binary data are the two most commonly encountered types. For example in image denoising (Luisier et al. 2010), the pixel values are typically integers and we may assume they follow Poisson distribution. Variational inference methods have been developed for non- Gan-Gassiaan likelihood by leveraging a Gaussian-based model, and our method can be easily adapted to handle these non-Gaussian data types. For example, see Seeger & Bouchard (2012): Xie (023) for Poisson data, and Durante Rigon (2019) Jaakkola & Jordan (1997) for binary data.
References
Artin Armagan, David B Dunson, and Jaeyong Lee. Generalized double Pareto shrinkage. Statistica Sinica, 23(1):119, 2013.
Taylor B Arnold and Ryan J Tibshirani. genlasso: Path algorithm for generalized lasso problems. R package version, 1(3), 2014.
Taylor B Arnold and Ryan J Tibshirani. Efficient implementations of the generalized lasso dual path algorithm. Journal of Computational and Graphical Statistics, 25(1):1-27, 2016.
David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: A review for statisticians. Journal of the American statistical Association, 112(518):859-81)) 877, 2017.
Philip J Brown and Jim E Griffin. Inference with normal gaama prior distributions in regression problems. Bayesian Analysis, 5(1):171-188 2010.
Carlos M Carvalho, Nicholas G Polson, and James G Scott. Handling sparsity via the horseshoe. In Artificial Intelligence and Statislics, pp. 73 80. PMI.R, 2009.
16

---

Published in Transactions on Machine Learning Research (03/2025)
Ismaäl Castillo and Étienne Roquain. On spike and slab empirical Bayes multiple testing. The Annals of Stalistics, 48(5):2548 - 2574, 2020. doi: 10.1214/19-A051897 URL https://doi.org:: org/10. 1214/
19- A0S1898
S Grace Chang, Bin Yu, and Martin Vetterli. Adaptive wavelet thresholding for image denoising and compression. IEEE Transactions on Image Processing, 9(9):1532 1546, 2000.
Hugh A Chipman, Eric D Kolaczyk, and Robert E McCulloch. Adaptive bayesian wavelet shrinkage. Journal of the American Slalistical Associalion, 92(440):1413 11111111111111,, 1997.
Dorota M Dabrowska. Non-parametric regression with censored survival time data. Scandinavian Journal of Statistics, pp. 181-197, 1987.
Kushal K Dey, Dongyue Xie, and Matthew Stephens. A new sequence logo plot to highlight enrichment and depletion. BMC bioinformatics, 19:1 9, 2018.
David L Donoho and Iain M Johnstone. Ideal spatial adaptation by wavelet shrinkage. Biometrika, 81(3): 425 455, 1994.
Daniele Durante and Tommaso Rigon. Conditionally conjugate mean-field variational bayes for logistic models. Slatistical Science, 34(3):472 485, 2019.
Bradley Efron. Large-scale sinultaneous hypothesis testing: the choice of null hypothesis. Journal of the American Statistical Association, 99(465):96-11) 104, 2004.
James R Faulkner and Vladimir Minin. Locally adaptive smoothing with markov random fields and shrinkage priors. Bayesian analysis, 13(1):225, 2018.
Hong-Ye Ye Gao. Wavelet shrinkage estimates for heteroscedastic regression models. In MathSoft. Citeseer, 1997.
Alan Gelfand and Erin M Schliep. Spatial statistics and gaussian processes: beautiful marriage. Spatial Statistics, 18:8--104, 2016.
Edward I George and Robert E McCulloch. Variable selection via gibbs sampling. Journal of the American Statistical Association, 88(423)8881 889, 1993.
Wolfgang Härdle. Applied nonparametric regression. Number 19. Cambridge University Press, 1990.
Robert J Hodrick and Edward C Prescott. Postwar us business cycles: empirical investigation. Journal of Money, credit, and Banking, pp. 1-16, 1997.
Nikolaos Ignatiadis and Stefan Wager. Confidence intervals for nonparametric empirical bayes analysis.
Journal of the American Statistical Association, 117(539):11499 1166, 2022.
Tommi S Jaakkola and Michael I Jordan. A variational approach to bayesian logistic regression models and their extensions. In Sirth International Workshop on Artificial Intelligence and Statistics, pp. 283-294. PMLR, 1997.
Iain M. Johnstone and Bernard W. Silverman. Needles and straw in haystacks: Empirical Bayes estimates of possibly sparse sequences. The Annals of Slalistics, 32(4):1594 - 1649, 2004.
Seung-Jean Kim, Kwangmoo Koh, Stephen Boyd, and Dimitry Gorinevsky. l trend filtering. SIAM Review, (1(2):339-360 2009.
Youngseok Kim, Wei Wang, Peter Carbonetto, and Matthew Stephens. A flexible empirical bayes approach to multiple lincar regression and connections with penalized regression. Journal of Machine Learning Research, 25 25(185):1 59, 2024.
Roger Koenker and Ivan Mizera. Convex optimization, shape constraints, compound decisions, and empirical bayes rules. Journal of the American Statistionl Association, 1596500):674 674-688 2014.
17

---

Published in Transactions on Machine Learning Research (03/2025)
Gary Koop and Dimitris Korobilis. Bayesian dynamic variable selection in high dimensions. International Economic Review, 64(3):1047-1074, 2023.
Daniel R Kowal, David S Matteson, and David Ruppert. Dynamic shrinkage processes. Journal of the Royal Statistical Society Series B: Statistical Methodology, (4):78811804, 2019.
Yusha Liu, Peter Carbonetto, Michihiro Takahama, Adam Gruenbaum, Dongyue Xie, Nicolas Chevrier, and Matthew Stephens. A flexible model for correlated count data, with application to multicondittin differential expression analyses of single-cell RNA sequencing data. The Annals of Applied Statislics, 18(3): 2551 2575, 2024.
Florian Luisier, Cédric Vonesch, Thierry Blu, and Michael Unser. Fast interscale wavelet denoising of poisson-corrupted images. Signal Processing, 90(2):415 427, 2010.
Ryan Martin, Raymond Mess, and Stephen Walker. Empirical Bayes posterior concentration in sparse high-dimensional linear models. Bernoulli, 23(3):1822 1847, 2017.
Eric Moulines, François Roueff, and Murad S Taqqu. On the spectral density of the wavelet coefficients of long-memory time series with application to the logrregression estimation of the memory parameter. Journal of Time Series Analysis, 282):1155-187, 2007.
Yu-Chien Bo Ning and Ning Ning. Spike and slab bayesian sparse principal component analysis. Statistics and Computing, 34(3):118, 2024.
Trevor Park and George Casella. The bayesian lasso. Journal of the American Statistical Association, 103 (482):681-686, 2008.
Matias Quiroz, Robert Kohn, Mattias Villani, and Minh-Ngoc Tran. Speeding up meme by efficient data subsampling. Journal of the American Statistical Association, (14(556))831-848 2019.
Aaditya Ramdas and Ryan J Tibshirani. Fast and flexible admm algorithms for trend filtering Journal of Computational and Graphical Statistics, 25(3):839 858, 2016.
Andrés Ramirez-Hassan Hassan. Dynamic variable selection in dynamic logistic regression: application to internet subscription. Empirical Economics, 59(():909-998 2020.
Kolyan Ray and Botond Szabó. Variational bayes for high high-dimensional linear regression with sparse priors. Journal of the American Statistical Association, 117(539)) (27-399):270-12555 1281, 2022.
Veronika Rockova and Kenichiro McAlinn. Dynamic variable selection with spike-and-slab process priors. Bayesian Analysis, (((1):233-69)) 2021.
Edward A Roualdes. Bayesian trend filtering. arXiv preprint arXiv:1505.07711 2015.
Matthias Seeger and Guillaune Bouchard. Fast variational bayesian inference for non-conjugate matrix factorization models. In Artificial Intelligence and Statistics, pp. 0012-1018. PMLR, 2012.
Bernhard W Silverman. Some aspects of the spline smoothing approach to non-parametric regression curve fitting. Journal of the Royal Statistical Society: Series B (Methodological), 47(1):1 21, 1985.
Gabriele Steidl, Stephan Didas, and Julia Neumann. Splines in higher order tv regularization. International Journal of Computer Vision, 70(3):241 255, 2006.
Matthew Stephens. False discovery rates: new deal. Biostatistics, 18(2):275 294, 2017.
Ryan J Tibshirani. Adaptive piecewisc polynomial estimation via trend filtering. The Annals of Statistics, 42 (1):285 323, 2014.
Ryan J Tibshirani and Jonathan Taylor. The solution path of the generalized lasso. The Annals of Statistics, 33335:133-1 2011.
18

---

Published in Transactions on Machine Learning Research (03/2025)
Riaz A Usmani. Inversion of tridiagonal jacobi matrix. Linear Algebra and its Applications, 212(213): 413-414, 1994.
Gerrit JJ Van den Burg and Christopher KI Williams. An evaluation of change point detection algorithms. arXiv preprint ar orx::200030000 06222, 2020.
Stefan Van der Walt, Johannes I, Schönberger, Juan Nunez-Iglesias, Frangois Boulogne, Joshua D Warner, Neil Yager, Emmanuelle Gouillart, and Tony Yu. scikit-image: image processing python. PeerJ, 2:e453, 2014.
Gao Wang, Abhishek Sarkar, Peter Carbonetto, and Matthew Stephens. A simple new approach to variable selection in regression, with application to genetic fine mapping. Journal of the Royal Statislical Sociely Series B: Statistical Methodology, 82(5):1273- 10012277-330, 2020.
Wei Wang and Matthew Stephens. Empirical bayes matrix factorization. Journal of Machine Learning Research, 22:120- 1, 2021.
Jason Willwerscheid and Matthew Stephens. ebnm: An r package for solving the empirical bayes normal means problem using variety of prior families. arXiv preprint arXh::111000012,, 2021.
Haixu Wu, Jiehui Xu, Jianmin Wang, and Mingsheng Long, Autoformer: Decomposition transformers with auto-correlation for long-term series forecasting. Advances in Neural Information Processing Systems, 34: 22419 22430, 2021.
Dongyue Xie. Empirical Bayes Melhods for Coant Dala. PhD thesis, The University of Chicago, 2023.
Dongyue Xic. Local false sign rate and the role of prior covariance rank in multivariate empirical bayes multiple testing. arXiv preprint 161188,2.16618, 2025.
Dongyue Xie and Matthew Stephens. Discussion of "confidence intervals for nonparametric empirical bayes analysis". Journal of the American Statistical Association, 117((39):1186--1194 2022.
Zhengrong Xing, Peter Carbonetto, and Matthew Stephens. Flexible signal denoising via flexible empirical bayes shrinkage. Journal of Machine Learning Research, 22(93): 2021.
Yuxin Zou, Peter Carbonetto, Dongyue Xie, Gao Wang, and Matthew Stephens. Fast and flexible joint fine-mapping of multiple traits via the sum of single effects model. bioRain, pp. 2023-04, 2024.
19

---

Published in Transactions on Machine Learning Research (03/2025)
A
Additional plots
blocks
bumps
heavi
linear
gauss
step
Figure 8: Six signal functions in the simulation study. The blocks, bumps and Heavisine functions are originally proposed by Donoho & Johnstone (1994) for evaluating wavelet denoising method. The blocks and step functions are piecewise constant; bumps have most signals at 0 but jump at certain locations; linear is a piecewise linear function; Heavisine is piecewise twice differentiable function; Gauss is the Gaussian density function.
20

---

Published in Transactions on Machine Learning Research (03/2025)
EBTF
genlasso-tf
Gauss
wave-Bayes
wave-hard
susie-tf-10
susie-tf-20
susie-tf-30
GP
Step
BTF-DHS
BTF-NIG
HeaviSine
Blocks
Bumps
Linear
0.0
0.1
0.2
0.4
0.5
MAE
Figure 9: Simulation results: MAE scores of the competing methods across six signal functions.
0.30
MAE
0.25
X
0.20
0.15
log time)
BTF-DHS
genlasso-tf
susie-tf-30
x BTF-NIG
susie-tf-10
wave-Bayes
EBTF
susle-tf-20
wave-hard
GP
Figure 10: Simulation results: plots of log-transformed runtime MAE.
21

---

Published in Transactions on Machine Learning Research (03/2025)
EBTF-wavelet
EBTF-cvtf
Gauss
Step
HeaviSine
Blocks
Bumps
Linear
0.0
0.1
0.2
0.4
0.5
RMSE
Figure 11: Simulation results: Comparison of the two posterior mean initialization methods for EBTF (wavelet and cross-validated trend filtering). Both methods yield similar performance, but we prefer the wavelet initialization due to its faster computation.
22

---

Published in Transactions on Machine Learning Research (03/2025)
1.00
0.75
0.50
0.25
sin(angle)
0.00
-0.25
-0.50
-0.75
11.0
0
100
200
300
400
500
600
position
Honey bee movement states. The 2-axis is the position and y-axis is sine of the head angle of a single bee. Black line is the EBTF fitted signal, and blue line is the genlasso-tf estimated signal.
Figure 12:
17.5
15.0
12.5
CO2 emissions
10.0
7.5
5.0
2.5
0.0
100
150
200
Year since 1800
Figure 13: CO2 emissions per person in Canada. Black line is the EBTF fitted signal, and blue line the genlasso-tf estimated signal.
23

---

Published in Transactions on Machine Learning Research (03/2025)
3000
2500
Home runs
2000
1500
1000
500
0
20
40
60
80
100
120
Year since 1900
Figure 14: Number of home runs in the American League of baseball since 1900. Black line is the EBTF fitted signal, and blue line is the genlasso-tf estimated signal.
1400
Volume of Nile River at Aswan
1200
1000
600
20
40
60
100
Year since 1870
Figure 15: The volume of the Nile river at Aswan over each year. There clear change point in 1898 due to built dam.
24

---

Published in Transactions on Machine Learning Research (03/2025)
1e6
1.08
Total kilometers of rail lines the world
1.06
1.04
1.02
1.00
0.98
0.96
0.94
10
15
Year since 1980
Figure 16: Total length rail lines in the world, in kilometers.
leg
5
Number of robocalls in the US
4
2
Month since 2015-05
Figure 17: The number of robot calls in US.
25

---

Published in Transactions on Machine Learning Research (03/2025)
250000
Number of applicants
200000
150000
100000
50000
25
50
75
100
125
150
175
200
Year since 2002
Figure 18: The number license plate applications in Shanghai since 2002. Two outliers were removed.
140000
135000
Nuclear magnetic response
130000
125000
120000
115000
110000
105000
100
200
300
400
500
600
Time
Figure 19: Well-log data. captures the nuclear magnetic responses over time. The length of the series has been reduced by sampling every 6 time points. Outliers were removed.
26

---

Published in Transactions on Machine Learning Research (03/2025)
B Empirical Bayes Gaussian Variance
We give details the general Gaussian variance problem. Consider the following model on Gaussian variance: for i = 1, 2, .... n4
2y/u ~ N(0, orw;), Wi ~ g()),
(14)
where o² is known. An empirical Bayes Gaussian variance procedure returns by maximizing the marginal log likelihood Σg logp(z;), and calculates the posterior lw; (-) (r, 9). The objective function of EBGV problem is
Σ Bgg log g(w:) gm, (w:)
FEBGV Ey log (ax 0, o
The procedure defines mapping from observations to the estimated prior and posterior distribution, and is denoted as
(0. fw) EBGV(2, o
In the EBGV problem equation 14 if wy has discrete prior
Wy Discrete of T),
T);
then the marginal distribution of x, is
p(x;)= = Ep(z)/wpp(u Σ [
The posterior distribution of ws
p(we(x,) II
[w:=oo]) i
where
And we have
B(w (
Σ 0ak
The marginal and posterior distribution are given by their definitions, and the Bayes formula. The expectation of 1/w; follows directly from the definition of discrete random variables.
C Algorithms
27

---

Published in Transactions on Machine Learning Research (03/2025)
Algorithm VEB algorithm for fitting EBTF
(outline only)
Input: Data yi, variances 82, for = 1,2,.... n.
Init: Posterior mean 8, posterior precision matrix diagonal d, and super-diagonal e, residual variance o². repeat
1. Update Ø(-) and qw by solving the EBGV problem equation 16
2. Update qa by updating the posterior precision matrix (its diagonal and super-diagonal elements only); then B by solving (tridiagonal) banded linear system;
3. Update residual variance o².
until converged
Derivations
D.1 Derivation of VEB updates for fitting model equation 2
Based on the model equation 2 and the posterior distribution equation Ö the evidence lower bound can be writlen in veclor-matrix form in B. as
n FEBTF = log 2TO
log $
2
2 1=1
T
2ny S
3 TS 1-131-65-0000
5-1 tr(S-
20
1
8 D D33
tr tr(D" W DV) + log V
2x
20
n-1
1
K
ayk (log & log 2To ox log azx);
2
where S = diag(sf) is the known diagonal variance matrix, D is the first order difference matrix equation ! and W- diag(W;;) is diagonal weight matrix, with Wyj Σ k= a
The update of cach parameter (denoted generally as 0) is given by solving the root-finding equation 0F/00 = 0.
D.2 Derivation of VEB updates for maximizing the ELBO FMNV
The ELBO for model equation
with posterior distribution being equation is
FNN E, log pgy, B, W) E, log (B, W),
Eg log p(B W)
Eg log g(W) qw (W)*
(15)
Eg log p()
(B)
We now show that the variational inference update for qs given fwr
-S-1 D' y.
V = D= D'
where
W-1 diag( w 1 )
and given qa, the update on g(-). dw is obtained by solving EBGV problem
((V3 ub .0?
(w)) = FBGV
(16)
28

---

Published in Transactions on Machine Learning Research (03/2025)
where
by = (D))s,
bby (DVD")>)
Given ys, the ELBO for updating tw is Fa
Eyg. (Egg (log (())()))))))) + Σ Beg log g(w,)
Qre (wy)
log Elg N( NV B3 + uby; 0, (wy) + Σ log g(wz) Ow (wy)
which is exactly the objective function for EBGV problem. And we have
T (. 103+ Uby 00y;0,07 0,02 Ex Te Nx N( 03 1 + vn, 0,
0sk =
Given Iw: the ELBO related to qa is
FMNV (ga) Elgg + EE (log (lgg)())))) (8)
The update formulas for B, are given by solving the root-finding equation 0F = 0, and (a 0.
Since the marginal distribution in EBGV problem is the Gaussian mixture distribution, the update of prior paramcters (aî) are the same as the ones in VEB updates for fitting model equation 2 with Φjk azk- It is obvious the update for o² is also the same since the ELBO related to o² in FMNV is the same as FFBTF(O²). Furthermore, the variational inference algorithm is the same as Algorithm
D.3 Derivation of VEB updates for maximizing the ELBO FML
The ELBO FMLR is
FALR = BB, ply, b, z) E.g log ((,,)) 2)
n 2
log 0²
2 Σ log s
1 (y'S S
1-1 y
-2y" Hb + -- -1- t
20
1
+
Oej (log
(83 + (b,+), log 0² 000 20
2
J.K
1
+ log [Vgl
Σ
ajı loga agk
2
j.k
"
log
1
Σ
log s
2
(yT
T
S y 2y 1T -i0Ts--H++- t 5-
20
1 W'b + (VVV)) + log |/
20
agggloo T
log 2o log 0.kk),
29

---

Published in Transactions on Machine Learning Research (03/2025)
where W diag(w,) and wy = Ex ayt/o
Let
b
B1 b
088 Vo.g
V3L.b
Vg
Vo
Given (b) N(b:b,Vb) Vi),, the inducing posterior distribution on / is also multivariate normal, denoted
(a(+) = (N: 3, Va); where B = Hb and Vi3 = HV, Based on the relationship between B, b, b, we also have
DB,
V. = DVgD"
Since H. is triangular matrix with diagonal elements all being 1, have [M] 1, 1, and VVl |H = =
Given the above equivalence, the ELBO can be written as
FMLR
11
log 2TO²²
log 2
2
2 1=1
TS 1-S-111 y
2n S 3+ (V3)
20
1 2
8" Dd W DB
tr(D"
+ log Vi
2
20
n-1 1 K
ast n.
log 2no ox log axx),
2
which is exactly the same as FEBTF-
D.4 Update formulas for sparse EBTF
Let Wo = diag(oso/o8), the updates for empirical Bayes sparse trend filtering are as follows:
1. The update for posterior probabilities are
1 log 2
(DB)² + (DVD")a )
log ask x log T&
2
20202
83 + Va
a,p log To
log 2 of
202
logan x log(1
no);
then
0.10
01
011 f-
a10
+
an
010 +m
and
Olik
Oik f-
Σ, Osl
for i 2......, and kt 1...,
The update for the posterior covariance matrix V and the posterior mean are
= 02 + D WD+ -1
Ø = Vu/o².
30

---

Published in Transactions on Machine Learning Research (03/2025)
3. The update for prior variance of for k 1,2, K and prior probabilities are o Ej ((D))) (DVD))))))
0² E Ask
n
0
Oik: for k 1,.., K,
T
for k = 0,1,.... K.
11=
4. Update 0²
Q S- + +D* Wo.
0² y" S 2y 5- + 2T + tr(2V)
2n + 0.00
The updates of each parameter (denoted generally as 0) is obtained by solving the root-finding equation
OF 00 = 0,
where is the ELBO for EBSTF, defined as
F log B, (g
2 log 2TTO² Σ log s/ 2
5 2y S 5 "0
20
+ *111(1o((1 - mo))= + Σ
aooglog Elog N (3,; 0, 02
1
rr DJ WD3
1
tr
log [
20
202
Σ
(logo n
log 2 0
loo axx) - a11 log Q1
2
2m2kkml
31


Reviews:
[{"review_id": "DLniI0ZZE8", "reviewer": "Reviewer_JeQM", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\n- The paper is generally well-written and easy to follow.\n- The proposed variational Bayes algorithm is more scalable than baseline methods, at least on synthetic data.\n\nWeaknesses:\n- The novelty appears limited. The probabilistic formulation of the trend filtering problem is well-established in the literature, and deriving a coordinate ascent algorithm for variational posteriors is relatively straightforward, primarily involving standard calculations rather than substantial technical innovations.\n- The real-data experiments lack quantitative evaluation metrics, making it difficult to assess the extent of improvement over baseline methods.\n- While the paper introduces multiple variants of the proposed model and discusses their trade-offs, the explanations remain somewhat vague. It is unclear from both the text and the experiments which variants are preferable under specific circumstances, raising questions about their necessity.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces a Bayesian trend filtering algorithm based on variational inference. Building on an existing probabilistic model, the authors formulate variational distributions and derive a coordinate ascent algorithm for inference. They further explore alternative parameterizations that yield the same variational objectives, arguing that these alternatives enhance modeling flexibility while preserving analytical tractability. The proposed method is evaluated on both synthetic and real-world datasets, demonstrating improved speed and accuracy compared to baseline approaches.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- The paper primarily derives a coordinate ascent variational inference algorithm, but beyond standard calculations such as writing down the ELBO, taking derivatives, and solving for updates, it is unclear what the core technical contribution is. Please clarify whether there are additional methodological innovations that go beyond these routine steps, as this is crucial for assessing the novelty of the work.\n- The necessity of the three model variants (the primary model, MNV, and MLR) is not entirely clear. While the paper states that Section 5 is intended to illustrate their differences and motivations, I found it difficult to grasp the key takeaways from that discussion. Could you further elaborate on why these variants are needed, in what scenarios they offer advantages, and how they compare empirically?", "broader_impact_concerns": "I don't see any potential broader impact concerns from this work.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "XodGjC8Dd1", "reviewer": "Reviewer_tkuc", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\n- The paper is generally well-written, and the theoretical results are sound.\n- On the datasets evaluated, the proposed method is indeed fast and accurate relative to the baselines. The empirical results suggest that EBTF underfits less than competing methods.\n- The sparse EBTF extension is a nice way to deal with periodically dropped signals, and shows the flexibility of the dynamic linear model formulation of trend filtering.\n\n## Weaknesses\n\n**Regarding Audience**\nI'm somewhat skeptical of the audience; personally, I don't feel that trend filtering is a common interest in the machine learning community, and it is more typically discussed in statistics and signal processing communities. While variational inference may be more generally of interest, the coordinate ascent algorithm derived in this work is a relatively straightforward application of the VEB framework.\n\n**`susie-tf-20` and Figure 1**\nI find Figure 1 a bit misleading. `susie-tf-20` has an outlier in performance on one dataset, which brings its average RMSE up dramatically. Aside from datasets with pathologically \"spiky\" data, `susie-tf` seems quite strong.  Plus, averaging RMSEs like this (ignoring the scale of different tasks) is not ideal in the first place.\n\n**Code Availability**\nThe paper mentions that \"all the code and analysis are available in the package `ebtfPy` on GitHub.\" I did not search for this package as I was concerned it might break the authors' anonymity, and I did not see any anonymized version provided, so I could not look at the code. I would be interested to look at the code, if an anonymized version can be provided.\n\n**On the Use of MCMC**\nIt's mentioned that \"while it's possible to use MCMC for sampling from the posterior, it's intractable for large-scale datasets.\" This is true, of course, but I don't think it is a major concern for datasets of the size considered in this paper ($n = 1024$ in the experiments), at least not to the level of \"intractability.\" For example, the listed reference (Quiroz et al., 2019) uses datasets several orders of magnitude larger in their experiments. I would thus be interested to see a comparison in run times between MCMC and VB.\n\n**On the Extent of Experiments**\nI find the experiments somewhat lacking. As an example, the datasets with quantitative results are somewhat simplistic. and finding more realistic datasets where qualitative performance may be evaluated would significantly strengthen the results in my opinion. As another example, using ash priors and variational Bayes introduces two moving parts that are not properly ablated: how much performance in terms of accuracy is lost by resorting to variational inference, and how much slower is MCMC? There is also mention that the theory in the paper holds for higher-order trend filtering problems, but there are no experiments. Finally, ash priors are used somewhat arbitrarily; it would be nice to investigate the effect of different shrinkage priors in trend filtering.\n\n## Questions\n- In the simulations, how is $\\hat{\\beta}$ chosen for EBTF? As the mean of the variational posterior?\n- The initializations in Remark 2.2 seem reasonable, but how much do they matter? For example, if $\\beta_i$ is intiialized as $y_i$, are results significantly worse? What if they are initialized using the results of `susie-tf-20` or `genlasso-tf`?\n- Kim et al. (2024) propose a VEB approach to multiple linear regression. Given the equivalence presented in Section 3.2, could the authors further comment on the relevance of Kim et al. (2024) to EBTF? Are the main algorithmic differences using banded solvers?\n- Why are $L=10$ and $L=20$ chosen for `susie-tf`?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces a variational empirical Bayes (VEB) approach to trend filtering. In particular, a variational solution to a first-order trend filtering problem is presented, using shrinkage priors based on mixtures of normals. The problem is originally formulated in terms of a dynamic linear model and a coordinate ascent algorithm presented; it is later shown that various other formulations of trend filtering lead to the same parameter updates. Experiments on toy datasets show strong quantitative results with faster runtime than similarly performing methods.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "No", "audience_explanation": "", "requested_changes": "## Critical\n- Please address the above concerns regarding the experiments.\n- Nearly all citations should be using `\\citep`, but are not.\n- Please include standard errors for all results in Table 1, even if the general scale is similar to the first three datasets.\n- Several papers cited reference the arXiv version of a paper, even though a more recent journal or conference version is available (e.g., Ning & Ning (2024), Castillo & Roquain (2020), and Kim et al. (2024)). If there is no particular reason to cite the arXiv version, I suggest the authors consider citing the latest version. Note that I did not check all references for this, and there may be more examples.\n\n## Minor\n- There are many places in the bibliography with missing capitalization, either in proper nouns or in journal/publisher titles; for example, \"Pareto\", \"Cambridge University Press\", and \"Signal Processing\". \"\\ell_1\" should also use math mode.\n- At the beginning of 2.3.2, it should read \"for the prior *in equation* 3 [...]\".\n- The acronym \"MLR\" (multiple linear regression) is not introduced.\n- The trend $\\beta$ is not introduced in Section 2.1.\n\n# References\n\nKim, Y., Wang, W., Carbonetto, P., & Stephens, M. (2024). A flexible empirical Bayes approach to multiple linear regression and connections with penalized regression. _Journal of Machine Learning Research_, _25_(185), 1-59.\n\nNing, Y. C. B., & Ning, N. (2024). Spike and slab Bayesian sparse principal component analysis. _Statistics and Computing_, _34_(3), 118.\n\nCarvalho, C. M., Polson, N. G., & Scott, J. G. (2009). Handling sparsity via the horseshoe. In _Artificial Intelligence and Statistics_ (pp. 73-80). PMLR.\n\nCastillo, I., & Roquain, É. (2020). On Spike and Slab Empirical Bayes Multiple Testing. _The Annals of Statistics_, _48_(5), 2548-2574.\n\nJaakkola, T. S., & Jordan, M. I. (1997, January). A variational approach to Bayesian logistic regression models and their extensions. In _Sixth International Workshop on Artificial Intelligence and Statistics_ (pp. 283-294). PMLR.\n\nSeeger, M., & Bouchard, G. (2012). Fast variational Bayesian inference for non-conjugate matrix factorization models. In _Artificial Intelligence and Statistics_ (pp. 1012-1018). PMLR.\n\nQuiroz, M., Kohn, R., Villani, M., & Tran, M. N. (2019). Speeding up MCMC by efficient data subsampling. _Journal of the American Statistical Association_.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "her8ehOYtH", "reviewer": "Reviewer_Qhzt", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper is clearly written. The only slight oddness to me in the presentation was the sparse version being described at the end after the results for the non-sparse version. It would be cleaner to have the sparse EBTF described along with the main approach and then added as a comparison method in Table 1 - presumably it should do well on the less smooth simulated data. As an additional point here: the sparse EBTF doesn't really feel like a different model, since as far as I can see it is equivalent to fixing one of the variances in EBTF to a small value. \n\nThe modeling and inference is all very reasonable. The simulated results are also fine. The \"real data\" feels very limited: this is a small dataset and there is no quantitative assessment, which seems relatively straightforward in this context: see \"Requested Changes\". \n\nIs there a reason smoothing approaches like Gaussian process regression aren't compared to? They would be a poor fit for data with jumps (although there are GP change point models) but would feel natural for something as smooth as the motorcycle data. \n\nMinor comments: \n\nThe authors say in the intro they will denote the vector beta as (beta_i) and W as diag(w_i) which is poor notation imo. i refers to one element. If you want to specify all the individual elements simply use beta, or specify {beta_i:i=1...n} or similar. Fortunately they don't actually seem to use this notation despite introducing it. \n\nWhat does the cdot between D^(1) and D^(k) denote? I believe this is matrix multiplication, in which case there should be no symbol there: matrix multiplication is implied. If it is Hadamard that should be specified. \n\nI'm not aware of people calling what is being done here \"variational empirical Bayes\" (VEB), I just call this variational Bayesian EM (VBEM), where the M steps are on the optimized parameters. I don't hate VEB though. \n\nsigma and pi are parameters, _not_ variational parameters (page 4)\n\nCalling the updates a \"Theorem\" is overselling. A \"computation\" maybe? \n\nTable 1 would be easier to grok as a barchart, and R^2 is easier to interpret than RMSE.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors consider a Bayesian approach to \"trend filtering\" which is smoothing of a one-dimensional trajectory. They construct a model with an ash-esque prior on the signal change between time-points, and develop a coordinate ascent variational inference procedure for fitting this model. They show promising performance, both computationally and statistically, on simulated data. They demonstrate the method on a small \"real world\" dataset.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "For the real data hold out some time points from the model and ask how well they can be imputed. I would like to see at least two more datasets with different characteristics - the motorcycle data looks pretty smooth, so having data with substantial jumps would be a nice addition. Ideally additionally having some larger datasets (financial time series maybe?) would be valuable. Especially since the VI is motivated by scalability it is disappointing to only see an application to a tiny dataset.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: See comments above.

One requirement I will impose prior to publication is to turn all theorem/lemma/etc statements into just regular mathematics interspersed with text both in the main text and appendix (appropriately formatted without an environment). These are all overclaimed, as they are all just routine derivations.
Decision: Accept (poster)