# An Alternating Direction Method of Multipliers Algorithm for the Weighted
Fused LASSO
Signal Approximator

Louis Dijkstra  Leibniz Institute for Prevention Research & Epidemiology –
BIPS, Bremen, Germany  Moritz Hanke  Leibniz Institute for Prevention Research
& Epidemiology – BIPS, Bremen, Germany  Niklas Koenen  Leibniz Institute for
Prevention Research & Epidemiology – BIPS, Bremen, Germany  Faculty of
Mathematics & Computer Science, University of Bremen, Bremen, Germany  Ronja
Foraita111Corresponding author. E-mail<EMAIL_ADDRESS>Leibniz
Institute for Prevention Research & Epidemiology – BIPS, Bremen, Germany

###### Abstract

We present an Alternating Direction Method of Multipliers (ADMM) algorithm
designed to solve the Weighted Generalized Fused LASSO Signal Approximator
(wFLSA). First, we show that wFLSAs can always be reformulated as a
Generalized LASSO problem. With the availability of algorithms tailored to the
Generalized LASSO, the issue appears to be, in principle, resolved. However,
the computational complexity of these algorithms is high, with a time
complexity of $\mathcal{O}(p^{4})$ for a single iteration, where $p$
represents the number of coefficients. To overcome this limitation, we propose
an ADMM algorithm specifically tailored for wFLSA-equivalent problems,
significantly reducing the complexity to $\mathcal{O}(p^{2})$. Our algorithm
is publicly accessible through the R package wflsa.
Keywords: ADMM algorithm $\cdot$ Generalized LASSO $\cdot$ Image denoising
$\cdot$ LASSO $\cdot$ Weighted Fused LASSO Signal Approximator

##  1 Introduction

Several penalized regression methods have been proposed to address
identifiability challenges and enforce sparsity in high-dimensional settings,
where the number of variables ($p$) exceeds the number of observations ($n$).
These methods impose penalties on regression coefficients, aiming to regulate
model complexity and promote more robust and interpretable results. Perhaps
the most prominent among these methods is the Least Absolute Shrinkage and
Selection Operator (LASSO; Tibshirani, (1996)), which is commonly defined as

$\widehat{\bm{\beta}}=\operatorname*{arg\,min}_{\bm{\beta}\in\mathds{R}^{p}}\frac{1}{2}\left\lVert\mathbf{y}-\mathbf{X}\bm{\beta}\right\rVert_{2}^{2}+\lambda_{1}\left\lVert\bm{\beta}\right\rVert_{1}.\qquad\text{{(LASSO)}}$
(1)

Here, $\bm{\beta}$ represents a $p$-dimensional vector of coefficients,
$\mathbf{y}$ stands for the $n$-dimensional response vector, $\mathbf{X}$ is a
$(n\times p)$-dimensional matrix of observations, and
$\left\lVert\cdot\right\rVert_{1}$ and $\left\lVert\cdot\right\rVert_{2}$
represent the $\ell_{1}$\- and $\ell_{2}$-norm, respectively. The term
$\lambda_{1}\left\lVert\bm{\beta}\right\rVert_{1}$ induces sparsity by
‘shrinking’ some of the coefficients to zero. The penalty term
$\lambda_{1}\geq 0$ regulates the sparseness of the solution, i.e., higher
values of $\lambda_{1}$ lead to sparser estimated regression vectors. We
assume throughout that the mean of $\mathbf{y}$ is zero.

A plethora of variations on the LASSO exist, one of which is the _Fused LASSO_
, introduced by Tibshirani et al., (2005). The Fused LASSO assumes that
consecutive regression coefficients (i.e., $\beta_{i}$ and $\beta_{i+1}$)
should be more similar than those that are ‘further apart’. It takes advantage
of this assumption by penalizing the absolute differences between consecutive
regression coefficients which encourages piecewise constant or piecewise
linear solutions. Formally, the Fused LASSO is defined as

$\widehat{\bm{\beta}}=\operatorname*{arg\,min}_{\bm{\beta}\in\mathds{R}^{p}}\frac{1}{2}\left\lVert\mathbf{y}-\mathbf{X}\bm{\beta}\right\rVert_{2}^{2}+\lambda_{1}\left\lVert\bm{\beta}\right\rVert_{1}+\lambda_{2}\sum_{i=1}^{p-1}|\beta_{i+1}-\beta_{i}|.\qquad\text{{(Fused
LASSO)}}$ (2)

The additional penalty term is referred to as the smoothness penalty, as it
fosters ‘smoother’ transitions from one coefficient to the next. The
associated tuning parameter, $\lambda_{2}\geq 0$, governs the tolerance for
differences; higher values of $\lambda_{2}$ lead to smaller differences
between consecutive coefficients.

The LASSO and Fused LASSO are commonly used in regression analysis to model
the relationship between a dependent variable ($\mathbf{y}$) and a set of
independent variables ($\mathbf{X}$). However, beyond their conventional use
in regression, these penalized techniques are also applied to transform the
original observed vector ($\mathbf{y}$) in cases where sparsity (various
values being zero) and smoothness (small differences between certain values)
are expected. Within this framework, the resulting transformed vector is
represented by the coefficient vector $\bm{\beta}$. Later in the text, we
discuss some applications in which such a transformation might be useful.

One method to achieve such a transformation is a particular case of the Fused
LASSO, where the design matrix $\mathbf{X}$ is equal to the identity matrix
$\mathbf{I}$. This specific instance is referred to as the _1-dimensional
Fused LASSO Signal Approximator_ (1-dimensional FLSA; Kim et al., (2009);
Condat, (2013)):

$\widehat{\bm{\beta}}=\underset{\bm{\beta}\in\mathbb{R}^{p}}{\text{argmin}}\frac{1}{2}\lVert\mathbf{y}-\bm{\beta}\rVert_{2}^{2}+\lambda_{1}\lVert\bm{\beta}\rVert_{1}+\lambda_{2}\sum_{i=1}^{p-1}|\beta_{i+1}-\beta_{i}|.\quad\text{{(1-dimensional
FLSA)}}$ (3)

Here, the vector $\bm{\beta}\in\mathds{R}^{p}$ represents the transformed
$\mathbf{y}\in\mathds{R}^{p}$. Similar to the regression scenario, the first
penalty term induces sparsity by driving certain coefficients to zero, while
the second term encourages similarity between consecutive coefficients.

The 1-dimensional FLSA finds applications in various domains such as
comparative genomic hybridization (Tibshirani and Wang,, 2008) and chromosomal
microarray analysis (Hoefling,, 2010). Higher dimensional versions of the FLSA
exist as well. The 2-dimensional version is, for example, used to denoise
images, where neighboring pixels are thought to be similar (Xu and Noo,,
2017), see Section 6. Another important application is change point detection,
which are abrupt changes of the values in the transformed regression
coefficient vector. This can be useful in a variety of fields in which
systematic changes in the system of interest need to be detected (Rojas and
Wahlberg,, 2014; Harris et al.,, 2022). Section 6 provides an example of how
the wFLSA can be applied to heterogeneous image smoothing and compares it to
other algorithms.

Hoefling, (2010) extended the definition of the Fused LASSO by introducing an
undirected graph $G=(V,E)$, where the vertices $V=\\{1,2,\ldots,p\\}$
correspond to the coefficients $\beta_{1},\beta_{2},\ldots,\beta_{p}$. The
adjacency of a pair of coefficients is determined by the graph’s edge set
$E\subseteq V\times V$, where $\\{i,j\\}\in E$ denotes an edge between
vertices $i$ and $j$ in the graph. Based on the graph, one can generalize the
FLSA as follows:

$\widehat{\bm{\beta}}=\operatorname*{arg\,min}_{\bm{\beta}\in\mathds{R}^{p}}\frac{1}{2}\left\lVert\mathbf{y}-\bm{\beta}\right\rVert_{2}^{2}+\lambda_{1}\left\lVert\bm{\beta}\right\rVert_{1}+\lambda_{2}\sum_{\\{i,j\\}\in
E}|\beta_{i}-\beta_{j}|.\qquad\text{{(Generalized FLSA)}}$ (4)

The Generalized FLSA encourages similarity between coefficients linked by an
edge in the graph. Let $\mathbf{A}$ denote the adjacency matrix of the graph
$G$. The 1-dimensional FLSA, see eq. (3), is a specific instance of the
Generalized FLSA, as it is equivalent when the adjacency matrix is defined as

$\mathbf{A}_{\text{1-dimensional}}=\begin{bmatrix}0&1&0&\dots&0&0\\\
1&0&1&\dots&0&0\\\ 0&1&0&\dots&0&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&0&\dots&0&1\\\
0&0&0&\dots&1&0\\\ \end{bmatrix}.$

Several efficient algorithms have been developed for solving the Generalized
FLSA (Hoefling,, 2010; Qian and Jia,, 2016; Son and Lim,, 2019).

In this context, we consider an even more general version: the _Weighted
Generalized Fused LASSO Signal Approximator_ (wFLSA) in which the edges in the
graph are weighted. Consider the graph $G=(V,E,w)$, where $w:V\times
V\rightarrow[0,\infty)$ is a weight function assigning non-negative values to
each possible edge $\\{i,j\\}$. Let the weighted adjacency matrix of the graph
be denoted as $\mathbf{W}=\\{w_{ij}\\}$, where $w_{ij}=w(\\{i,j\\})$. The
wFLSA is defined as

$\widehat{\bm{\beta}}=\operatorname*{arg\,min}_{\bm{\beta}\in\mathds{R}^{p}}\frac{1}{2}\left\lVert\mathbf{y}-\bm{\beta}\right\rVert_{2}^{2}+\lambda_{1}\left\lVert\bm{\beta}\right\rVert_{1}+\lambda_{2}\sum_{i<j}w_{ij}|\beta_{i}-\beta_{j}|.\quad\text{{(Weighted
Generalized FLSA)}}$ (5)

Not only does this approach allow for determining which coefficients should be
similar, but it also provides the flexibility to regulate the degree of
similarity by assigning varying weights to the edges. It is evident that the
1-dimensional and Generalized FLSA are special cases of this weighted version.

To the best of our knowledge, no algorithm for the wFLSA exists in the
literature. Our approach to estimating a wFLSA follows two main steps. First,
we establish that the wFLSA can be reformulated as a special case of the
Generalized LASSO framework (Tibshirani and Taylor,, 2011), see Section 2.
Throughout, we refer to Generalized LASSO problems of this kind as _wFLSA-
equivalent problems_. By demonstrating their equivalence, we essentially
resolve the problem at hand. Existing path algorithms for solving the
Generalized LASSO, see Arnold and Tibshirani, (2016) and Zhao and Bondell,
(2020), theoretically possess the capability to address problems equivalent to
wFLSAs. However, the computational complexity associated with these path
algorithms becomes prohibitively high for wFLSA-equivalent problems, i.e., in
the order of $\mathcal{O}(N_{\text{GL}}\cdot p^{4})$, where $N_{\text{GL}}$ is
the number of iterations needed for the algorithm to converge (Arnold and
Tibshirani,, 2016; Zhao and Bondell,, 2020). This makes it ill-suited for
situations with high-dimensional data characterized by a large number of
coefficients.

To overcome this computational limitation, we proceed with the second step, in
which we introduce a tailored Alternating Direction Method of Multipliers
(ADMM; Boyd and Vandenberghe, (2004)) algorithm designed specifically to
tackle wFLSA-equivalent problems, see Section 3. The ADMM algorithm, an
optimization method used in a variety of fields (Forero et al.,, 2010;
Wahlberg et al.,, 2012; Yang et al.,, 2018; Glowinski et al.,, 2022), divides
the original optimization problem into three subproblems (Boyd and
Vandenberghe,, 2004). At each iteration, the ADMM updates 1) the primal
variables, 2) the dual variables, and 3) a penalty parameter (hence its name
‘alternating direction’). The primal variables correspond to the variables of
interest in the original optimization problem. Updating them involves
minimizing an augmented Lagrangian function, consisting of components of the
original objective function, alongside Lagrange multipliers terms to enforce
the constraints. Subsequently, the algorithm updates the dual variables based
on the differences between the current value of the primal variables and the
constraints. The ADMM’s penalty parameter regulates the trade-off between
staying close to the original problem on the one hand and adhering to the
constraints on the other. The ADMM proves particularly advantageous for
handling large-scale data, as it can be parallelized (Boyd and Vandenberghe,,
2004).

The proposed ADMM algorithm notably reduces the complexity to
$\mathcal{O}(N_{\text{ADMM}}\cdot p^{2})$, where $N_{\text{ADMM}}$ is the
number of iterations for the algorithm to converge (see Section 4), rendering
it more scalable to real-world problems.

The paper is structured as follows. In Section 2, we demonstrate how the wFLSA
can be represented as a special instance of the Generalized LASSO. We continue
by describing the ADMM algorithm tailored to these wFLSA-equivalent
Generalized LASSO problems in Section 3. We show in Section 4 that, by
leveraging the unique structure, the computational complexity of the ADMM
algorithm is $\mathcal{O}(N_{\text{ADMM}}\cdot p^{2})$. Section 5 discusses
various methods for choosing the weight matrix and the two tuning parameters
based on data. Section 6 provides an example of how the wFLSA can be used for
heterogeneous image smoothing. Finally, in Section 7, a conclusion and
suggestions for future research are provided. Appendix A includes an example
of how the wFLSA can be formulated as a Generalized LASSO problem when $p=4$.
In addition, Appendix B provides information on how to choose matrix $Q$ used
in the ADMM efficiently. An implementation of the algorithm is publicly
available as an R package under https://github.com/bips-hb/wflsa. The code for
the example in Section 6 can be found at https://github.com/bips-hb/wFLSA-
paper.

##  2 The Weighted Fused LASSO Signal Approximator as a Special Case of the
Generalized LASSO

First, we use a fundamental result in the work of Friedman et al., (2007), in
which they introduce a pathwise coordinate optimization algorithm used to
solve the 1-dimensional FLSA, among others. Let us denote the estimate for the
regression coefficients of the wFLSA for specific tuning parameters and a
given weight matrix as

$\begin{split}\widehat{\bm{\beta}}(\lambda_{1},\lambda_{2},\mathbf{W})&=\left[\widehat{\beta}_{1}(\lambda_{1},\lambda_{2},\mathbf{W}),\widehat{\beta}_{2}(\lambda_{1},\lambda_{2},\mathbf{W}),\ldots,\widehat{\beta}_{p}(\lambda_{1},\lambda_{2},\mathbf{W})\right]^{\top}\\\
&=\arg\min_{\bm{\beta}\in\mathbb{R}^{p}}\frac{1}{2}\|\mathbf{y}-\bm{\beta}\|_{2}^{2}+\lambda_{1}\|\bm{\beta}\|_{1}+\lambda_{2}\sum_{i<j}w_{ij}|\beta_{i}-\beta_{j}|.\end{split}$

Friedman et al., (2007) demonstrated that once an estimate
$\widehat{\bm{\beta}}(0,\lambda_{2},\mathbf{W})$ is obtained (the wFLSA
without the sparsity penalty, i.e., $\lambda_{1}=0$), one can readily derive
the estimate $\widehat{\bm{\beta}}(\lambda_{1},\lambda_{2},\mathbf{W})$ for
any $\lambda_{1}>0$. This is achieved by applying a soft thresholding function
to each of the values of $\widehat{\bm{\beta}}(0,\lambda_{2},\mathbf{W})$. The
soft thresholding function is defined as

$\text{ST}_{\lambda_{1}}(x):=\begin{cases}x-\lambda_{1}&\text{if
}x>\lambda_{1},\\\ 0&\text{if }-\lambda_{1}\leq x\leq\lambda_{1},\\\
x+\lambda_{1}&\text{otherwise}.\end{cases}$ (6)

The estimate $\widehat{\bm{\beta}}(\lambda_{1},\lambda_{2},\mathbf{W})$ is
then given by

$\widehat{\bm{\beta}}(\lambda_{1},\lambda_{2},\mathbf{W})=\left[\text{ST}_{\lambda_{1}}\left(\widehat{\beta}_{1}(0,\lambda_{2},\mathbf{W})\right),\text{ST}_{\lambda_{1}}\left(\widehat{\beta}_{2}(0,\lambda_{2},\mathbf{W})\right),\ldots,\text{ST}_{\lambda_{1}}\left(\widehat{\beta}_{p}(0,\lambda_{2},\mathbf{W})\right)\right]^{\top}.$

This result is important for three reasons. First, it allows us to focus
solely on solving the simpler problem
$\widehat{\bm{\beta}}(0,\lambda_{2},\mathbf{W})$, which we will do throughout
the rest of the paper. Second, determining the estimate for a different
$\lambda_{1}>0$ is computationally efficient, as applying soft thresholding to
a $p$-dimensional vector requires $\mathcal{O}(p)$ steps. Third, and perhaps
most importantly, we can identify the specific $\lambda_{1}$ values at which
the $i$-th coefficient enters the active set, i.e., becomes non-zero. This
stems from the fact that
$\text{ST}_{\lambda_{1}}(\widehat{\beta}_{i}(0,\lambda_{2},\mathbf{W}))$ only
becomes non-zero if
$|\widehat{\beta}_{i}(0,\lambda_{2},\mathbf{W})|>\lambda_{1}$. Consequently,
for a fixed value of $\lambda_{2}$, we know the entire _solution path_ of
$\lambda_{1}$. This greatly simplifies the problem at hand and substantially
reduces the computation time required.

Now, we show that the problem of estimating
$\widehat{\bm{\beta}}(0,\lambda_{2},\mathbf{W})$ can be reformulated as a
Generalized LASSO problem. The Generalized LASSO, as introduced by Tibshirani
and Taylor, (2011), is formulated as follows:

$\widehat{\bm{\beta}}=\operatorname*{arg\,min}_{\bm{\beta}\in\mathds{R}^{p}}\leavevmode\nobreak\
\frac{1}{2}\left\lVert\mathbf{y}-\mathbf{X}\bm{\beta}\right\rVert_{2}^{2}+\gamma\left\lVert\mathbf{D}\bm{\beta}\right\rVert_{1}.\qquad\text{{(Generalized
LASSO)}}$ (7)

Here, $\mathbf{D}$ represents a real-valued matrix with $m$ rows and $p$
columns, and $\gamma\geq 0$ serves as the tuning parameter controlling the
penalty’s strength. Since we are concerned with transformation rather than
classic regression analysis, see Section 1, the design matrix $\mathbf{X}$ is
equal to $\mathbf{I}$. We demonstrate that it is possible to choose the matrix
$\mathbf{D}$ such that the sparsity penalty term from the Generalized LASSO is
equal to the smoothness penalty term of the wFLSA. We do this by using the
concept of the _oriented incidence matrix_ (Johnson et al.,, 1973; Tibshirani
and Taylor,, 2011).

Let $\widetilde{G}$ denote a directed version of the undirected graph $G$.
Specifically, let $\widetilde{w}$ represent the weight function assigning
weights to each edge in the graph $\widetilde{G}$ as follows:

$\widetilde{w}(\\{i,j\\})=\begin{cases}w(\\{i,j\\})=w_{ij}&\text{ if
}i<j\text{ and }\\\ 0&\text{otherwise},\end{cases}$ (8)

where $w$ is the weight function of the original graph $G$. In other words,
the edge $\\{i,j\\}$ in the original graph translates to the directed edge
$i\rightarrow j$ in $\widetilde{G}$ for $i<j$. We now define
$\mathbf{D}\in\mathds{R}^{\binom{p}{2}\times p}$ to be the oriented incidence
matrix of the graph $\widetilde{G}$, where each row represents a potential
edge $\\{i,j\\}\in V\times V$. The columns correspond to the $p$
vertices/coefficients. Let $\mathbf{d}_{ij}$ denote the $p$-dimensional row
vector related to the edge $\\{i,j\\}$, such that

$\mathbf{d}_{ij}=\Big{[}d_{ij}(1),d_{ij}(2),\ldots,d_{ij}(p)\Big{]},\qquad\text{where}\quad
d_{ij}(l)=\begin{cases}w_{ij}&\text{if }l=i,\\\ -w_{ij}&\text{if }l=j,\text{
and }\\\ 0&\text{otherwise}.\end{cases}$ (9)

The oriented incidence matrix $\mathbf{D}$ is then a matrix of these row
vectors, one for each potential edge, i.e.,

$\mathbf{D}=\begin{bmatrix}\mathbf{d}_{12}\\\ \mathbf{d}_{13}\\\
\mathbf{d}_{14}\\\ \vdots\\\ \mathbf{d}_{p-2,p-1}\\\
\mathbf{d}_{p-1,p}\end{bmatrix}.$ (10)

For an example demonstrating the construction of an oriented incidence matrix
$\mathbf{D}$ when there are four coefficients ($p=4$), see Appendix A. By
defining the matrix $\mathbf{D}$ as such, we obtain

$\left\lVert\mathbf{D}\bm{\beta}\right\rVert_{1}=\sum_{i<j}w_{ij}|\beta_{i}-\beta_{j}|.$
(11)

Combining these results, the Generalized LASSO problem is equivalent to a
wFLSA problem when the following conditions are met:

  1. 1.

the sparsity tuning parameter $\lambda_{1}$ equals $0$;

  2. 2.

the design matrix $\mathbf{X}$ is the identity matrix $\mathbf{I}$;

  3. 3.

the matrix $\mathbf{D}$ is the oriented incidence matrix of the weighted
directed graph $\widetilde{G}$ associated with the original wFLSA problem, and

  4. 4.

the tuning parameter of the Generalized LASSO $\gamma$ equals the smoothness
penalty term of the wFLSA, i.e., $\gamma=\lambda_{2}$.

In conclusion,

$\begin{split}\widehat{\bm{\beta}}(0,\lambda_{2},\mathbf{W})&=\operatorname*{arg\,min}_{\bm{\beta}\in\mathds{R}^{p}}\leavevmode\nobreak\
\frac{1}{2}\left\lVert\mathbf{y}-\bm{\beta}\right\rVert_{2}^{2}+\lambda_{2}\sum_{i<j}w_{ij}|\beta_{i}-\beta_{j}|\\\
&=\operatorname*{arg\,min}_{\bm{\beta}\in\mathds{R}^{p}}\leavevmode\nobreak\
\frac{1}{2}\left\lVert\mathbf{y}-\bm{\beta}\right\rVert_{2}^{2}+\lambda_{2}\left\lVert\mathbf{D}\bm{\beta}\right\rVert_{1}.\end{split}$
(12)

As mentioned earlier, showing that a wFLSA can be reformulated as a
Generalized LASSO problem in principle resolves the issue, given the existence
of algorithms for solving the Generalized LASSO (Arnold and Tibshirani,, 2016;
Zhao and Bondell,, 2020). Nonetheless, their computational complexity is
$\mathcal{O}(N_{\text{GL}}\cdot p^{4})$. To overcome this limitation, we
introduce an ADMM tailored to efficiently tackle wFLSA-equivalent problems in
the following section.

##  3 An Alternating Direction Method of Multipliers for wFLSA-equivalent
Problems

The results presented in this section rely heavily on the work done by Zhu,
(2017), wherein an ADMM tailored for the Generalized LASSO is developed. In an
ADMM algorithm, the original problem is decomposed into three smaller, more
manageable steps, see Section 1. However, Zhu, (2017) demonstrates that in the
context of the Generalized LASSO, the algorithm can be simplified, reducing
the number of update steps to two. We refer to these two steps as the
$\bm{\beta}$\- and $\bm{\alpha}$-update step.

For each update step indexed by $k=0,1,2,\ldots$, we consider two vectors:
$\bm{\beta}(k)$ and $\bm{\alpha}(k)$. The former is a $p$-dimensional vector,
while the latter is a vector whose dimensionality equals the number of rows of
matrix $\mathbf{D}$, i.e., $\binom{p}{2}$.

#### The $\bm{\beta}$-update step:

The $\bm{\beta}$-update step, as proposed by Zhu, (2017), is expressed as

$\bm{\beta}(k+1)=\big{(}\rho\mathbf{Q}+\mathbf{I}\big{)}^{-1}\Big{[}\rho\mathbf{Q}\bm{\beta}(k)+\mathbf{y}-\lambda_{2}\cdot\mathbf{D}^{\top}\big{(}2\bm{\alpha}(k)-\bm{\alpha}(k-1)\big{)}\Big{]}\quad\text{for
}k=1,2,\ldots.$

Here, $\mathbf{Q}$ denotes a $(p\times p)$-dimensional matrix, $\rho>0$ stands
for the ADMM’s dual update step length (Boyd and Vandenberghe,, 2004; Zhu,,
2017), $\mathbf{y}$ is the $p$-dimensional vector of observations,
$\lambda_{2}$ is the smoothness penalty term and $\mathbf{D}$ is the oriented
incidence matrix corresponding to the wFLSA problem under consideration.

The choice of matrix $\mathbf{Q}$ plays a critical role, as it can
significantly accelerate the algorithm when chosen appropriately (Zhu,, 2017).
It must satisfy the condition that the matrix
$\mathbf{Q}-\lambda_{2}^{2}\cdot\mathbf{D}^{\top}\mathbf{D}\succ 0$ is
positive definite. For simplicity, we choose $\mathbf{Q}=q\mathbf{I}$, i.e., a
matrix with the constant $q>0$ along the diagonal. To ensure positive
definiteness, $q$ must exceed the maximum eigenvalue of
$\lambda_{2}^{2}\cdot\mathbf{D}^{\top}\mathbf{D}$, which can be computed
numerically. However, performing the matrix multiplication
$\mathbf{D}^{\top}\mathbf{D}$ can be challenging due to the potentially large
number of rows in $\mathbf{D}$. For instance, with $p=1{,}000$, the row count
of $\mathbf{D}$ is close to half a million; for $p=10{,}000$, it approaches 50
million, and so on. Fortunately, $\mathbf{D}^{\top}\mathbf{D}$ possesses a
distinct structure that substantially reduces the computational time required
(see Appendix B). By defining $\mathbf{Q}$ as $q\mathbf{I}$, the update step
for $\bm{\beta}$ simplifies to

$\bm{\beta}(k+1)=\frac{1}{\rho q+1}\Big{[}\rho
q\bm{\beta}(k)+\mathbf{y}-\lambda_{2}\cdot\mathbf{D}^{\top}\big{(}2\bm{\alpha}(k)-\bm{\alpha}(k-1)\big{)}\Big{]}\quad\text{for
}k=1,2,\ldots$ (13)

#### The $\bm{\alpha}$-update step:

Following the $\bm{\beta}$ update, we proceed by updating the value of
$\bm{\alpha}$. As in the work by Zhu, (2017), the $\bm{\alpha}$-update step is
expressed as:

$\bm{\alpha}(k+1)=\text{T}\big{(}\bm{\alpha}(k)+\rho\cdot\lambda_{2}\cdot\mathbf{D}\bm{\beta}(k+1)\big{)}.$
(14)

Here, $\rho$ is again the dual update step length. The function
$\text{T}(\cdot):\mathds{R}^{\binom{p}{2}}\rightarrow[-1,1]^{\binom{p}{2}}$
denotes a threshold function defined as

$\text{T}\big{(}x_{1},x_{2},\ldots,x_{m}\big{)}=\big{(}\widetilde{x}_{1},\widetilde{x}_{2},\ldots,\widetilde{x}_{m}\big{)}^{\top}$
(15)

where $\widetilde{x}_{i}$ is determined as

$\widetilde{x}_{i}=\begin{cases}-1&\text{if }x_{i}<-1,\\\ 1&\text{if
}x_{i}>1\text{ and }\\\ x_{i}&\text{otherwise},\end{cases}$

for $i=1,2,\ldots,m$. This function restricts each component of the input
vector to the range $[-1,1]$, setting values below $-1$ to $-1$, values above
$1$ to $1$, and leaving values within the range unchanged. We iterate through
the update steps for $\bm{\beta}$ and $\bm{\alpha}$ until a stopping condition
is met. The algorithm terminates when the sum of the absolute differences
between $\bm{\beta}(k+1)$ and $\bm{\beta}(k)$ falls below a predefined
threshold $\epsilon>0$:

$\left\lVert\bm{\beta}(k+1)-\bm{\beta}(k)\right\rVert_{1}<\epsilon.$ (16)

Combining the findings from this section with those from the previous
sections, we arrive at the following pseudo-algorithm:

  1. 1.

Let $\mathbf{y}\in\mathds{R}^{p}$ be the vector containing the observations,
and let $G=(V,E,w)$ denote the undirected graph associated with the wFLSA,
characterized by the weighted adjacency matrix $\mathbf{W}$. It is necessary
to select the tuning parameters $\lambda_{2}\geq 0$, the dual update step
length $\rho>0$ (by default $\rho=1$), and $\epsilon>0$, which is the
convergence threshold for the algorithm;

  2. 2.

Determine the directed graph $\widetilde{G}$, as per eq. (8), along with its
oriented incidence matrix $\mathbf{D}$ as defined in eq. (10);

  3. 3.

Compute numerically the largest eigenvalue $\kappa_{\text{max}}$ (or an upper
bound) of the matrix $\lambda_{2}^{2}\cdot\mathbf{D}^{\top}\mathbf{D}$. See
Appendix B for further details. Select $q$ to be greater than the largest
eigenvalue, i.e., $q>\kappa_{\text{max}}$. This choice guarantees that
$q\mathbf{I}-\lambda_{2}^{2}\cdot\mathbf{D}^{\top}\mathbf{D}$ is strictly
positive definite;

  4. 4.

Initialize the vectors $\bm{\beta}(0)=\bm{\beta}(1)=\mathbf{0}$ and
$\bm{\alpha}(0)=\bm{\alpha}(1)=\mathbf{0}$.

  5. 5.

Execute the $\bm{\beta}$-update step according to eq. (13);

  6. 6.

Execute the $\bm{\alpha}$-update step according to eq. (14), and

  7. 7.

If $\left\lVert\bm{\beta}(k+1)-\bm{\beta}(k)\right\rVert_{1}<\epsilon$,
terminate the algorithm and return $\bm{\beta}(k+1)$ as the coefficient vector
$\widehat{\bm{\beta}}(0,\lambda_{2},\mathbf{W})$. Otherwise, repeat steps 5
and 6.

The estimates $\widehat{\bm{\beta}}(\lambda_{1},\lambda_{2},\mathbf{W})$ with
$\lambda_{1}>0$ can be easily obtained by applying soft thresholding, see eq.
(6). An implementation of this algorithm is available at
https://github.com/bips-hb/wflsa.

##  4 Computational Complexity

To assess the computational complexity of ADMM, we analyze the update steps
for both $\bm{\beta}$ and $\bm{\alpha}$, separately.

###  4.1 The $\bm{\beta}$-Update Step

The most computationally intensive aspect of the $\bm{\beta}$-update step, see
eq. (13), is the matrix-vector multiplication

$\mathbf{D}^{\top}\big{(}2\bm{\alpha}(k)-\bm{\alpha}(k-1)\big{)}.$

Here, $\mathbf{D}^{\top}$ has $p$ rows and $\binom{p}{2}$ columns, and
$\bm{\alpha}(k)$ and $\bm{\alpha}(k-1)$ are $\binom{p}{2}$-dimensional
vectors. A naive approach would entail multiplying two
$\binom{p}{2}$-dimensional vectors $p$ times (once for every row in
$\mathbf{D}^{\top}$), resulting in a complexity of $\mathcal{O}(p^{3})$.

However, we can show that each row of the transposed oriented incidence matrix
$\mathbf{D}^{\top}$ contains at most $p-1$ non-zero entries. Hence, instead of
multiplying $\binom{p}{2}$ values, only $p-1$ multiplications are needed.
Given that this operation must be performed for each of the $p$ rows, the
complexity is thereby reduced to $\mathcal{O}(p^{2})$.

To confirm this assertion, let us start by denoting the $l$-th column vector
of matrix $\mathbf{D}$ by $\mathbf{e}_{l}\in\mathds{R}^{\binom{p}{2}}$:

$\mathbf{e}_{l}^{\top}=\Big{[}d_{12}(l),d_{13}(l),d_{14}(l),\ldots,d_{p-2,p}(l),d_{p-1,p}(l)\Big{]}^{\top}\quad\text{for
}l=1,2,\ldots,p.$ (17)

By using the definition of $d_{ij}(l)$ from eq. (9), it is evident that at
most $p-1$ entries can be non-zero (an entry is non-zero only if the
corresponding entry in the weight matrix is non-zero as well). This
observation leads us to conclude that each row of $\mathbf{D}^{\top}$ has at
most $p-1$ non-zero entries. Consequently, the expression
$\mathbf{e}_{l}^{\top}(2\bm{\alpha}(k)-\bm{\alpha}(k-1))$ involves at most
$p-1$ multiplications. Given there are $p$ rows, the complexity of the
$\bm{\beta}$-update step is indeed $\mathcal{O}(p^{2})$.

###  4.2 The $\bm{\alpha}$-Update Step

The most computationally intensive aspect of the $\bm{\alpha}$-update step,
see eq. (14), involves the matrix-vector multiplication
$\mathbf{D}\bm{\beta}(k+1)$. If approached naively, this operation would yield
a complexity of $\mathcal{O}(p^{3})$. The structural properties of the matrix
$\mathbf{D}$ can be utilized in this context as well. By employing the
definition of $\mathbf{D}$ in eq. (10), we can express

$\mathbf{D}\bm{\beta}(k+1)=\begin{bmatrix}\mathbf{d}_{12}\cdot\bm{\beta}(k+1)\\\
\mathbf{d}_{13}\cdot\bm{\beta}(k+1)\\\ \mathbf{d}_{14}\cdot\bm{\beta}(k+1)\\\
\vdots\\\ \mathbf{d}_{p-2,p}\cdot\bm{\beta}(k+1)\\\
\mathbf{d}_{p-1,p}\cdot\bm{\beta}(k+1)\end{bmatrix}$

where

$\mathbf{d}_{ij}\cdot\bm{\beta}(k+1)=w_{ij}\beta_{i}(k+1)-w_{ij}\beta_{j}(k+1)=w_{ij}\left(\beta_{i}(k+1)-\beta_{j}(k+1)\right)$

for all $i<j$. Each entry of the vector $\mathbf{D}\bm{\beta}(k+1)$ thus
requires a constant number of operations, specifically a summation and a
multiplication. Since the vector has $\binom{p}{2}$ entries, the complexity is
of the order $\mathcal{O}(p^{2})$. Consequently, computing
$\mathbf{D}\bm{\beta}(k+1)$ requires $\mathcal{O}(p^{2})$ computations.

As mentioned earlier, selecting an appropriate value for $q>0$ is crucial to
ensure that the matrix $\lambda_{2}^{2}\cdot\mathbf{D}^{\top}\mathbf{D}$ is
positive definite. Computing the largest eigenvalue of this matrix numerically
requires $\mathcal{O}(p^{3})$ operations, which is manageable for small values
of $p$. However, as $p$ increases, this computational cost can become
prohibitive. Fortunately, it suffices to establish an upper bound for the
largest eigenvalue, which can be achieved using the Gershgorin circle theorem
(Gloub and Van Loan,, 1996). Importantly, this process can be accomplished
with a complexity of $\mathcal{O}(p^{2})$. Further information is provided in
Appendix B.

In conclusion, given that both the $\bm{\beta}$\- and $\bm{\alpha}$-update
steps have a complexity of $\mathcal{O}(p^{2})$, a complete single update step
has a complexity of $\mathcal{O}(p^{2})$ as well. Let $N_{\text{ADMM}}$ be the
number of iterations required to satisfy the stopping criterion defined in eq.
(16), the overall computational complexity of the ADMM algorithm is
$\mathcal{O}(N_{\text{ADMM}}\cdot p^{2})$.

##  5 Tuning Parameter and Weight Matrix Selection

The above-described ADMM algorithm requires the selection of appropriate
values for:

  1. 1.

the tuning parameter $\lambda_{1}$, which controls the sparsity of the
coefficient vector;

  2. 2.

the tuning parameter $\lambda_{2}$, which manages the level of smoothness or
similarity among different coefficients, and

  3. 3.

a weight matrix $\mathbf{W}$, determining which coefficients are to be similar
and to what degree.

Choosing suitable values for these parameters is challenging, and only a few
research papers have addressed this issue. For instance, Gao, (2021) proposed
using the ‘classical’ approaches AIC, BIC, and General Cross Validation (GCV)
for their Least Absolute Derivations Fused LASSO Signal Approximator, but did
not provide a comprehensive comparison. Son and Lim, (2019) compared the BIC
and the extended BIC (eBIC) for the FLSA, where the latter imposes a stronger
penalty on larger models. As a result, eBIC achieved a better false discovery
rate with only a slightly lower true positive rate. Recently, Son et al.,
(2023) proposed a two-step selection procedure based on a General Information
Criterion (GIC). In the first step, GIC is used to select $\lambda_{1}$, and
in the second step, a preliminary test statistic with a piece-wise mean
structure is used to select $\lambda_{2}$. Their numerical results show a good
false discovery rate under the assumption of low (Gaussian) noise. For the
wFLSA, we also need to search over a finite number of weight matrices in
$\mathbf{W}\in\mathcal{W}$. This makes approaches such as cross-validation or
stability selection quite challenging due to their subsampling related
algorithms. Moreover, even for information criteria such as eBIC, the
computational cost seems to be quiet high.

##  6 Heterogeneous Image Smoothing

In addition to its originally intended applications for regression problems,
the smoothing property can also be utilized for image processing. Generally,
LASSO or Fused LASSO have already been employed for sparse coding or image
denoising (Li and Shouno,, 2014; Agarwal et al.,, 2007). However, these
methods, along with many traditional Bayesian approaches, often assume the
presence of additive zero-mean and homogeneous Gaussian noise in the image
(Elad et al.,, 2023). This means the noise intensity is uniform across all
pixels, commonly referred to as additive white Gaussian noise. In reality,
images often exhibit pixel-dependent heterogeneous noise structures, which
also pose challenges for the deep learning methods that have recently become
prevalent (Zhou et al.,, 2019).

By using the weight matrix $\mathbf{W}$, prior knowledge can be incorporated
to apply the smoothing effect only to certain regions rather than uniformly
throughout the image. This approach prevents the blurring of completely clear
regions in the presence of heterogeneous or spatially variant noise, a common
issue with methods like the median filter. However, applying wFLSA to high-
resolution images is not feasible, as the algorithm’s complexity grows
quadratically with the number of pixels (see Section 4). Given the plausible
assumption of locality that pixels only interact with their local neighbors,
this problem can be reduced to a patch-wise evaluation, which is a common
technique in image processing, especially for deep learning-based denoisers.
This technique divides an image into smaller sub-images (called patches) of a
predefined size (e.g., $5\times 5$) and applies an operation (usually called a
filter or kernel) to each patch. By sliding the patch frame across the image,
e.g., by moving it one pixel to the right, the results can be averaged for the
overlapping pixels, providing an additional smoothing effect (Szeliski,,
2011). In this context, wFLSA is applied to all image patches with the
corresponding weight patches of $\mathbf{W}$ and then convoluted over the
entire high-dimensional image.

(a) Original${}^{\text{1}}$

(b) wFLSA ($\lambda_{2}=.04$)

(c) Median filter

(d) Radial noise

(e) wFLSA ($\lambda_{2}=.1$)

(f) NLM

Figure 1: Application of wFLSA for image smoothing with prior knowledge on the
noise structure. Image (a) is the original (a photo of Taylor Swift), and (d)
is the same image with added radial noise. Images (b) and (e) are the
reconstructions using wFLSA with different smoothing parameters $\lambda_{2}$.
Image (c) uses a $5\times 5$ median filter, and (f) employs the non-local
means (NLM) method. 11footnotetext: Image source:
https://de.m.wikipedia.org/wiki/Datei:Taylor_Swift_2_-_2019_by_Glenn_Francis_(cropped).jpg

To illustrate this property, a grayscale image of size $456\times 400$ is
taken, and additive noise is added, which increases radially from the center
towards the corners, resulting in higher noise levels at the edges. This
intensity matrix is then used as the weight matrix $\mathbf{W}$ to apply a
stronger smoothing effect in the corners compared to the central region of the
image. A patch size of $5\times 4$ is chosen for the convolution.
Additionally, the denoised results are visually compared with those obtained
using the median filter and the non-local means (NLM) denoising algorithm
(Buades et al.,, 2011). The tuning parameter $\lambda_{2}$ was set to $.04$
and $.1$, varying the level of smoothness. The sparsity parameter
$\lambda_{1}$ was fixed to $.001$ in both cases. The results are presented in
Figure 1, showing the trade-off between reconstructing and smoothing the
image. The figure can be reproduced using the code available at
https://github.com/bips-hb/wFLSA-paper.

##  7 Conclusions and Discussion

We proposed an ADMM algorithm tailored for the Weighted Generalized Fused
Lasso Signal Approximator (wFLSA). Our approach involved two key steps. First,
we have shown that a wFLSA problem could be reformulated as a Generalized
LASSO by leveraging the oriented incidence matrix of the directed version of
the underlying graph of the wFLSA, see Section 2. This step theoretically
resolves the issue, as path algorithms exists for solving the Generalized
LASSO (Arnold and Tibshirani,, 2016; Zhao and Bondell,, 2020). However, the
computational complexity of these algorithms is prohibitively high, i.e.,
$\mathcal{O}(N_{\text{GL}}\cdot p^{4})$ where $N_{\text{GL}}$ is the number of
iterations needed for the algorithm to converge, which makes real-world
applications challenging. To address this, we developed an ADMM approach,
relying on the groundwork laid by Zhu, (2017). We established that our ADMM-
based solution drastically reduces complexity to
$\mathcal{O}(N_{\text{ADMM}}\cdot p^{2})$, see Section 4.

We discuss various methodologies for selecting the model’s tuning parameters
and weight matrix in Section 5. This is far from trivial since the
computational complexity of the ADMM algorithm makes a broad search space
rather challenging. Here, we propose a different approach that combines a
heuristic method with the use of an information criterion. Specifically, we
first select a set of ‘reasonable’ weight matrices ($\mathcal{W}$), possibly
based on prior knowledge about a network structure. Next, for each weight
matrix $\mathbf{W}\in\mathcal{W}$, a $\lambda_{1}$-value is chosen using the
General Information Criterion (GIC). Finally, based on the $\mathbf{W}$ and
$\lambda_{1}$-value that yield the smallest GIC, $\lambda_{2}$ is determined
using a preliminary test statistic as suggested by Son et al., (2023). It
would be interesting for future research to compare this approach with the
existing methodologies discussed in Section 5.

In Section 6, we demonstrated how wFLSA can efficiently be used for spatially
varying noise reduction by using prior knowledge on the noise intensities in
the weight matrix. Our approach relaxes the very common assumption of white
additive Gaussian noise and does not require an extensive training procedure,
as for many deep learning techniques.

Hoefling, (2010) suggested the possibility of extending his path algorithm
designed for solving the Generalized FLSA to the weighted version. However,
how to adapt the algorithm remains unclear.

Considering the wFLSA problem as a Generalized LASSO problem offers a
particular advantage. The algorithms proposed by Arnold and Tibshirani,
(2016); Zhao and Bondell, (2020) are _path algorithms_ , providing the entire
solution path for the smoothness penalty term $\lambda_{2}$, i.e.,
$\widehat{\bm{\beta}}(0,\lambda_{2},\mathbf{W})$ for all $\lambda_{2}\geq 0$.
This is not the case for the ADMM presented in this paper; the value of
$\lambda_{2}$ is fixed and chosen beforehand. The proof of Friedman et al.,
(2007), as discussed in Section 2, shows that once the estimate of the
coefficients for a specific value of $\lambda_{2}$ is known, the solution path
for $\lambda_{1}\geq 0$ is known as well. Combining these results means that
when employing the Generalized LASSO, one can obtain the entire _solution
plane_ (!), encompassing solutions for all possible values of $\lambda_{1}$
and $\lambda_{2}$. However, due to the computational complexity associated
with the Generalized LASSO, using this approach may not always be feasible.

The implementation of the application example from Section 6 can be found at
https://github.com/bips-hb/wFLSA-paper. The ADMM algorithm has been
implemented as an R package and is publicly available under
https://github.com/bips-hb/wflsa.

## Acknowledgments

We gratefully acknowledge the financial support of the German Research
Foundation (DFG – Project FO 1045/2-1).

##  Appendix A An Example: wFLSA as Generalized LASSO with $p=4$

In Section 2, we demonstrate that the wFLSA problem in eq. (5) can be
reformulated as the Generalized LASSO problem as defined in eq. (7). Let us
illustrate how a wFLSA problem can be translated into an equivalent
Generalized LASSO problem when the number of coefficients is $p=4$. We argue
that the penalty terms of the wFLSA can be expressed as

$\lambda_{2}\sum_{i=1}^{3}\sum_{j=i+1}^{4}w_{ij}|\beta_{i}-\beta_{j}|=\lambda_{2}\left\lVert\mathbf{D}\bm{\beta}\right\rVert_{1}.$

Here, $\lambda_{2}$ represent the smoothness penalty term,
$\bm{\beta}=(\beta_{1},\beta_{2},\beta_{3},\beta_{4})^{\top}$ denotes the
vector of coefficients, $w_{ij}\geq 0$ indicates the weight of the edge
between the $i$-th and $j$-th node in the graph $G$, and $\mathbf{D}$
represents the oriented incidence matrix. By using eq. (10), we find that the
oriented incidence matrix equals

$\mathbf{D}=\begin{bmatrix}w_{12}&-w_{12}&0&0\\\ w_{13}&0&-w_{13}&0\\\
w_{14}&0&0&-w_{14}\\\ 0&w_{23}&-w_{23}&0\\\ 0&w_{24}&0&-w_{24}\\\
0&0&w_{34}&-w_{34}\end{bmatrix}.$

Expanding the expression $\mathbf{D}\bm{\beta}$, we observe that

$\mathbf{D}\bm{\beta}=\left[\begin{array}[]{c c c c}w_{12}&-w_{12}&0&0\\\
w_{13}&0&-w_{13}&0\\\ w_{14}&0&0&-w_{14}\\\ 0&w_{23}&-w_{23}&0\\\
0&w_{24}&0&-w_{24}\\\ 0&0&w_{34}&-w_{34}\\\
\end{array}\right]\begin{bmatrix}\beta_{1}\\\ \beta_{2}\\\ \beta_{3}\\\
\beta_{4}\end{bmatrix}=\left[\begin{array}[]{c}w_{12}(\beta_{1}-\beta_{2})\\\
w_{13}(\beta_{1}-\beta_{3})\\\ w_{14}(\beta_{1}-\beta_{4})\\\
w_{23}(\beta_{2}-\beta_{3})\\\ w_{24}(\beta_{2}-\beta_{4})\\\
w_{34}(\beta_{3}-\beta_{4})\end{array}\right].$

The $\ell_{1}$-norm of $\mathbf{D}\bm{\beta}$ is the summation of the absolute
values of the entries of the vector on the right-hand side, which yields

$\displaystyle\lambda_{2}\sum_{i<j}|w_{ij}(\beta_{i}-\beta_{j})|=\lambda_{2}\sum_{i<j}w_{ij}|\beta_{i}-\beta_{j}|.$

This equality holds as the weights are non-negative values.

##  Appendix B Choice of Matrix $\mathbf{Q}$ for the ADMM Algorithm

The ADMM algorithm introduced in this paper requires selecting a $(p\times
p)$-dimensional matrix $\mathbf{Q}$ such that
$\mathbf{Q}-\lambda_{2}^{2}\cdot\mathbf{D}^{\top}\mathbf{D}$ is positive
definite, see Section 3 and, more specifically, Zhu, (2017). To make the
problem more tractable, we define $\mathbf{Q}$ as a diagonal matrix with a
constant $q>0$ on its diagonal, i.e., $\mathbf{Q}=q\mathbf{I}$. Determining
suitable values of $q$ can be done numerically; if $\kappa_{\text{max}}$ is
the largest eigenvalue of $\lambda_{2}^{2}\cdot\mathbf{D}^{\top}\mathbf{D}$,
then any $q>\kappa_{\text{max}}$ satisfies the condition. However, computing
the product $\mathbf{D}^{\top}\mathbf{D}$ can be challenging due to the
potentially large number of rows of $\mathbf{D}$. Fortunately,
$\mathbf{D}^{\top}\mathbf{D}$ exhibits a distinct structure.

Recall that we can express the incidence matrix
$\mathbf{D}=[\mathbf{e}_{1},\mathbf{e}_{2},\ldots,\mathbf{e}_{p}]$ in terms of
its $p$ column vectors, see eq. (17),

$\mathbf{e}_{l}=\begin{bmatrix}d_{12}(l)\\\ d_{13}(l)\\\ d_{14}(l)\\\
\vdots\\\ d_{p-2,p}(l)\\\ d_{p-1,p}(l)\end{bmatrix}\qquad\text{for
}l=1,2,\ldots,p.$

By doing so, the premultiplication of $\mathbf{D}$ by its transpose can be
written as

$\mathbf{D}^{\top}\mathbf{D}=\begin{bmatrix}\mathbf{e}_{1}^{\top}\mathbf{e}_{1}&\mathbf{e}_{1}^{\top}\mathbf{e}_{2}&\ldots&\mathbf{e}_{1}^{\top}\mathbf{e}_{p}\\\
\mathbf{e}_{2}^{\top}\mathbf{e}_{1}&\mathbf{e}_{2}^{\top}\mathbf{e}_{2}&\ldots&\mathbf{e}_{2}^{\top}\mathbf{e}_{p}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{e}_{p}^{\top}\mathbf{e}_{1}&\mathbf{e}_{p}^{\top}\mathbf{e}_{2}&\ldots&\mathbf{e}_{p}^{\top}\mathbf{e}_{p}\end{bmatrix}.$

Let us start by examining the diagonal elements. For $l=1,2,\ldots,p$, the
diagonal entry is given by

$\bar{\mathbf{w}}_{l}=\mathbf{e}_{l}^{\top}\mathbf{e}_{l}=\sum_{i<j}d_{ij}(l)d_{ij}(l)=\sum_{i<j}d_{ij}^{2}(l).$

Using the definition of $d_{ij}(l)$ from eq. (9), we deduce that

$d_{ij}^{2}(l)=\begin{cases}w_{il}^{2}&\text{if }l=j,\\\ w_{lj}^{2}&\text{if
}l=i,\text{ and}\\\ 0&\text{otherwise}.\end{cases}$

This leads to a simpler expression for the diagonal value:

$\bar{\mathbf{w}}_{l}=\mathbf{e}_{l}^{\top}\mathbf{e}_{l}=\sum_{i=1}^{l-1}w_{il}^{2}+\sum_{j=l+1}^{p}w_{lj}^{2}.$

Due to the symmetry of the weight matrix ($w_{ij}=w_{ji}$) and its zero
diagonal ($w_{ll}=0$), we can express this as

$\bar{\mathbf{w}}_{l}=\mathbf{e}_{l}^{\top}\mathbf{e}_{l}=\sum_{i=1}^{p}w_{il}^{2}.$

In other words, the $l$-th entry on the diagonal of the matrix
$\mathbf{D}^{\top}\mathbf{D}$ equals the sum of squared values of the $l$-th
column (or row) of the weighted adjacency matrix $\mathbf{W}$.

The off-diagonal values are rather straightforward. For $s\neq t$, we have:

$\mathbf{e}_{s}^{\top}\mathbf{e}_{t}=\sum_{i<j}d_{ij}(s)d_{ij}(t)=w_{st}\cdot(-w_{ts})=-w_{st}^{2}.$

Combining these results, the matrix $\mathbf{D}^{\top}\mathbf{D}$ can be
expressed as

$\mathbf{D}^{\top}\mathbf{D}=\text{diag}(\bar{\mathbf{w}}_{1},\bar{\mathbf{w}}_{2},\ldots,\bar{\mathbf{w}}_{p})-\mathbf{W}\circ\mathbf{W},$

where $\circ$ denotes element-wise multiplication. Numerically determining
this matrix is computationally much more efficient than directly
premultiplying $\mathbf{D}$ by $\mathbf{D}^{\top}$.

Consider that computing the largest eigenvalue of a $(p\times p)$-dimensional
matrix is computationally intensive, scaling as $\mathcal{O}(p^{3})$. While
this poses no significant challenge for smaller values of $p$, it becomes a
notable bottleneck for larger dimensions. Instead of directly estimating the
largest eigenvalue, we can also employ an upper bound.

According to the Gershgorin circle theorem (Gloub and Van Loan,, 1996), the
largest eigenvalue is upper bounded by the maximum absolute row sum of the
matrix $\mathbf{D}^{\top}\mathbf{D}$. Using our previous result, we find that
the sum of the absolute values of the $i$-th row is simply
$2\bar{\mathbf{w}}_{i}$. An upper bound for the largest eigenvalue is,
therefore,

$\kappa_{\text{max}}\leq 2\cdot\text{max}\\{\bar{\mathbf{w}}_{i}\text{ for
}i=1,2,\ldots,p\\}.$

The computational complexity involved in determining this upper bound is of
the order $\mathcal{O}(p^{2})$.

## References

  * Agarwal et al., (2007) Agarwal, V., Gribok, A. V., and Abidi, M. A. (2007).  Image restoration using $l_{1}$ norm penalty function.  Inverse Problems in Science and Engineering, 15(8):785–809.
  * Arnold and Tibshirani, (2016) Arnold, T. B. and Tibshirani, R. J. (2016).  Efficient implementations of the generalized Lasso dual path algorithm.  Journal of Computational and Graphical Statistics, 25(1):1–27.
  * Boyd and Vandenberghe, (2004) Boyd, S. and Vandenberghe, L. (2004).  Convex optimization.  Cambridge University Press.
  * Buades et al., (2011) Buades, A., Coll, B., and Morel, J.-M. (2011).  Non-local means denoising.  Image Processing On Line, 1:208–212.
  * Condat, (2013) Condat, L. (2013).  A direct algorithm for 1-D total variation denoising.  IEEE Signal Processing Letters, 20(11):1054–1057.
  * Elad et al., (2023) Elad, M., Kawar, B., and Vaksman, G. (2023).  Image denoising: The deep learning revolution and beyond – A survey paper.  SIAM Journal on Imaging Sciences, 16(3):1594–1654.
  * Forero et al., (2010) Forero, P. A., Cano, A., and Giannakis, G. B. (2010).  Consensus-based distributed linear support vector machines.  In Proceedings of the 9th ACM/IEEE International Conference on Information Processing in Sensor Networks, pages 35–46.
  * Friedman et al., (2007) Friedman, J., Hastie, T., Höfling, H., and Tibshirani, R. (2007).  Pathwise coordinate optimization.  The Annals of Applied Statistics, 1:302–332.
  * Gao, (2021) Gao, X. (2021).  Estimation and selection properties of the LAD fused Lasso signal approximator.
  * Gloub and Van Loan, (1996) Gloub, G. H. and Van Loan, C. F. (1996).  Matrix computations.  Johns Hopkins Universtiy Press, 3rd edtion.
  * Glowinski et al., (2022) Glowinski, R., Song, Y., Yuan, X., and Yue, H. (2022).  Application of the alternating direction method of multipliers to control constrained parabolic optimal control problems and beyond.  Annals of Applied Mathematics, 38(2):115–158.
  * Harris et al., (2022) Harris, T., Li, B., and Tucker, J. D. (2022).  Scalable multiple changepoint detection for functional data sequences.  Environmetrics, 33(2):e2710.
  * Hoefling, (2010) Hoefling, H. (2010).  A path algorithm for the fused Lasso signal approximator.  Journal of Computational and Graphical Statistics, 19(4):984–1006.
  * Johnson et al., (1973) Johnson, D. E., Johnson, J. R., and Chen, W.-K. (1973).  Graph theory with engineering applications.  IEEE Transactions on Systems, Man, and Cybernetics, (5):530–530.
  * Kim et al., (2009) Kim, S.-J., Koh, K., Boyd, S., and Gorinevsky, D. (2009).  $\ell_{1}$ trend filtering.  SIAM Review, 51(2):339–360.
  * Li and Shouno, (2014) Li, A. and Shouno, H. (2014).  Dictionary-based image denoising by fused-Lasso atom selection.  Mathematical Problems in Engineering, 2014(1):368602.
  * Qian and Jia, (2016) Qian, J. and Jia, J. (2016).  On stepwise pattern recovery of the fused Lasso.  Computational Statistics & Data Analysis, 94:221–237.
  * Rojas and Wahlberg, (2014) Rojas, C. R. and Wahlberg, B. (2014).  On change point detection using the fused Lasso method.  arXiv preprint arXiv:1401.5408.
  * Son and Lim, (2019) Son, W. and Lim, J. (2019).  Modified path algorithm of fused Lasso signal approximator for consistent recovery of change points.  Journal of Statistical Planning and Inference, 200:223–238.
  * Son et al., (2023) Son, W., Lim, J., and Yu, D. (2023).  Tuning parameter selection in fused Lasso signal approximator with false discovery rate control.  Brazilian Journal of Probability and Statistics, 37(3).
  * Szeliski, (2011) Szeliski, R. (2011).  Computer vision: Algorithms and applications.  Springer London.
  * Tibshirani, (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the Lasso.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 58(1):267–288.
  * Tibshirani et al., (2005) Tibshirani, R., Saunders, M., Rosset, S., Zhu, J., and Knight, K. (2005).  Sparsity and smoothness via the fused Lasso.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 67(1):91–108.
  * Tibshirani and Wang, (2008) Tibshirani, R. and Wang, P. (2008).  Spatial smoothing and hot spot detection for CGH data using the fused Lasso.  Biostatistics, 9 1:18–29.
  * Tibshirani and Taylor, (2011) Tibshirani, R. J. and Taylor, J. (2011).  The solution path of the generalized Lasso.  The Annals of Statistics, 39(3):1335 – 1371.
  * Wahlberg et al., (2012) Wahlberg, B., Boyd, S., Annergren, M., and Wang, Y. (2012).  An ADMM algorithm for a class of total variation regularized estimation problems.  IFAC Proceedings Volumes, 45(16):83–88.
  * Xu and Noo, (2017) Xu, J. and Noo, F. (2017).  A sequential solution for anisotropic total variation image denoising with interval constraints.  Physics in Medicine & Biology, 62(18):N428.
  * Yang et al., (2018) Yang, Y., Sun, J., Li, H., and Xu, Z. (2018).  ADMM-CSNet: A deep learning approach for image compressive sensing.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 42(3):521–538.
  * Zhao and Bondell, (2020) Zhao, Y. and Bondell, H. (2020).  Solution paths for the generalized Lasso with applications to spatially varying coefficients regression.  Computational Statistics & Data Analysis, 142:106821.
  * Zhou et al., (2019) Zhou, Y., Jiao, J., Huang, H., Wang, Y., Wang, J., Shi, H., and Huang, T. S. (2019).  When AWGN-based denoiser meets real noises.  In AAAI Conference on Artificial Intelligence.
  * Zhu, (2017) Zhu, Y. (2017).  An augmented ADMM algorithm with application to the generalized Lasso problem.  Journal of Computational and Graphical Statistics, 26(1):195–204.