# Ordinary Differential Equation-based
Sparse Signal Recovery

Tadashi Wadayama and Ayano Nakai-Kasai
Part of this research was presented at the IEICE International Symposium on
Information Theory and Its Applications 2022 (ISITA2022) [1]. 1Nagoya
Institute of Technology, Gokiso, Nagoya, Aichi 466-8555, Japan,
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This study investigates the use of continuous-time dynamical systems for
sparse signal recovery. The proposed dynamical system is in the form of a
nonlinear ordinary differential equation (ODE) derived from the gradient flow
of the Lasso objective function. The sparse signal recovery process of this
ODE-based approach is demonstrated by numerical simulations using the Euler
method. The state of the continuous-time dynamical system eventually converges
to the equilibrium point corresponding to the minimum of the objective
function. To gain insight into the local convergence properties of the system,
a linear approximation around the equilibrium point is applied, yielding a
closed-form error evolution ODE. This analysis shows the behavior of
convergence to the equilibrium point. In addition, a variational optimization
problem is proposed to optimize a time-dependent regularization parameter in
order to improve both convergence speed and solution quality. The deep
unfolded-variational optimization method is introduced as a means of solving
this optimization problem, and its effectiveness is validated through
numerical experiments.

##  I Introduction

Continuous-time dynamical systems are ubiquitously used in signal processing.
Electronic circuits are a common example where the behavior is represented by
nonlinear ordinary differential equations (ODEs). These circuits typically
handle low-dimensional signals and simple signal processing tasks. From both a
practical and theoretical point of view, continuous-time systems offer a
number of advantages when processing high-dimensional signals. One such
advantage is the potential for analog-domain implementations of high-
dimensional signal processing tasks that are typically performed in the
digital domain.

Recently, the trend in Moore’s law, which asserts that the number of
transistors in an integrated circuit will double every 24 months, has not been
reflected in reality. A major limiting factor is the power consumption of the
integrated circuit, which hampers the ability to handle large-scale signal
processing problems that require high throughput. For instance, base stations
(BSs) in wireless network systems [4][5] need to perform a significant number
of signal processing tasks, including multiple-input and multiple-output
(MIMO) signal detection and beamforming. In next-generation systems (e.g.,
those beyond 5G or 6G), the limitations of the central processing unit (CPU)
are likely to become a major obstacle to achieving desired specifications. A
potential solution to this challenge would be to use an analog-domain signal
processing unit as a specialized processing unit with lower power consumption.

Analog computing [6][7][8] can be effectively employed in a variety of
domains, including the electrical and optical domains. An analog computer
capable of simulating any nonlinear ODE is composed of analog adders,
multipliers, integrators, and other nonlinear devices. Operational amplifiers
can be used to implement such analog devices required for analog computing.
Prior studies have been conducted on the use of specialized analog domain
processing units [9, 10] based on these electrical devices to achieve better
energy efficiency. In the field of optical computing, optical adders and
multipliers have been used to solve ODEs [11]. Recent examples of optical
computing include a neural network implemented in the optical domain [12] and
an optical Ising machine [13] for combinatorial optimization. Programmable
integrated optical circuits [14] are becoming an active research topic in
optical computing. While it is too early to make a definitive statement, it
appears that analog computing has the potential to perform large-scale signal
processing tasks with high energy efficiency and extremely fast computation in
the future.

Another advantage of continuous-time dynamical systems is that they can
provide deeper insights into discrete-time algorithms, which are the
counterpart of continuous-time systems. They also offer a way to generate new
discrete-time algorithms. A numerical method for solving ODEs, such as the
Euler method or the Runge-Kutta method [16], can be used to discretize a high-
dimensional ODE corresponding to a signal processing task. This principle has
been proposed in the study of neural ODEs [2]. The discrete-time algorithms
derived in such a manner inherit properties of the original ODE. This
correspondence between an ODE and a discrete-time algorithm opens an
insightful perspective for analyzing discrete-time algorithms.

There have been several studies into ODE-based signal processing, including a
recent work [15] that proposed an ODE-based method for MMSE signal detection
in MIMO systems and presented a concise analytical formula for the mean
squared error. This work [15] deals with a linear ODE that is relatively easy
to analyze. Despite the usefulness, there remains a plethora of opportunities
for further research in nonlinear ODE-based signal processing.

The focus of this paper is on an ODE-based method for sparse signal recovery
[27, 28]. Our work is motivated by the fact that most tasks in sparse signal
recovery (or compressed sensing) are high-dimensional signal processing tasks
that are not trivial to solve. Furthermore, there is a growing demand for
high-speed and large-scale sparse signal recovery tasks in a compressed
sensing context [29].

A number of discrete time sparse signal recovery algorithms [29] have been
developed based on the Lasso formulation [30, 31, 32, 33]. The Iterative
Shrinkage Thresholding Algorithm (ISTA) [34, 35] is one of the best-known
algorithms for solving the Lasso problem. ISTA is an iterative algorithm
comprising two processes: a linear estimation process and a shrinkage process
based on a soft shrinkage function. ISTA can be seen as a proximal gradient
descent algorithm [36] and can be directly derived from the Lasso formulation.
Another powerful iterative algorithm is Approximate Message Passing (AMP) [37,
38] which provides much faster convergence than ISTA. Ma and Ping proposed
Orthogonal AMP (OAMP) [39], which can handle various classes of sensing
matrices including unitary invariant matrices. Rangan et al. proposed Vector
AMP [40] for right-rotationally invariant matrices and provided a theoretical
justification for its state evolution.

Although the advancement of these sparse signal recovery algorithms is
significant, all of them are discrete-time algorithms. The extension to
continuous-time signal processing is not straightforward. Our approach is
based on the Lasso formulation [30], which can be viewed as a typical
regularized least-squares (LS) objective function. The proposed ODE
corresponds to the gradient flow based on the Lasso objective function. To
gain insight into the local convergence properties of the system, a linear
approximation around the equilibrium point is applied, which yields a closed-
form error evolution ODE. The analysis tracks system behavior as the solution
converges to the equilibrium point. In addition, a variational optimization
problem is proposed to optimize a time-dependent regularization parameter in
order to improve both convergence speed and solution quality. A deep
unfolding-based method is presented for solving the variational problem. The
methodology presented in the paper is directly applicable to any type of
regularized LS problem.

The remainder of the paper is organized as follows: In Section II, the system
model and the Lasso formulation used in the paper are introduced. Section III
presents the proposed ODE-based method for sparse signal recovery and includes
numerical demonstrations. In Section IV, an analysis of the local convergence
of the proposed ODE-based method is given, with the key idea being the use of
a linear approximation around the equilibrium point. Numerical examples are
provided to support the theory presented in this section. In Section V, a
parametric ODE, which includes a continuous-time function to be optimized, is
introduced. To improve the convergence behavior and solution quality, a
variational optimization problem including the parametric ODE must be solved.
A deep unfolding-based optimization method is presented for solving this
problem. Section VI concludes the discussion.

##  II Preliminaries

###  II-A Notation

The following mathematical notation is used throughout the paper: The symbols
$\mathbb{R}$ and $\mathbb{R}_{+}$ represent the set of real numbers and the
set of positive real numbers, respectively. The one-dimensional Gaussian
distribution with mean $\mu$ and variance $\sigma^{2}$ is denoted by ${\cal
N}(\mu,\sigma^{2})$. The multivariate Gaussian distribution with mean vector
$\bm{\mu}$ and covariance matrix $\bm{\Sigma}$ is represented by ${\cal
N}(\bm{\mu},\bm{\Sigma})$. Assume that a function
$f:\mathbb{R}\rightarrow\mathbb{R}$ is given. The same function $f$ can be
applied to $\bm{x}\equiv(x_{1},x_{2},\ldots,x_{n})$ as
$f(\bm{x})=(f(x_{1}),f(x_{2}),\ldots,f(x_{n}))$. The expectation operator is
denoted by ${\sf E}[\cdot]$. The notation $\mbox{diag}(\bm{x})$ denotes the
diagonal matrix whose diagonal elements are given by
$\bm{x}\in\mathbb{R}^{n}$. The matrix exponential
$\exp(\bm{X})(\bm{X}\in\mathbb{R}^{n\times n})$ is defined by

$\displaystyle\exp(\bm{X})\equiv\sum_{k=0}^{\infty}\frac{1}{k!}\bm{X}^{k}.$
(1)

The spectral norm of $\bm{X}\in\mathbb{R}^{n\times n}$ is denoted by
$\|\bm{X}\|_{2}$. The notation $[n]$ denotes the set of consecutive integers
from $1$ to $n$.

###  II-B System model

In this paper, we follow a common system model for compressed sensing [27,
28]. A sparse vector $\bm{s}\in\mathbb{R}^{n}$ is given, and it is assumed
that the original sparse signal follows the Bernoulli-Gaussian distribution,
where a non-zero element follows a Bernoulli distribution with probability $p$
and a non-zero element follows a Gaussian distribution ${\cal N}(0,1)$. A
linear observation vector $\bm{y}\in\mathbb{R}^{m}$ is generated by

$\displaystyle\bm{y}=\bm{A}\bm{s}+\bm{n},$ (2)

where $\bm{A}\in\mathbb{R}^{m\times n}$ is the sensing matrix. The vector
$\bm{n}\in\mathbb{R}^{m}$ is an additive Gaussian noise vector following
$\bm{n}\sim{\cal N}(\bm{0},\sigma^{2}\bm{I})$. The goal of the sparse signal
recovery problem is to estimate the original sparse vector $\bm{s}$ from the
knowledge of $\bm{y}$ and $\bm{A}$ as accurately as possible.

###  II-C Lasso formulation of sparse signal recovery

The Lasso objective function is defined as

$\displaystyle
f(\bm{x})\equiv\|\bm{y}-\bm{A}\bm{x}\|^{2}_{2}+\lambda\|\bm{x}\|_{1},$ (3)

where $\lambda(>0)$ is the regularization parameter. The L1 regularization
term is included in order to promote a sparse solution in the regularized LS
estimation. The optimization problem is then defined as

$\displaystyle\hat{\bm{x}}=\mbox{arg min}_{\bm{x}\in\mathbb{R}^{n}}\
f(\bm{x}).$ (4)

This optimization problem is convex and can be solved by any convex
optimization algorithm. A common approach to solving the Lasso problem is to
use a proximal gradient descent method, such as ISTA [35].

##  III ODE-Based Sparse Signal Recovery

###  III-A Gradient flow for sparse signal recovery

In the following argument, we introduce a gradient flow [41] based on the
Lasso objective function.

Gradient flow is a dynamical system of the form

$\displaystyle\frac{d\bm{x}(t)}{dt}=-\nabla{\cal E}(\bm{x}(t)),$ (5)

where ${\cal E}$ is an energy function. The value of the energy ${\cal
E}(\bm{x}(t))$ decreases monotonically as time $t$ increases. The behavior of
many physical systems can be described by an energy minimization process.

The Lasso objective function $f(\bm{x})$ in (3) is not differentiable at
$\bm{x}$ if $\bm{x}$ includes zero elements. This means that the gradient
$\nabla f(\bm{x})$ is not well-defined at some $\bm{x}$. Instead of the
objective function $f$ defined in (3), we here use the differentiable function
$g:\mathbb{R}^{n}\rightarrow\mathbb{R}$ defined by

$\displaystyle
g(\bm{x})\equiv\|\bm{y}-\bm{A}\bm{x}\|_{2}^{2}+\lambda\xi_{\alpha}(\bm{x})$
(6)

as an energy function in the gradient flow.

The function $\xi_{\alpha}:\mathbb{R}\rightarrow\mathbb{R}$ is given by

$\displaystyle\xi_{\alpha}(x)\equiv\frac{\log(\exp(\alpha x)+\exp(-\alpha
x))}{\alpha},$ (7)

where $\alpha\in\mathbb{R}_{+}$ is called the proximity parameter. The
function $\xi_{\alpha}(x)$ can be seen as a proxy function of the absolute
value function, i.e., $\xi_{\alpha}(x)\simeq|x|$ for sufficiently large
$\alpha$. At the limit $\alpha\rightarrow\infty$, the function $g$ converges
to the function $f$. The derivative function of $\xi_{\alpha}$ is given by

$\displaystyle\xi_{\alpha}^{\prime}(x)=\tanh(\alpha x).$ (8)

Figure 1 shows the shapes of $\xi_{\alpha}(x)$ and
$\xi_{\alpha}^{\prime}(x)=\tanh(\alpha x)$ for $\alpha=1,5,100$. It can be
observed that $\xi_{\alpha}(x)$ approaches the absolute value function $|x|$
as $\alpha$ increases.

Figure 1: Proxy function for $|x|$: Left: $\xi_{\alpha}(x)$, Right:
$\xi_{\alpha}^{\prime}(x)=\tanh(\alpha x)$

Note that we thus have the gradient of the energy function:

$\displaystyle\nabla
g(\bm{x})=\bm{A}^{T}(\bm{A}\bm{x}(t)-\bm{y})+\lambda\tanh(\alpha\bm{x}(t)).$
(9)

In this paper, we focus on the gradient flow defined by the following
nonlinear ODE:

$\displaystyle\frac{d\bm{x}(t)}{dt}$
$\displaystyle=-\left(\bm{A}^{T}(\bm{A}\bm{x}(t)-\bm{y})+\lambda\tanh(\alpha\bm{x}(t))\right),$
(10) $\displaystyle\bm{x}(0)$ $\displaystyle=\bm{x}_{0},$ (11)

where $\bm{x}(t)\in\mathbb{R}^{n}$. The variable $t$ represents continuous
time in our context. The vector $\bm{x}_{0}$ is the initial state providing a
boundary condition. Because the above ODE contains a nonlinear term
$\tanh(\alpha\bm{x}(t))$, it is not possible to obtain a succinct closed-form
solution. Therefore, in order to understand the dynamics of the system, we
must resort to using numerical methods to solve the differential equation.

###  III-B Numerical solution via Euler method

The Euler method is the simplest numerical method for solving simultaneous
nonlinear differential equations [16]. Although the convergence order of the
Euler method is inferior to that of higher-order methods such as the Runge-
Kutta methods [16], the Euler method is simple to use and can provide
sufficiently precise solutions if sufficient discretization of the time
interval is used. Thus, in this work, we use the Euler method to solve (10).

Consider the ODE in (10). Assume that we need an approximation of the
numerical solution of the above ODE in the time interval $0\leq t\leq T$. This
interval is first divided into $N$ bins. The discrete-time ticks
$t_{k}=k\eta(k=0,1,\ldots,N)$ define the boundaries of the bins, where the
width of a bin, $\eta$, is given by

$\displaystyle\eta\equiv\frac{T}{N}.$ (12)

It should be noted that the choice of the bin width $\eta$ is crucial in order
to ensure the stability and the accuracy of the Euler method. A small width
leads to a more accurate solution, but requires more computational time. A
large width may be computationally efficient but may lead to instability in
the solution. Let us define a discretized sample $\bm{x}^{(k)}$ as
$\bm{x}^{(k)}\equiv\bm{x}(t_{k})$.

By using the Euler method, the solution of (10) can be approximated by the
following recursive formula:

$\displaystyle\bm{x}^{(k+1)}$ $\displaystyle=\bm{x}^{(k)}-\eta\nabla
g(\bm{x}^{(k)})$ (13)
$\displaystyle=\bm{x}^{(k)}-\eta\left(\bm{A}^{T}(\bm{A}x^{(k)}-\bm{y})+\lambda\tanh(\alpha\bm{x}^{(k)})\right),$
(14) $\displaystyle\bm{x}^{(0)}$ $\displaystyle=\bm{x}_{0},$ (15)

for $k=0,1,2,\ldots,N$. The initial value is set to $\bm{x}^{(0)}=\bm{x}_{0}$.

###  III-C Numerical experiments

To confirm the behavior of the ODE-based sparse signal recovery (10),
numerical experiments were conducted. The parameter settings for the
experiments were as follows. The length of the sparse signal $\bm{s}$ was set
to $n=128$. The length of the observation vector $\bm{y}$ was $m=64$. The
standard deviation of the Gaussian noises was set to $\sigma=0.1$. The sensing
matrix $\bm{A}$ was randomly generated. Each element in $\bm{A}$ follows
${\cal N}(0,1)$. In the following part of this paper, the same $\bm{A}$ is
used for the experiments.

Figure 2 shows the state trajectories of $\bm{x}(t)$, which can be considered
as a solution of the ODE (10). Each component of $\bm{x}(t)$ corresponds a
curve in Fig. 2. The Euler method with $N=500$ was used to solve the ODE. The
observation vector $\bm{y}$ was generated randomly according to the system
model (2). We can see that the values of many components approach zero and
that only a faction of the curves corresponding to non-zero elements in
$\bm{s}$ deviate from zero. Namely, we can see that $\bm{x}(t)$ becomes a
sparse vector with time.

Figure 2: A solution of the ODE (10). Each component of $\bm{x}(t)$ is
depicted as a function of time $t$. The parameters are:
$n=128,m=64,p=0.1,\sigma=0.1,\lambda=5,\alpha=50$, $\bm{x}_{0}=\bm{0}$,
$T=0.5$.

Figure 3 shows sparse signal recovery processes based on ODE (10) with
regularization parameter $\lambda=5$ and proximity parameter $\alpha=50$. For
the numerical evaluation, the Euler method with $N=5000$ and $T=0.5$ was used
for each $T\in\\{0.01,0.05,0.1,1.0\\}$. The squared error
$\|\bm{s}-\bm{x}(T)\|_{2}^{2}$ is represented as SE in Fig. 3. At $T=0.01$,
the reconstructed vector appears to be a noise vector. However, we can observe
that the reconstructed vector approaches the original signal as $T$ increases.
These results confirm that the ODE-based sparse signal recovery can be
successful with appropriate parameter settings.

Figure 3: Example of the sparse signal recovery process. The horizontal axis
represents the index of a vector of length 128. The value of $x_{i}(T)$ is
plotted at the $i$th index position. The parameters are:
$n=128,m=64,p=0.1,\sigma=0.1,\lambda=5,\alpha=50$, $\bm{x}_{0}=\bm{0}$.

To confirm the accuracy of the Euler method, we conducted an experiment in
which we examined the relationship between the mean squared error (MSE) and
the number of bins $N$. We first generated a set of sparse signals
$\bm{s}_{1},\bm{s}_{2},\ldots,\bm{s}_{M}$ according to the Bernoulli-Gaussian
distribution with $p=0.1$ and non-zero elements following ${\cal N}(0,1)$. The
MSE at time $t$ is defined by ${\sf E}[\|\bm{x}(t)-\bm{s}\|_{2}^{2}]$, which
can be estimated by

$\displaystyle{\sf MSE}(t)\equiv{\sf
E}[\|\bm{x}(t)-\bm{s}\|_{2}^{2}]\simeq\frac{1}{M}\sum_{i\in[M]}\|\bm{x}^{(k)}_{i}-\bm{s}_{i}\|_{2}^{2},$
(16)

where $\bm{x}^{(k)}_{i}$ is the output of the Euler method for the $i$th trial
at index $k$ corresponding to time $t$. The parameter settings are nearly the
same as the settings in the previous experiment; that is,
$n=128,m=64,\sigma=0.1,\lambda=5,\alpha=50$. The number of samples was set to
$M=1000$.

Table I shows the estimates of MSE for several different discretizations.
Parameter $T$ for the discretizations was set to 4. From Table I, we can see
that the values of the MSE estimates converge as $N$ increases. This result
provides evidence that the Euler method works adequately for this ODE. In the
following experiments, we use $N=5000$ to ensure good accuracy and a
reasonable computational time.

TABLE I: Number of bins $N$ and estimated MSE ($T=4$) $N$ | ${\sf MSE}(T)$
---|---
1000 | 0.503796
2000 | 0.393782
5000 | 0.408937
10000 | 0.408701
Figure 4: Estimates of mean squared error ${\sf MSE}(t)$ as a function of time
$t$ ($n=128,m=64,p=0.1,\sigma=0.1,\alpha=50$).

Figure 4 plots the estimates of ${\sf MSE}(t)$ as a function of time $t$. The
MSE estimates in Fig. 4 are based on 100 trials of sparse signal recovery
processes, i.e., $M=100$. The Euler method with $N=5000$ was used for $T=4$ to
estimate ${\sf MSE}(t)(0\leq t\leq T)$. The parameter settings are nearly the
same as in the previous experiment. We examined three cases:
$\lambda=0.5,1.5,3$. From Fig. 4, it is apparent that the MSE curve for
$\lambda=0.5$ has the largest MSE in the range $0\leq t\leq 3$; however, the
MSE value continues to decrease after $t=3$. On the other hand, the MSE curve
for $\lambda=3$ shows the fastest convergence, but the value of ${\sf MSE}(t)$
saturates to a relatively high constant value in the range $t>1.0$. The MSE
curve for $\lambda=1.5$ provides a slightly slower convergence compared with
the case of $\lambda=3$, but it shows a lower floor for the MSE values. These
experimental results imply that the regularization constant $\lambda$ has a
strong influence on convergence behavior, similar to discrete time sparse
signal recovery algorithms such as ISTA.

##  IV Local Convergence Analysis

###  IV-A Linear approximation around equilibrium point

In the previous section, we saw that the convergence behavior of the ODE (10)
is highly dependent on the choice of the regularization parameter $\lambda$.
In this section, the mechanism by which the regularization parameter affects
the convergence behavior is examined. The core of the analysis presented here
is the use of a linear approximation around the equilibrium point [41] of the
ODE to study the local convergence behavior.

In the following analysis, we consider the objective function (6). The
gradient of $g$ is thus given by

$\displaystyle\nabla
g(\bm{x})=\bm{A}^{T}(\bm{A}\bm{x}-\bm{y})+\lambda\tanh(\alpha\bm{x}).$ (17)

We can now examine the local convergence behavior of the continuous-time
dynamical system defined by

$\displaystyle\frac{d\bm{x}(t)}{dt}$
$\displaystyle=-\left(\bm{A}^{T}(\bm{A}\bm{x}(t)-\bm{y})+\lambda\tanh(\alpha\bm{x}(t))\right),$
(18) $\displaystyle\bm{x}(0)$ $\displaystyle=\bm{x}_{0},$ (19)

where the initial vector $\bm{x}(0)=\bm{x}_{0}$ is assumed to be given as a
boundary condition.

The equilibrium point $\bm{x}^{*}$ of the above ODE is the point satisfying
$\nabla g(\bm{x})=\bm{0}$. Thus, the equilibrium point $\bm{x}^{*}$ satisfies
the equality

$\displaystyle\bm{A}^{T}(\bm{A}\bm{x}^{*}-\bm{y})+\lambda\tanh(\alpha\bm{x}^{*})=\bm{0},$
(20)

which is called equilibrium equality.

To establish the dynamics of the residual error, we change the coordinate of
the ODE according to

$\displaystyle\bm{x}(t)=\bm{x}^{*}+\bm{e}(t),$ (21)

where $\bm{e}(t)\equiv\bm{x}(t)-\bm{x}^{*}$ represents the residual error
vector. In the following analysis, $\bm{e}(t)$ is abbreviated to $\bm{e}$ for
simplicity. Due to the change in the coordinate, the right-hand side of ODE
(18) can be transformed to

$\displaystyle-$
$\displaystyle\left(\bm{A}^{T}(\bm{A}\bm{x}-\bm{y})+\lambda\tanh(\alpha\bm{x})\right)$
(22)
$\displaystyle=-\left(\bm{A}^{T}(\bm{A}(\bm{x}^{*}+\bm{e})-\bm{y})+\lambda\tanh(\alpha(\bm{x}^{*}+\bm{e}))\right)$
(23)
$\displaystyle=-\left(\bm{A}^{T}(\bm{A}\bm{x}^{*}-\bm{y})+\bm{A}^{T}\bm{A}\bm{e}+\lambda\tanh(\alpha(\bm{x}^{*}+\bm{e}))\right)$
(24)
$\displaystyle=-\left(\bm{A}^{T}\bm{A}\bm{e}+\lambda\tanh(\alpha(\bm{x}^{*}+\bm{e}))-\lambda\tanh(\alpha\bm{x}^{*})\right).$
(25)

The last equality is due to the equilibrium equality.

By expanding $\tanh(\alpha\bm{x})$ using the Taylor expansion and ignoring the
higher-order terms, we can obtain the linear approximation around the
equilibrium point:

$\displaystyle\tanh(\alpha\bm{x})\simeq\tanh(\alpha\bm{x}^{*})+\bm{J}(\bm{x}^{*})(\bm{x}-\bm{x}^{*}),$
(26)

where $\bm{J}(\bm{x}^{*})$ is the Jacobian matrix of $\tanh(\alpha\bm{x})$ at
$\bm{x}^{*}$. The Jacobian matrix is given by

$\displaystyle J(\bm{x}^{*})=\mbox{diag}\left(\frac{\alpha}{\cosh^{2}(\alpha
x_{1}^{*})},\frac{\alpha}{\cosh^{2}(\alpha
x_{2}^{*})},\ldots,\frac{\alpha}{\cosh^{2}(\alpha x_{n}^{*})}\right),$ (27)

where $\bm{x}^{*}=(x_{1}^{*},x_{2}^{*},\ldots,x_{n}^{*})^{T}$. Note that
$\cosh^{2}(x)>0$ for any $x\in\mathbb{R}$. This implies that $J(\bm{x}^{*})$
is positive definite for any $\bm{x}^{*}\in\mathbb{R}^{n}$.

By using the linear approximation, we immediately have the approximation:

$\displaystyle\lambda\tanh(\alpha(\bm{x}^{*}+\bm{e}))-\lambda\tanh(\alpha\bm{x}^{*})$
$\displaystyle\simeq\lambda J(\bm{x}^{*})\bm{e}$ (28)

In the following argument, we will treat this approximated equality as an
equality for the sake of simplicity. Substituting this equation into (25), we
obtain the linear approximation of the right-hand side of (18):

$\displaystyle-\left(\bm{A}^{T}(\bm{A}\bm{x}-\bm{y})+\lambda\tanh(\alpha\bm{x})\right)$
$\displaystyle=-\left(\bm{A}^{T}\bm{A}+\lambda\bm{J}(\bm{x}^{*})\right)\bm{e}$
(29)

Finally, we have the ODE representing the evolution of the residual error
vector:

$\displaystyle\frac{d\bm{e}(t)}{dt}$
$\displaystyle=-\left(\bm{A}^{T}\bm{A}+\lambda\bm{J}(\bm{x}^{*})\right)\bm{e}(t),$
(30) $\displaystyle\bm{e}(0)$ $\displaystyle=\bm{e}_{0},$ (31)

where $\bm{e}_{0}$ represents the initial error vector. This is a linear ODE,
which can be solved in a concise form [16]:

$\displaystyle\bm{e}(t)=\exp\left(-\left(\bm{A}^{T}\bm{A}+\lambda\bm{J}(\bm{x}^{*})\right)t\right)\bm{e}(0).$
(32)

We can expect the above solution to reflect the actual error behavior well if
the initial error $\bm{e}(0)$ is sufficiently close to the zero vector, i.e.,
in the case where the linear approximation is appropriate.

From the error evolution equation (32), it is clear that the eigenvalues of
the matrix $\left(\bm{A}^{T}\bm{A}+\lambda\bm{J}(\bm{x}^{*})\right)t$ dominate
the behavior of the error evolution. Due to the positive definiteness of the
Jacobian matrix $\bm{J}(\bm{x}^{*})$ and the semi-positive definiteness of the
Gram matrix $\bm{A}^{T}\bm{A}$, all the eigenvalues of the matrix
$\left(\bm{A}^{T}\bm{A}+\lambda\bm{J}(\bm{x}^{*})\right)t$ are positive real
numbers. This means that the time evolution of error vector (32) is locally
asymptotically stable; that is, $\lim_{t\rightarrow\infty}\bm{e}(t)=\bm{0}$
for any initial error $\bm{e}(0)$ sufficiently close to the zero vector.

###  IV-B Upper bound on convergence speed

From the time evolution or error vector (32), we immediately have the
following theorem, under the assumption that the linear approximation is
appropriate.

###### Theorem 1

If the initial vector $\bm{x}(0)$ is sufficiently close to $\bm{x}^{*}$, then
the following inequality holds:

$\displaystyle\frac{\|\bm{e}(t)\|_{2}}{\|\bm{e}(0)\|_{2}}$
$\displaystyle\leq\exp(-\omega_{1}t),$ (33)

where $\\{\omega_{i}\\}_{i=1}^{n}$ is the set of eigenvalues of
$\bm{A}^{T}\bm{A}+\lambda\bm{J}(\bm{x}^{*}),$ with the order
$0<\omega_{1}\leq\omega_{2}\leq\cdots\leq\omega_{n}$.

(Proof) Taking the norm of both sides of (32), we have an inequality:

$\displaystyle\|\bm{e}(t)\|_{2}$
$\displaystyle=\|\exp\left(-\left(\bm{A}^{T}\bm{A}+\lambda\bm{J}(\bm{x}^{*})\right)t\right)\bm{e}(0)\|_{2}$
(34)
$\displaystyle\leq\|\exp\left(-\left(\bm{A}^{T}\bm{A}+\lambda\bm{J}(\bm{x}^{*})\right)t\right)\|_{2}\|\bm{e}(0)\|_{2},$
(35)

where the matrix norm is the spectral norm. By dividing both sides by
$\|\bm{e}(0)\|_{2}$, we obtain

$\displaystyle\frac{\|\bm{e}(t)\|_{2}}{\|\bm{e}(0)\|_{2}}$
$\displaystyle\leq\|\exp\left(-\left(\bm{A}^{T}\bm{A}+\lambda\bm{J}(\bm{x}^{*})\right)t\right)\|_{2}.$
(36)

Note that the above matrix exponential can be rewritten as

$\displaystyle\exp\left(-\left(\bm{A}^{T}\bm{A}+\lambda\bm{J}(\bm{x}^{*})\right)t\right)=\bm{U}\mbox{diag}(e^{-\omega_{1}t},\ldots,e^{-\omega_{n}t})\bm{U}^{T},$
(37)

where $\bm{U}$ is an orthogonal matrix. Since the spectral norm of a semi-
positive definite symmetric matrix $\bm{X}$ coincides with the largest
eigenvalue of $\bm{X}$, we have

$\displaystyle\|\exp\left(-\left(\bm{A}^{T}\bm{A}+\lambda\bm{J}(\bm{x}^{*})\right)t\right)\|_{2}=\exp(-\omega_{1}t).$
(38)

Substituting this equality into $(\ref{fund_ineq})$, we obtain the claim of
the theorem.

The above theorem indicates that the term $\exp(-\omega_{1}t)$ dominates the
convergence speed around an equilibrium point. The average behavior is
described by the following corollary, which can be derived directly from
Theorem 1.

###### Corollary 1

If the initial vector $\bm{x}(0)$ is sufficiently close to $\bm{x}^{*}$, then
the following inequality holds:

$\displaystyle{\sf E}\left[\frac{\|\bm{e}(t)\|_{2}}{\|\bm{e}(0)\|_{2}}\right]$
$\displaystyle\leq{\sf E}[\exp(-\omega_{1}t)].$ (39)

Let $l_{min}(\bm{X})$ represent the minimum eigenvalue of a symmetric matrix
$\bm{X}$. The minimum eigenvalue of
$\bm{A}^{T}\bm{A}+\lambda\bm{J}(\bm{x}^{*})$ satisfies

$\displaystyle l_{\sf min}(\bm{A}^{T}\bm{A}+\lambda\bm{J}(\bm{x}^{*}))$
$\displaystyle\geq l_{\sf min}(\bm{A}^{T}\bm{A})+l_{\sf
min}(\lambda\bm{J}(\bm{x}^{*}))$ (40) $\displaystyle=l_{\sf
min}(\bm{A}^{T}\bm{A})+\lambda\min_{i\in[n]}\frac{\alpha}{\cosh^{2}(\alpha
x_{i}^{*})}.$ (41)

Therefore, increasing the value of $\lambda$ leads to larger $\omega_{1}$,
which increases the convergence speed.

###  IV-C MSE behavior

The value on the right-hand side of equation (39) can be reduced by choosing a
larger $\lambda$. This implies that a larger $\lambda$ will achieve faster
convergence to the equilibrium point. This qualitative argument explains the
numerical results shown in Fig. 4; in particular, a larger $\lambda$ results
in a faster decrease of ${\sf MSE}(t)$.

On the other hand, from the equilibrium equation (20), we have

$\displaystyle\|\bm{A}^{T}(\bm{A}\bm{x}^{*}-\bm{y})\|_{2}=\lambda\|\tanh(\alpha\bm{x}^{*})\|_{2}.$
(42)

Consider a scenario where the noise is negligible. Recall that the original
sparse vector $\bm{s}$ should satisfy

$\displaystyle\|\bm{A}^{T}(\bm{A}\bm{s}-\bm{y})\|_{2}\simeq 0$ (43)

under the assumption of negligible noise. In such a case, an increase in
$\lambda$ can lead to an increase in the discrepancy
$\|\bm{x}^{*}-\bm{s}\|_{2}$. From this observation, it can be concluded that a
larger $\lambda$ leads to higher MSE floors, which is indeed observed in Fig.
4.

###  IV-D Numerical experiments

The analysis presented in the previous sections is based on a linear
approximation around the equilibrium point. Therefore, the accuracy of this
approximation should be verified by numerical experiments.

The error norm ratio $\rho(t)$ is defined by

$\displaystyle\rho(t)\equiv\frac{\|\bm{e}(t)\|_{2}}{\|\bm{e}(0)\|_{2}}=\frac{\|\bm{x}(t)-\bm{x}^{*}\|_{2}}{\|\bm{x}(0)-\bm{x}^{*}\|_{2}}.$
(44)

In the following experiment, we will evaluate the norm ratio $\rho(t)$ for the
two cases $\lambda=1.5$ and $\lambda=5$.

The parameter settings are nearly the same as in the previous experiment:
$n=128,m=64,\sigma=0.1,\alpha=50$. The equilibrium point $\bm{x}^{*}$ is
approximated by the value of $\bm{x}(T)$ where $T=4$.

We consider two types of initial points, $\bm{x}(0)=\bm{0}$ and
$\bm{x}(0)=\hat{\bm{x}}\equiv\bm{x}^{*}+\bm{\epsilon}$, where
$\bm{\epsilon}\sim{\cal N}(\bm{0},0.1^{2}\bm{I})$. Note that
$\bm{x}(0)=\bm{0}$ may not be in the region where the linear approximation is
valid but $\hat{\bm{x}}$ can be in such a region because the norm of
$\bm{\epsilon}$ is so small.

Figure 5: Error norm ratio $\rho(t)$ as a function of time $t$. Left
$\lambda=1.5$, Right $\lambda=5$ ($n=128,m=64,p=0.1,\sigma=0.1,\alpha=50$).

Figure 5 displays the error norm ratio $\rho(t)$ as a function of time $t$.
The solid line labeled “Theory” corresponds to the value of
$\exp(-\omega_{1}t)$. The left panel shows the case where $\lambda=1.5$. The
error norm ratio curve of $\rho(t)$ with the initial condition
$\bm{x}(0)=\bm{0}$ has a shallower slope than that of the curve of
$\exp(-\omega_{1}t)$ when $t<0.7$. This is due to the invalidity of the linear
approximation for small $t$. It can be observed that the curve of $\rho(t)$
with $\bm{x}_{0}=\bm{0}$ has almost the same slope as that of
$\exp(-\omega_{1}t)$ when $t>0.7$. This can be explained by the fact that the
state point $\bm{x}(t)$ gradually enters the region where the linear
approximation is valid. Once the state point enters such a region, the linear
approximation theory can be applied. On the other hand, in the case of
$\bm{x}(0)=\hat{\bm{x}}$, the slope of $\rho(t)$ is almost the same as that of
$\exp(-\omega_{1}t)$ from the beginning. The result is consistent with the
claim of Theorem 1. Since $\hat{\bm{x}}$ is sufficiently close to
$\bm{x}^{*}$, the linear approximation is valid from the start.

The right panel shows the numerical results for $\lambda=5$. Nearly the same
observations can be made as in the left panel. Compared to the left panel, the
value of $\rho(t)$ is much smaller for each $t$. This observation is
consistent with the argument concerning the relationship between $\lambda$ and
the convergence rate in the previous section.

Note that in all the cases, the asymptotic slope of $\rho(t)$ can be
accurately estimated by the slope of $\exp(-\omega_{1}t)$. This observation
supports our claim that the smallest eigenvalue $\omega_{1}$ determines the
convergence rate around an equilibrium point. The numerical results presented
in Fig. 5 can be seen as a numerical validation of the argument in the
previous sections.

Figure 6 shows ${\sf MSE}(\infty)$, ${\sf
E}[|\lambda\tanh(\alpha\bm{x}^{*})\|_{2}^{2}]$, and ${\sf E}[\omega_{1}]$
estimated from 500 trials. For estimating the asymptotic MSE ${\sf
MSE}(\infty)$, we used $T=4$ rather than $T=\infty$. The numerical results of
${\sf MSE}(\infty)$ indicate that the MSE floor increases as $\lambda$
increases. We can also observe that the values ${\sf
E}[|\lambda\tanh(\alpha\bm{x}^{*})\|_{2}^{2}]$ (center panel) is an increasing
function of $\lambda$. This observation is consistent with the argument in
Subsection IV-C. The local convergence rate is determined by $\omega_{1}$. The
tendency of ${\sf E}[\omega_{1}]$ is shown in the third panel of Fig. 6. It
can be confirmed that a larger $\lambda$ achieves faster convergence.

Figure 6: Values of ${\sf MSE}(\infty)$, ${\sf
E}[|\lambda\tanh(\alpha\bm{x}^{*})\|_{2}^{2}]$, and ${\sf E}[\omega_{1}]$ from
left to right ($n=128,m=64,p=0.1,\sigma=0.1,\alpha=50$).

Finally, we will experimentally confirm the relationship between $\alpha$ and
the MSE convergence behavior. Figure 7 shows the values of ${\sf MSE}(\infty)$
and ${\sf E}[\omega_{1}]$ where the value of $\lambda$ is fixed to 3. From the
results of ${\sf MSE}(\infty)$, we can say that a larger $\alpha$ results in a
smaller MSE floor. The right panel showing ${\sf E}[\omega_{1}]$ indicates
that a larger $\alpha$ leads to faster convergence. In summary, $\alpha$
should be large in order to achieve faster convergence and a high-quality
solution. We therefore use a relatively large constant $\alpha=50$ throughout
this paper.

Figure 7: Values of ${\sf MSE}(\infty)$ and ${\sf E}[\omega_{1}]$
($n=128,m=64,p=0.1,\sigma=0.1,\lambda=3$).

##  V Deep Unfolded-Variational Optimization

###  V-A Parametric ODE and its optimization

A summary of the theoretical analysis and numerical results presented in the
previous section is illustrated in Fig. 8, which shows the typical behavior of
the convergence of ${\sf MSE}(t)={\sf E}[\|\bm{x}(t)-\bm{s}\|_{2}^{2}]$. As
discussed earlier, there is a trade-off between the MSE floor value (the
asymptotic MSE, ${\sf MSE}(\infty)$) and the convergence rate. If we use a
large $\lambda$, we can expect fast convergence, but we will need to
compromise on the quality of the solution. On the other hand, if we use a
small $\lambda$, we must allow for slow convergence.

An important consideration is the time period during which sub-linear
convergence occurs. If the state point $\bm{x}(t)$ is close enough to the
equilibrium point, the linear approximation can be applied and linear
convergence can be expected, as discussed in a previous section. However, if
$\bm{x}(t)$ is not close enough to the equilibrium point, linear convergence
cannot be expected. The left panel of Fig. 5 illustrates such a regime,
specifically for $0\leq t\leq 0.6$. As shown in Fig. 8, the linear convergence
regime appears after a sub-linear regime.

Figure 8: Typical tendency of the convergence behavior of ${\sf MSE}(t)$

To improve both convergence speed and solution quality, one strategy is to
introduce a time-dependent regularization parameter, i.e., varying the value
of $\lambda$. Starting with a large $\lambda$ shortens the sub-linear regime
and increases the convergence rate. Changing $\lambda$ judiciously keeps
$\bm{x}(t)$ in the linear convergence regime. By adjusting $\lambda$ to take
on the smallest value at the end of the process, a high-quality solution can
be obtained.

This idea is closely related to the homotopy or continuation method, where the
value of the regularization parameter is modulated so that the state vector
progresses along the solution path. For example, Xiao and Zhang [48] proposed
a technique for solving the Lasso problem by using a proximal algorithm, which
incorporates a sequence of decreasing values of the regularization parameter.
The concept of the continuation method is shown in Fig. 9.

Figure 9: Concept of continuation method: $x^{*}(\lambda_{i})(i\in[5])$
represents the equilibrium point when $\lambda=\lambda_{i}$. Using an
appropriately chosen set of parameters $\\{\lambda_{i}\\}$, the solution path
is fully covered by the linear convergence regions. We can anticipate both
faster convergence and a high-quality solution in such a case.

The above idea naturally leads to the following parametric ODE:

$\displaystyle\frac{d\bm{x}(t)}{dt}$
$\displaystyle=-\left(\bm{A}^{T}(\bm{A}\bm{x}(t)-\bm{y})+\lambda(t)\tanh(\alpha\bm{x}(t))\right),$
(45) $\displaystyle\bm{x}(0)$ $\displaystyle=\bm{x}_{0},$ (46)

where the regularizing constant is replaced with the function
$\lambda:\mathbb{R}\rightarrow\mathbb{R}$.

In this section, we examine the optimization of $\lambda(t)$ in this
parametric ODE. The optimization problem at hand can be formulated as

$\displaystyle\mbox{minimize }{\sf E}[\|\bm{x}(T^{*})-\bm{s}\|_{2}^{2}]\
\mbox{subject to }\lambda:\mathbb{R}\rightarrow\mathbb{R},$ (47)

where $T^{*}$ is the predetermined target time. In essence, we aim to find the
optimal continuous-time schedule for $\lambda$ that produces a high-quality
solution at a given time $T^{*}$. This optimization problem falls into the
category of variational optimization problems. Unfortunately, we cannot expect
to derive a compact analytical solution to this problem, and, thus, it must be
solved numerically.

###  V-B Deep unfolded-variational optimization (DU-VO) method

The development of deep neural networks has also had a significant impact on
the design of algorithms for communications and signal processing [17, 18,
19]. Deep unfolding, as described in works such as [22, 25, 26], can be seen
as a highly effective method for improving the convergence of iterative
algorithms. Gregor and LeCun introduced the Learned ISTA (LISTA) [22], which
uses learnable matrices. LISTA achieves a recovery performance that is much
superior to that of the original ISTA. Borgerding et al. also presented
variants of AMP and VAMP with learnable capability [23][24]. Trainable ISTA
(TISTA) [26] is a recent learnable sparse signal recovery algorithm with fast
convergence. TISTA requires a small number of trainable parameters, which
provides a fast and stable training process.

The concept of deep unfolding is simple: Embed trainable parameters in the
original iterative algorithm, followed by the unfolding of the signal-flow
graph of the original algorithm. The standard supervised training techniques
used in deep learning, such as stochastic gradient descent (SGD) [20] and back
propagation [21], can then be applied to the unfolded signal-flow graph to
optimize the trainable parameters.

In the following, we will introduce the deep unfolded-variational optimization
(DU-VO) method for the numerical solution of variational optimization
problems. The combination of deep unfolding and the Euler method for
differential equation solvers [44] is a current area of active research in
scientific machine learning. However, it should be noted that the technique is
not limited to applications within scientific machine learning; rather, it
seems to be a powerful tool for optimizing continuous-time dynamical systems.

This subsection deals with variational optimization problems in a general
context. We will return to the specific optimization problem defined in (47)
later.

Assume that we have an ODE

$\displaystyle\frac{d\bm{x}(t)}{dt}=h(\bm{x}(t),u(t),t)$ (48)

and it defines the dynamics of $\bm{x}(t)$. For a given integrable function
$F$, a functional $J$ to be minimized is defined by

$\displaystyle J[u(t)]\equiv\int_{0}^{T}F(\bm{x}(t),u(t),t)dt.$ (49)

In the following argument, we assume the case where $u(t)$ is one-dimensional,
i.e., $u:\mathbb{R}\rightarrow\mathbb{R}$, and can control this function. We
here consider a variational optimization problem including $\bm{x}(t)$ defined
by an ODE (48):

$\displaystyle\mbox{minimize }_{u(t)}\ J[u(t)].$ (50)

Namely, we want to find a one-dimensional function $u$ minimizing the value of
the functional $J$.

This type of variational optimization problem often appears in the field of
optimal control [46]. The optimal solution of (50) can be obtained by solving
a Hamilton-Jacobi-Bellman (HJB) equation [46]. However, solving the HJB
equation is generally difficult because it requires solving a non-linear
partial differential equation. In a stochastic context, solving a stochastic
HJB equation [46] becomes even more complex and challenging.

The DU-VO method is a numerical approach that aims to solve the variational
problem (50). The first step in the implementation of the DU-VO method is to
introduce a function approximation using a radial basis function (RBF), which
is used to approximate $u(t)$.

Figure 10: Set of Gaussian RBFs: $\phi(t-c_{i})=\exp(-\beta(t-c_{i})^{2})$.
$\beta=20$, $c_{1}=-0.5,c_{2}=-0.45,\ldots,c_{m}=2$.

Let $\phi:\mathbb{R}\rightarrow\mathbb{R}$ be an RBF, which satisfies
$\phi(x)=\phi(|x|).$ A continuous function $u(t)$ can be approximated by

$\displaystyle\tilde{u}(t)\equiv\sum_{i=1}^{S}w_{i}\phi(t-c_{i}),$ (51)

where the weight vector

$\bm{w}\equiv(w_{1},w_{2},\ldots,w_{S})^{T}\in\mathbb{R}^{S}$

is a trainable parameter that can be updated in an optimization process. The
above function approximation is known as the RBF approximation. The shift
parameter

$\displaystyle\bm{c}\equiv(c_{1},c_{2},\ldots,c_{S})^{T}\in\mathbb{R}^{S}$
(52)

is treated as a hyperparameter. Thus, we will use the parametric model of
$u(t)$ defined by (51) in the following discussion. Figure 10 shows a set of
Gaussian RBFs defined by $\phi(t-c_{i})=\exp(-\beta(t-c_{i})^{2}).$ We use the
Gaussian RBF in the numerical computation below.

The advantage of using the RBF is that it requires only a small number of
trainable parameters to approximate a one-dimensional continuous function,
resulting in a stable and efficient optimization process.

As in Subsection III-B, the closed interval $\\{t\mid 0\leq t\leq T\\}$ is
divided into $N$ bins, and the discrete time indices are defined as
$t^{(k)}\equiv\eta k(k=0,1,\ldots,N)$, where $\eta\equiv T/N$ is the width of
the bins. We also introduce the discretized state vector
$\bm{x}^{(k)}\in\mathbb{R}^{n}$ corresponding to $\bm{x}(t^{(k)})$.

Algorithm 1 DU-VO method

0: $\bm{x}_{0}$ (initial state), $\bm{w}_{0}$ (initial weight)

0: $\bm{w}$ (optimized weight)

1: Set the initial values: $\bm{x}^{(0)}\equiv\bm{x}_{0}$, $J^{(0)}\equiv 0$,
$\bm{w}=\bm{w}_{0}$.

2: for $i=1$ to $I$ do

3: for $i=0$ to $N-1$ do

4: Compute the Euler step:  $\bm{x}^{(k+1)}\equiv\bm{x}^{(k)}+\eta
h(\bm{x}^{(k)},\tilde{u}(t^{(k)}),t^{(k)}).$

5: Compute the numerical integration step:  $J^{(k+1)}\equiv J^{(k)}+\eta
F(\bm{x}^{(k)},\tilde{u}(t^{(k)}),t^{(k)}).$

6: end for

7: Compute the gradient of the loss function by using back propagation:
$\displaystyle\bm{g}\equiv\nabla J(\bm{w}).$

8: The weight parameter $\bm{w}$ is updated by using $\bm{g}$ (an SGD
optimizer is used).

9: end for

10: Output the weight vector $\bm{w}$.

In order to approximate the solution of (48), we have a number of numerical
methods to choose from, such as Runge-Kutta methods; however, for simplicity,
the simplest of these, the Euler method, will be used in the following
argument. Similarly, for the numerical integration of the functional (49), we
will exploit the simplest rectangular method. If necessary, the accuracy of
the following DU-VO method can be improved by replacing the ODE solver
algorithm and the numerical integration method.

The core of the DU-VO method is the Euler method for numerically solving
differential equation (48),

$\displaystyle\bm{x}^{(k+1)}=\bm{x}^{(k)}+\eta
h(\bm{x}^{(k)},\tilde{u}(t^{(k)}),t^{(k)}),\quad k=0,1,\ldots,N-1,$ (53)

and the numerical integration step for the functional (49),

$\displaystyle J^{(k+1)}=J^{(k)}+\eta
F(\bm{x}^{(k)},\tilde{u}(t^{(k)}),t^{(k)}),\quad k=0,1,\ldots,N-1.$ (54)

If the width of the bins is small enough, we can expect that $J^{(N)}$ will
become an approximation of $J[\tilde{u}(t)]$. The concept of deep unfolding
can be naturally applied to (53) and (54) to optimize the weight vector
$\bm{w}$ that controls the shape of $u$.

The DU-VO method uses a loss function $J(\bm{w})\equiv J^{(N)}$ to approximate
the solution of the variational problem (50). The gradient of the loss
function with respect to the weight vector $\bm{w}$, $\nabla J(\bm{w})$, is
evaluated by back propagation. The gradient is then used to update the weight
vector $\bm{w}$. This process can be optimized using common optimizers such as
Stochastic Gradient Descent (SGD) or Adam. The entire process is summarized in
Algorithm 1.

Even if the functional to be optimized contains random variables, the DU-VO
method can minimize the expected functional ${\sf E}[J(\bm{w})]$ by using
mini-batch training. This is the one of major advantages of the DU-VO method.

###  V-C Numerical example

To illustrate the application of the DU-VO method, we present a small
numerical example of an optimal control problem.

The variational optimization problem is described below. The state evolution
equation is given by

$\displaystyle\frac{dx(t)}{dt}=u(t).$ (55)

The functional is defined by

$\displaystyle J[u(t)]=\int_{0}^{1}(x(t)^{2}+u(t)^{2})dt,$ (56)

where $u:\mathbb{R}\rightarrow\mathbb{R}$. The boundary condition is $x(0)=1$,
and $x(1)$ can take any value. The exact solution of the problem is known
[47]:

$\displaystyle u^{*}(t)=-\frac{\sinh(1-t)}{\cosh(1)}.$ (57)

Figure 11: Optimal input $u(t)$: Comparison between the output of the DU-VO
method and the exact solution.

Figure 11 shows the output of the DU-VO method and the exact solution. The
parameter settings for the DU-VO method were as follows. The time
discretization with $N=200,T=1$ was used. The number of optimization loops $I$
was set to 200. The Adam optimizer with learning rate $0.1$ was used. The
parameter settings for the RBF approximation were as follows:
$|c_{i}-c_{i+1}|=0.05,\beta=20,S=50$. All the components of the weight vector
$\bm{w}$ were initialized to 1.

It can be seen here that there is no obvious visible difference between the
numerical results and the exact solution. This example thus supports the
validity of the DU-VO method.

###  V-D DU-VO method for ODE-based sparse signal recovery

In the following, we consider the parametric ODE (45) and the optimization
problem (47) once again. In order to optimize the convergence behavior with
respect to $\lambda(t)$, we need a functional of $\lambda$ to be minimized
that is consistent with the optimization problem (47). We here introduce a
functional $J$ defined by

$\displaystyle J[\lambda(t)]\equiv{\sf
E}\left[\int_{0}^{T^{*}}\delta(t-T^{*})\|\bm{x}(t)-\bm{s}\|^{2}dt\right],$
(58)

where $\bm{s}$ represents the original sparse vector and $\bm{x}(t)$ in the
functional is the solution of the ODE (45). The function $\delta$ represents
the Dirac’s delta function. The functional (58) measures closeness of the
solution $\bm{x}(T^{*})$ to the original vector. The time $T^{*}$ is the
predetermined target time.

The Euler recursive equation can be obtained by discretizing the ODE (45):

$\displaystyle\bm{x}^{(k+1)}$
$\displaystyle=\bm{x}^{(k)}-\eta\left(\bm{A}^{T}(\bm{A}\bm{x}^{(k)}-\bm{y})+\lambda(\eta
k)\tanh(\alpha\bm{x}^{(k)})\right),$ (59) $\displaystyle\bm{x}^{(0)}$
$\displaystyle=\bm{0},$ (60)

where $\eta=T^{*}/N$ is the width of the discretized bins. The initial vector
$\bm{x}^{(0)}$ is assumed to be the zero vector. This recursive equation
becomes the basis of the following DU-based optimization. The function
$\lambda:\mathbb{R}\rightarrow\mathbb{R}$ is replaced with a trainable
Gaussian RBF approximation [45] defined by

$\displaystyle\lambda(t)\equiv\sum_{i=1}^{S}w_{i}\exp(-\beta(t-\Delta
i+\theta)^{2}),$ (61)

where $\\{w_{i}\\}_{i=1}^{S}$ are the trainable weight parameters. The
parameters $\beta$, $\Delta$, and $\theta$ are treated as hyperparameters
given before an optimization process.

In an optimization process, we randomly generate multiple mini-batches. A
mini-batch is composed of $K$ pairs of vectors, namely

$\displaystyle{\cal
D}\equiv\\{(\bm{s}_{1},\bm{y}_{1}),\ldots,(\bm{s}_{K},\bm{y}_{K})\\},$ (62)

where $\bm{s}_{i}(i\in[K])$ is a sparse random vector following the Bernoulli-
Gaussian distribution, and $\bm{y}_{i}(i\in[K])$ follows our system model
$\bm{y}_{i}=\bm{A}\bm{s}_{i}+\bm{n}_{i}.$ The noise vector $\bm{n}_{i}$
follows ${\cal N}(\bm{0},\sigma^{2}\bm{I})$.

The functional $J$ can be approximated as

$\displaystyle J[\lambda(t)]$ $\displaystyle={\sf
E}\left[\int_{0}^{T^{*}}\delta(t-T^{*})\|\bm{x}(t)-\bm{x}^{*}\|^{2}dt\right]$
$\displaystyle\simeq\frac{1}{K}\sum_{i=1}^{K}\|\bm{x}^{(N)}_{i}-\bm{s}_{i}\|^{2},$
(63)

where the vector $\bm{x}_{i}^{(N)}$ is the state vector corresponding to the
data $(\bm{s}_{i},\bm{y}_{i})$. The state vector $\bm{x}_{i}^{(k)}$ is evolved
according to the Euler recursive equation (59) with the initial condition
$\bm{x}_{i}(0)=\bm{0}$.

From this approximation above, it is reasonable to define the loss function:

$\displaystyle{\sf Loss}({\cal
D})\equiv\frac{1}{K}\sum_{i=1}^{K}\|\bm{x}^{(N)}_{i}-\bm{s}_{i}\|^{2},$ (64)

which gives an approximation of MSE:

$\displaystyle{\sf Loss}({\cal D})\simeq{\sf E}[{\sf MSE}(T^{*})].$ (65)

For each given mini-batch, the trainable parameters $\\{w_{i}\\}_{i=1}^{S}$
are updated based on the gradient of ${\sf Loss}({\cal D})$. The optimization
process is illustrated in Fig. 12

Figure 12: Block diagram of Deep Unfolded-Variational Optimization (DU-VO)
method for optimizing the parametric ODE (45).

###  V-E Numerical experiments

Figure 13: Left: MSE as a function of $t$. Right: optimized $\lambda(t)$.
Target time is $T^{*}=3$. ($n=128,m=64,p=0.1,\sigma=0.1,\alpha=50$).

In the previous sections, we introduced the DU-based optimization method for
solving variational optimization problems. This section presents the results
of numerical experiments conducted to evaluate the performance of the proposed
method. These experiments were performed using the automatic differentiation
mechanism provided by the Flux.jl library [42] in the Julia programming
language [43].

The problem setup for the first experiment is similar to the previous ones,
with $n=128,m=64,p=0.1,\sigma=0.1$. The Euler method with $N=5000$ was used
for the optimization process and for estimating the MSE. The following
parameters were used for the optimization process. The mini-batch size was set
to $K=10$. The Adam optimizer with a learning rate of $10^{-2}$ was used. The
number of mini-batches, or iterations, used for training was $100$. The
parameter settings for the RBF approximation of $\lambda(t)$ were
$\Delta=0.25,\beta=20,S=20,\theta=0.5$.

The left panel of Fig. 13 shows the MSE curve of the parametric ODE (45) with
the optimized $\lambda(t)$. The MSE curves corresponding to constant
regularization parameters $\lambda=0.5,1.5,3.0$ are also shown. As can be seen
here, the optimized $\lambda(t)$ attains the smallest MSE value among the four
curves at $T^{*}=3$. Additionally, the convergence rate, i.e., the slope of
the optimized MSE curve, is nearly identical to that of the MSE curve with
$\lambda=0.5$ over the range $2\leq t\leq 3$. The right panel of Fig. 13 shows
the shape of the curve for the optimized $\lambda(t)$. The curve starts at a
value of approximately 2.2 and rapidly decreases over the range of $1\leq
t\leq 2$. This result supports the argument in Subsection V-A and suggests
that the curve of the optimized $\lambda(t)$ provides an appropriate schedule
for decreasing the value of the regularization parameter.

To confirm that the optimized $\lambda(t)$ provides the best MSE among the ODE
systems with constant $\lambda$, the MSE at $T^{*}=3$ was evaluated
numerically for several values of $\lambda$. The results are plotted in Fig.
14. As shown, the MSE values for the constant $\lambda$ are above $10^{-1}$,
while the optimized $\lambda(t)$ achieves an MSE below $10^{-1}$. This result
thus demonstrates the effectiveness of the proposed method in finding an
optimal $\lambda(t)$ that improves the MSE.

Figure 14: MSE at $T^{*}=3$ as a function of $\lambda$
($n=128,m=64,p=0.1,\sigma=0.1,\alpha=50$).

We now examine another problem setup that assumes a sparser original signal.
Here, the probability is set to $p=0.05$. The parameters and hyperparameters
remain unchanged from the previous experiments, i.e., $n=128,m=64,\sigma=0.1$.
The target time is set to $T^{*}=2$. The numerical results of this case are
shown in Fig. 15. It is clear that the ODE system with the optimized
$\lambda(t)$ consistently gives the lowest MSE. It is also noteworthy that the
optimized system consistently achieves close to the optimal value in the range
of $0\leq t\leq 2$. Both the convergence rate and the quality of the solution
are improved by the optimization.

In summary, for both cases, the proposed DU-VO method successfully finds an
optimized schedule for $\lambda(t)$ that improves the convergence speed and
the quality of the solution.

Figure 15: Left: MSE as a function of $t$, Right: Optimized $\lambda(t)$.
Target time is $T^{*}=2$. ($n=128,m=64,p=0.05,\sigma=0.1,\alpha=50$).

##  VI Conclusion

In this study, we presented a novel approach for solving sparse signal
recovery problems using continuous-time dynamical systems. The method is based
on a nonlinear ODE corresponding to the gradient flow of the Lasso objective
function. We presented the local convergence analysis of the proposed ODE
using a linear approximation around the equilibrium point. In addition, we
introduced a variational optimization problem to optimize the regularization
schedule and applied deep unfolding techniques [22, 25] to solve this problem.
To the best of our knowledge, this is the first work to propose the use of
deep unfolding techniques to solve variational optimization problems.

It should be noted that the methodology proposed in this paper can be applied
to various forms of regularized LS problems such as binary quadratic
minimization problems. Moreover, with the future advances in modern analog
computing, especially in the optical domain, it is likely that ODE-based
signal processing will become increasingly viable for achieving fast and
energy-efficient signal processing.

## Acknowledgment

This work was supported by JSPS KAKENHI Grant Number JP22H00514.

## References

  * [1] T. Wadayama and A. Nakai-Kasai, “Ordinary Differential Equation-based Sparse Signal Recovery,” International Symposium on Information Theory and Its Applications (ISITA), 2022.
  * [2] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud, “Neural ordinary differential equations,” 32nd Conference on Neural Information Processing Systems (NeurIPS), 2018.
  * [3] C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, D. Skinner, A. Ramadhan, and A. Edelman, “Universal differential equations for scientific machine learning,” arXiv:2001.04385, 2020.
  * [4] E. G. Larsson, O. Edfors, F. Tufvesson, and T. L. Marzetta, “Massive MIMO for next generation wireless systems,” IEEE Comm. Magazine, vol. 52, no. 2, pp. 186-195, 2014.
  * [5] S. Yang and L. Hanzo, “Fifty years of MIMO detection: the road to large-scale MIMOs,” IEEE Comm. Surveys and Tutorials, vol. 17, no. 4, pp. 1941-1988, 2015.
  * [6] B. Ulmann, “Analog and hybrid computer programming,” Walter de Gruyter, 2020.
  * [7] W. Haensch, T. Gokmen, and R. Puri, “The Next Generation of Deep Learning Hardware: Analog Computing,” Proc. IEEE, vol. 107, no. 1, pp. 108–122, Jan. 2019.
  * [8] J. Welser, J. W. Pitera, and C. Goldberg, “Future Computing Hardware for AI,” 2018 IEEE International Electron Devices Meeting (IEDM), pp. 1.3.1–1.3.6, Dec. 2018.
  * [9] Y. Huang, N. Guo, M. Seok, Y. Tsividis, and S. Sethumadhavan, “Analog Computing in a Modern Context: A Linear Algebra Accelerator Case Study,” in IEEE Micro, vol. 37, no. 3, pp. 30-38, 2017.
  * [10] N. Guo, Y. Huang, T. Mai, S. Patil, C. Cao, M. Seok, S. Sethumadhavan, and Y. Tsividis, “Continuous-time hybrid computation with programmable nonlinearities,” ESSCIRC Conference 2015 - 41st European Solid-State Circuits Conference (ESSCIRC), pp. 279-282, 2015
  * [11] L. Lu, J. Wu, T. Wang, and Y. SU, “Compact all-optical differential-equation solver based on silicon microring resonator,” Frontiers of Optoelectronics, vol. 5, pp. 99–106, 2012.
  * [12] H. Zhang, M. Gu, X. D. Jiang, et al. “An optical neural chip for implementing complex-valued neural network,” Nature Commun. 12, 457, 2021.
  * [13] M. Prabhu, C. Roques-Carmes, Y. Shen, et al. “Accelerating recurrent Ising machines in photonic integrated circuits,” Optica 7, pp.551-558, 2020.
  * [14] J. Capmany and D. Peŕez , “Programmable Integrated Photonics,” Oxford, 2020.
  * [15] A. Nakai-Kasai and T. Wadayama, “MMSE Signal Detection for MIMO Systems based on Ordinary Differential Equation,” IEEE Global Communications Conference (GLOBECOM), 2022
  * [16] D. F. Griffiths, and D. J. Higham, “Numerical methods for ordinary differential equations, ” Springer, 2010.
  * [17] B. Aazhang, B. P. Paris and G. C. Orsak, “Neural networks for multiuser detection in code-division multiple-access communications,” IEEE Trans. Comm., vol. 40, no. 7, pp. 1212-1222, Jul. 1992.
  * [18] E. Nachmani, Y. Beéry and D. Burshtein, “Learning to decode linear codes using deep learning,” 2016 54th Annual Allerton Conf. Comm., Control, and Computing, 2016, pp. 341-346.
  * [19] T. O’Shea and J. Hoydis, “An introduction to deep learning for the physical layer,” IEEE Trans. Cog. Comm. Net., vol. 3, no. 4, pp. 563-575, Dec. 2017.
  * [20] Y. A. LeCun, L. Bottou, G. B. Orr, and K. R. Müller, “Efficient backprop,” in Neural networks: Tricks of the trade, G. B. Orr and K. R. Müller, Eds. Springer-Verlag, London, UK, 1998, pp. 9-50.
  * [21] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, “Learning representations by back-propagating errors,” Nature, vol. 323, no. 6088, pp. 533-536, Oct. 1986.
  * [22] K. Gregor and Y. LeCun, “Learning fast approximations of sparse coding,” Proc. 27th Int. Conf. Machine Learning, pp. 399-406, 2010.
  * [23] M. Borgerding and P. Schniter, “Onsager-corrected deep learning for sparse linear inverse problems,” 2016 IEEE Global Conf. Signal and Inf. Process. (GlobalSIP), Washington, DC, Dec. 2016, pp. 227-231.
  * [24] M. Borgerding, P. Schniter, and S. Rangan, “AMP-inspired deep networks for sparse linear inverse problems, ” IEEE Trans, Sig. Process. vol 65, no. 16, pp. 4293-4308 Aug. 2017.
  * [25] A. Balatsoukas-Stimming and C. Studer, “Deep Unfolding for Communications Systems: A Survey and Some New Directions,” 2019 IEEE International Workshop on Signal Processing Systems (SiPS), pp. 266-271, 2019.
  * [26] D. Ito, S. Takabe, and T. Wadayama, ”Trainable ISTA for sparse signal recovery,” IEEE Transactions on Signal Processing, vol. 67, no. 12, pp. 3113-3125, 2019.
  * [27] D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, vol. 52, no. 4, pp. 1289-1306, Apr. 2006.
  * [28] E. J. Candes and T. Tao, “Near-optimal signal recovery from random projections: Universal encoding strategies?” IEEE Trans. Inf. Theory, vol. 52, no. 12, pp. 5406-5425, Dec. 2006.
  * [29] Z. Zhang, Y. Xu, J. Yang, X. Li, and D. Zhang, “A survey of sparse representation: Algorithms and applications,” IEEE Access, vol. 3, pp. 490-530, May. 2015.
  * [30] R. Tibshirani, “Regression shrinkage and selection via the lasso,” J. Royal Stat. Society, Series B, vol. 58, pp. 267-288, 1996.
  * [31] G. Davis, S. Mallat, and M. Avellaneda, “Adaptive greedy approximations,” Constructive Approximation, vol. 13, no. 1, pp. 57-98, Mar. 1997.
  * [32] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani, “Least angle regression,” Ann. Stat., vol. 32, no. 2, pp. 407-499, Apr. 2004.
  * [33] T. T. Wu and K. Lange, “Coordinate descent algorithms for lasso penalized regression,” Ann. Appl. Stat., vol. 2, no. 1, pp. 224-244, 2008.
  * [34] A. Chambolle, R. A. DeVore, N. Lee, and B. J. Lucier, “Nonlinear wavelet image processing: Variational problems, compression, and noise removal through wavelet shrinkage,” IEEE Trans. Image Process., vol. 7, no. 3, pp. 319-335, Mar, 1998.
  * [35] I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Comm. Pure and Appl. Math., vol. 57, no. 11, pp. 1413-1457, Aug. 2004.
  * [36] N. Parikh and S. Boyd, “Proximal algorithms,” Foundations and Trends in Optimization, vol. 1, no. 3, pp. 123-231, 2014.
  * [37] Y. Kabashima, “A CDMA multiuser detection algorithm on the basis of belief propagation,” J. Phys. A: Math. Gen., vol. 36 pp. 11111-11121, Oct. 2003.
  * [38] D. L. Donoho, A. Maleki, and A. Montanari, “Message-passing algorithms for compressed sensing,” Proceedings of the National Academy of Sciences, vol. 106, no. 45, pp. 18914-18919, Nov. 2009.
  * [39] J. Ma and L. Ping, “Orthogonal AMP,” IEEE Access, vol. 5, pp. 2020-2033, Jan. 2017.
  * [40] S. Rangan, P. Schniter, and A. K. Fletcher, “Vector approximate message passing,” 2017 IEEE Int. Symp. Inf. Theory, Jun. 2017, pp. 1588-1592.
  * [41] S. H. Strogatz, “Nonlinear dynamics and chaos: with applications to physics, biology chemistry and engineering,” Addison Wesley, 1994.
  * [42] M. Innes, “Flux: Elegant machine learning with Julia,” Journal of Open Source Software, 2018.
  * [43] J. Bezanson, S. Karpinski, B. Viral, and A. Edelman, “Julia: A fast dynamic language for technical computing,” arXiv preprint arXiv:1209.5145, 2012.
  * [44] C. Rackauckas and Q. Nie, “DifferentialEquations.jl: a performant and feature-rich ecosystem for solving differential equations in Julia, ” Journal of Open Research Software, vol. 5, no. 1, 2017.
  * [45] M. D. Buhmann, “Radial basis functions: theory and Implementations, ” Cambridge University Press, 2003.
  * [46] J. Yong and X. Y. Zhou, “Stochastic controls: Hamiltonian systems and HJB equations,” Springer, 1999.
  * [47] M. Keyanpour, and M. Azizsefat, “Numerical solution of optimal control problems by an iterative scheme,” _Advanced Modeling and Optimization_ , vol. 13, no. 1, pp. 25 – 37, 2011.
  * [48] L. Xiao and T. Zhang, “A proximal-gradient homotopy method for the sparse least-squares problem,” SIAM Journal on Optimization, vol.23, number 2, pp.1062–1091, 2013. year = 2013,