# Taylor-Lagrange Neural Ordinary Differential Equations:
Toward Fast Training and Evaluation of Neural ODEs

Franck Djeumou1 These authors contributed equally.
The code is available at https://github.com/wuwushrek/TayLaNets. Cyrus Neary
11footnotemark: 11 Eric Goubault2 Sylvie Putot2 Ufuk Topcu1
1The University of Texas at Austin, United States
2LIX, CNRS, École Polytechnique, Institut Polytechnique de Paris, France
{fdjeumou, cneary<EMAIL_ADDRESS>{goubault<EMAIL_ADDRESS>
###### Abstract

Neural ordinary differential equations (NODEs) – parametrizations of
differential equations using neural networks – have shown tremendous promise
in learning models of unknown continuous-time dynamical systems from data.
However, every forward evaluation of a NODE requires numerical integration of
the neural network used to capture the system dynamics, making their training
prohibitively expensive. Existing works rely on off-the-shelf adaptive step-
size numerical integration schemes, which often require an excessive number of
evaluations of the underlying dynamics network to obtain sufficient accuracy
for training. By contrast, we accelerate the evaluation and the training of
NODEs by proposing a data-driven approach to their numerical integration. The
proposed Taylor-Lagrange NODEs (TL-NODEs) use a fixed-order Taylor expansion
for numerical integration, while also learning to estimate the expansion’s
approximation error. As a result, the proposed approach achieves the same
accuracy as adaptive step-size schemes while employing only low-order Taylor
expansions, thus greatly reducing the computational cost necessary to
integrate the NODE. A suite of numerical experiments, including modeling
dynamical systems, image classification, and density estimation, demonstrate
that TL-NODEs can be trained more than an order of magnitude faster than
state-of-the-art approaches, without any loss in performance.

##  1 Introduction

Neural ordinary differential equations (NODEs) have recently shown tremendous
promise as a means to learn unknown continuous-time dynamical systems from
trajectory data Chen et al. (2018). By parametrizing differential equations as
neural networks, as opposed to directly fitting the available trajectory data,
NODEs provide compact representations of continuous-time systems that are
memory-efficient and that are well understood; they allow the user to harness
an existing wealth of knowledge from applied mathematics, physics, and
engineering. For example, recent works have used NODEs as a means to
incorporate physics-based knowledge into the learning of dynamical systems
Djeumou et al. (2022); Menda et al. (2019); Gupta et al. (2020); Cranmer et
al. (2020); Greydanus et al. (2019); Finzi et al. (2020); Zhong et al. (2021).
Furthermore, NODEs have been used to define continuous normalizing flows – a
class of invertible density models – to learn complex probability
distributions over data Chen et al. (2018); Grathwohl et al. (2019); Mathieu
and Nickel (2020); Salman et al. (2018).

AccuracyTrain TimeEval. TimeNFE$0$$1$$2$Relative ScoreVanilla NODETL-NODE
(Ours)TayNODET-NODE Figure 1: The proposed TL-NODE (magenta) achieves
evaluation and training times that are more than an order of magnitude faster
than state-of-the-art methods, without compromising any accuracy. The plot
illustrates the results of using TL-NODE for a classification task on the
MNIST dataset Deng (2012). All scores are relative to those obtained by a
vanilla NODE (grey), which uses an adaptive timestep numerical integrator. We
compare against TayNODE Kelly et al. (2020) (green), and T-NODE (red) which
uses a Taylor expansion for integration without the proposed correction
employed by TL-NODE. The number of function evaluations (NFE) measures the
regularity of the learned NODE (lower is better). System State Dynamics
Network$f_{\theta}(\cdot)$ Automatic Differentiation Truncated $p^{th}$ Order
Taylor Expansion of $x_{t}$ $x_{t_{i}}+\sum\nolimits_{l=1}^{p-1}\Delta
t^{l}f_{\theta}^{[l]}(x_{t_{i}})$$+$ TL-NODE Output Midpoint Prediction
Network $\Gamma_{\phi}(\cdot)$ Lagrange Form of Remainder $\Delta
t^{l}f_{\theta}^{[p]}(\Gamma)$Remainder Estimation $x_{t_{i}}$
$f_{\theta}(x_{t_{i}})$
$f^{[1]}_{\theta}(x_{t_{i}}),\ldots,f^{[p-1]}_{\theta}(x_{t_{i}})$ $\Delta t$
$\Gamma$ $f^{[p]}_{\theta}(x_{t_{i}})$ Correction for approximation error of
Taylor expansion $\hat{x}_{t_{i+1}}$ Figure 2:  Taylor-Lagrange NODE (TL-
NODE): an illustration of forward model evaluations. Given the system state
$x_{t_{i}}$ at time $t_{i}$, TL-NODE outputs a prediction of the state
$\hat{x}_{t_{i+1}}$ at future time $t_{i+1}=t_{i}+\Delta t$. The dynamics
network (yellow) parametrizes the differential equation being modeled. We use
a truncated Taylor expansion (brown) of the state dynamics $x_{t}$ to predict
the future state. A separate midpoint prediction network (blue) is trained to
estimate the remainder of the expansion, which is used as a correction for the
model’s prediction.

However, the training of NODEs can become prohibitively expensive Grathwohl et
al. (2019); Kelly et al. (2020); Finlay et al. (2020). In particular, every
forward evaluation of the NODE requires the numerical integration of the
underlying neural network parametrizing the system dynamics. Existing methods
for the training of NODEs use off-the-shelf adaptive step-size numerical
integration schemes for this purpose. However, in order to obtain sufficient
accuracy, such integration schemes have been shown in practice to require an
excessive number of evaluations of the underlying neural network. Furthermore,
the severity of this problem has been shown to grow as the training of the
neural ODE progresses; while the neural ODE learns to fit the available data,
it does not learn a representation of the dynamics that is easy to integrate
Finlay et al. (2020). These computational issues render the training of neural
ODEs on large datasets intractable, and they also prevent neural ODEs from
being deployed in applications requiring repeated fast online predictions;
such as for model-predictive control of robotic systems.

To address the above issues, we present Taylor-Lagrange NODEs (TL-NODEs); we
use a truncated Taylor expansion of the underlying neural network to predict
the future system state, and we train a separate neural network to correct
this prediction according to the Lagrange form of the expansion’s remainder.
By training the second corrector network, the approach significantly reduces
the computational cost necessary for accurate numerical integration, while
ensuring little-to-no-loss in the accuracy of the model. Figure 2 illustrates
the major components of TL-NODE, which are discussed below.

##### (1) Taylor expansions for the numerical integration of the NODE.

In order to integrate the NODE, we use a fixed-order Taylor expansion of the
dynamical system in time. We take advantage of Taylor-mode automatic
differentiation to efficiently compute the higher-order terms of the expansion
on a GPU Bettencourt et al. (2019). By specifying the number of terms to
include in the expansion, we ensure that only a fixed number of evaluations of
the underlying dynamics network are required per training step.

##### (2) Correcting the expansion’s approximation error.

We use the Lagrange form of the expansion’s remainder to define a correction
term, which greatly improves the accuracy of the predictions obtained from the
truncated expansion. To estimate this approximation error efficiently, we
propose to use a neural network to predict the so-called midpoint value – a
point near the center of the expansion at which the approximation error can be
evaluated exactly. While learning this midpoint value may, in general, be as
complex as learning the neural ODE itself, we derive explicit formulas for the
midpoint using assumptions on the regularity of the dynamics. These
expressions reduce the complexity of the learning problem; only one unknown
term in the expression must be learned. We provide upper bounds on the error
of the proposed Taylor-Lagrange expansion, in terms of the error in the
predicted midpoint value and the order of the expansion.

We demonstrate the effectiveness of the proposed approach through a suite of
numerical experiments. The experimental tasks include the integration of known
dynamics, learning to predict unknown dynamics, supervised classification, and
density estimation. Figure 1 illustrates the result of applying TL-NODE to a
classification task. Across all experiments we observe that the training and
the evaluation of TL-NODEs is more than an order of magnitude faster than
existing NODE methods, while incurring no loss in performance.

##### Related Work.

Similarly to our work, a number of papers study how to reduce the
computational cost of training neural ODEs. After Chen et al. (2018) initially
presented the neural ODE, Grathwohl et al. (2019) proposed a stochastic
estimator of the likelihood to reduce computational cost when using neural
ODEs for continuous normalizing flows. Kelly et al. (2020); Finlay et al.
(2020); Pal et al. (2021) propose additional regularization terms to learn
neural ODEs that are easy to integrate. Ghosh et al. (2020) propose to
regularize the neural ODE by randomly sampling the end time of the ODE during
training. However, all of these works use off-the-shelf numerical integration
algorithms for the forward evaluation of the NODE. By contrast, our work
suggests a novel data-driven integration scheme, resulting in training times
that are an order of magnitude faster than the current state-of-the-art.

Meanwhile, Poli et al. (2020) also suggest training an additional corrector
neural network to speed up the numerical integration of NODEs. However, they
do not present a technique that is able to apply such a corrector network
during the training of the NODE. By contrast, we propose algorithms for the
simultaneous training of the dynamics network and the remainder estimation
network; this simultaneous training results not only in a speedup of the NODE
evaluations, but also in a speedup of the NODE’s training. Furthermore, we
propose a method to simplify the learning of the remainder term by taking
advantage of regularity assumptions on the system dynamics. This
simplification leads to more efficient and generalizable learning of the
correction term.

##  2 Background

We begin by introducing necessary background on ordinary differential
equations (ODEs) and on neural ODEs.

##### Ordinary Differential Equations (ODEs).

Let the function
$f(x,t):\mathbb{R}^{n}\times\mathbb{R}_{+}\mapsto\mathbb{R}^{n}$ be
_Lipschitz-continuous_ in $x$ and $t$. An ordinary differential equation
specifies the instantaneous change of a vector-valued signal of time
$x:\mathbb{R_{+}}\mapsto\mathbb{R}^{n}$.

$\displaystyle\dot{x}(t)=f(x(t),t)$ (1)

We note that in general, the explicit dependence of $f(x,t)$ on $t$ can be
removed by adding a dimension to the state variable $x$. As such, throughout
the remainder of the paper, we consider autonomous systems of the form
$\dot{x}(t)=f(x(t))$. Furthermore, for notational simplicity we use the
subscript notation $x_{t}$ in place of $x(t)$.

Given some initial state $x_{0}\vcentcolon=x_{t_{0}}\in\mathbb{R}^{n}$ at
initial time $t_{0}\geq 0$, we wish to compute a solution to (1). In this
work, we are specifically interested in predicting the value of the state
$x_{t}$, at an arbitrary future point in time $T\geq t_{0}$. The value of
$x_{T}$ can be found by integrating or solving the ODE:
$x_{T}=x_{0}+\int_{t_{0}}^{T}f(x_{s})\,ds.$

##### Neural Ordinary Differential Equations (NODEs).

Neural ODEs (NODEs) are a class of deep learning models that use a neural
network $f_{\theta}(x)$ to parametrize an ODE, and that subsequently use
numerical integration algorithms to evaluate the model’s outputs. More
specifically, given a NODE input tuple $(x_{0},t_{0},T)$ consisting of an
initial state $x_{0}$, an initial time $t_{0}$, and a prediction time $T$, the
output $\hat{x}_{T}\vcentcolon=\mathrm{NODE}_{\theta}(x_{0},t_{0},T)$ of the
NODE with parameter vector $\theta$ is given by
$\mathrm{ODESolve}(f_{\theta},x_{0},t_{0},T)\approx
x_{0}+\int_{t_{0}}^{T}f_{\theta}(x_{s})ds$. Here,
$\mathrm{ODESolve}(f_{\theta},x_{0},t_{0},T)$ is a numerical approximation of
the solution to the dynamics parametrized by $f_{\theta}(x)$.

We note that the particular algorithm used in place of
$\mathrm{ODESolve}(\cdot)$ influences the model’s accuracy, the computational
cost of forward evaluations of the model, and the computational cost of
training the model. Existing implementations of neural ODEs have typically
relied on existing adaptive-step numerical integration algorithms for this
purpose. By contrast, this work proposes the use of a novel Taylor expansion
based method that uses data-drive estimations of the expansion’s truncation
error to evaluate and train neural ODEs efficiently and accurately.

##  3 Taylor-Lagrange Neural Ordinary Differential Equations (TL-NODE)

In this section we propose TL-NODEs for the efficient and accurate training
and evaluation of neural ODEs.

###  3.1 Data-Driven Taylor-Lagrange Numerical Integration of the NODE
Dynamics

Our objective is to efficiently and accurately evaluate neural ODEs through
numerical integration of $f_{\theta}(\cdot)$. Toward this end, we propose to
make direct use of the Taylor-Lagrange expansion of $x_{t}$ in time. Taylor
expansions are techniques used to approximate a function near a particular
expansion point using polynomials. We use the term Taylor-Lagrange expansion
to refer to the truncated Taylor expansion including the Lagrange form of the
remainder. In order to obtain highly accurate expansions of $x_{t}$ without
requiring an excessive number of evaluations of $f_{\theta}(\cdot)$, we
propose to train a separate neural network with parameters $\phi$ to estimate
the remainder term in the expansion. We give a step-by-step description of the
proposed methodology below.

##### Partitioning the prediction interval.

We begin by partitioning the prediction interval $[t_{0},T]$ into a collection
of $H$ sub-intervals $[t_{i},t_{i+1}]$. For notational convenience, we assume
the value of the integration time step is fixed, i.e. $\Delta t=t_{i+1}-t_{i}$
for all $i=0,\ldots,H-1$. Given the value of the state $x_{t_{i}}$ at time
$t_{i}$, we compute an approximation of the value of the state at the next
time step $\hat{x}_{t_{i+1}}$ via a $p^{th}$ order Taylor-Lagrange expansion
of $x_{t}$ about time $t_{i}$, which we denote by $TL^{\Delta
t}_{\theta,\phi}(x_{t_{i}})$. The output
$\mathrm{TLNODE}_{\theta,\phi}(x_{0},t_{0},T)$ of the neural ODE is computed
by iteratively using $TL^{\Delta t}_{\theta,\phi}(\cdot)$ to approximate the
integral of $f_{\theta}(\cdot)$ over each sub-interval.

##### Expressing $\boldsymbol{TL^{\Delta t}_{\theta,\phi}(\cdot)}$ in Terms of
$\boldsymbol{f_{\theta}(\cdot)}$.

Note that because $f_{\theta}(x)$ estimates the time derivative of the system
state $x$, the Taylor-Lagrange expansion of $x_{t}$ may be expressed in terms
of the Taylor coefficients $f_{\theta}^{[l]}(x)$, which are recursively
defined through the equations $f^{[1]}_{\theta}(x)=f_{\theta}(x)$ and
$f^{[l+1]}_{\theta}(x)=\frac{1}{l+1}[\frac{\partial f^{[l]}_{\theta}}{\partial
x}f_{\theta}](x)$. Equation (2) accordingly presents the Taylor-Lagrange
expansion of $x_{t}$ about the point in time $t_{i}$, evaluated at the future
time $t_{i+1}=t_{i}+\Delta t$.

$\begin{split}TL_{\theta,\phi}^{\Delta
t}(x_{t_{i}})\vcentcolon=x_{t_{i}}+\sum\nolimits_{l=1}^{p-1}&\Delta
t^{l}f_{\theta}^{[l]}(x_{t_{i}})\\\
&+\mathcal{R}_{\phi}(f_{\theta},x_{t_{i}},\Delta t)\end{split}$ (2)

The first two terms on the right hand side of (2) make up the truncated Taylor
expansion of $x_{t}$, while $\mathcal{R}_{\phi}(f_{\theta},x_{t_{i}},\Delta
t)$ denotes an estimation of the remainder of this truncated expansion. More
specifically, $\mathcal{R}_{\phi}(f_{\theta},x_{t_{i}},\Delta t)$ estimates
the approximation error of the $p^{th}$ order expansion; if we could known
this value exactly, (2) would provide an exact evaluation of the integral
$x_{t_{i}}+\int_{t_{i}}^{t_{i+1}}f_{\theta}(x_{s})ds$. Below we propose a
methodology to learn to accurately estimate the value of this remainder term,
given $f_{\theta}$, $x_{t_{i}}$, and $\Delta t$ as inputs.

##### Estimating the Remainder Term $\boldsymbol{\mathcal{R}_{\phi}(\cdot)}$.

To obtain accurate and generalizable estimations of the remainder term
$\mathcal{R}_{\phi}(f_{\theta},x_{t_{i}},\Delta t)$, we begin by using
Taylor’s theorem to express it as
$\mathcal{R}_{\phi}(f_{\theta},x_{t_{i}},\Delta t)=f_{\theta}^{[p]}(\Gamma)$.
Here, $\Gamma\in\mathbb{R}^{n}$ denotes the midpoint of the Taylor-Lagrange
expansion. More specifically, there exists some point in time $\xi$ with
$t_{i}\leq\xi\leq t_{i+1}$ such that when we define
$\Gamma\vcentcolon=x_{\xi}$, then $\Delta t^{p}f_{\theta}^{[p]}(\Gamma)$
provides the exact value of the approximation error of the expansion.

While no closed form expression for the midpoint $\Gamma$ exists, we propose
to learn to predict its value given the state $x_{t_{i}}$ and the time step
$\Delta t$. Learning to predict $\Gamma$ directly as a function of these
inputs is a challenging problem in general. We instead propose to use the
result of Theorem 1, which provides a closed-form expression for $\Gamma$ in
terms of some unknown term $\bar{\Gamma}\in\mathbb{R}^{n}$. By taking
advantage of this expression for $\Gamma$, we greatly simplify the task of
learning to predict its value.

######  Theorem 1 (Simplified Midpoint Expression).

If $f_{\theta}$ is a Lipschitz-continuous function, then there exists a
function
$\bar{\Gamma}:\mathbb{R}^{n}\times\mathbb{R}_{+}\to\mathbb{R}^{n\times n}$
such that the midpoint value $\Gamma$ of the Taylor-Lagrange expansion
$TL^{\Delta t}_{\theta,\phi}(x_{t_{i}})$ is related to $x_{t_{i}}$,
$f_{\theta}(x_{t_{i}})$, and $\bar{\Gamma}(x_{t_{i}},\Delta t)$ through

$\Gamma=x_{t_{i}}+\bar{\Gamma}(x_{t_{i}},\Delta t)\odot
f_{\theta}(x_{t_{i}}),$ (3)

where $\odot$ denotes matrix-vector multiplication.

Theorem 1 is obtained using tools from interval Taylor-Langrange based
reachability analysis Djeumou et al. (2021). More specifically, the equation
given in (3) is derived from explicit formulas for the so-called _apriori_
enclosure – a set derived from the local Lipschitzness of $f_{\theta}$ that is
guaranteed to contain the value of $TL^{\Delta t}_{\theta,\phi}(x_{t_{i}})$. A
proof of Theorem 1 is provided in Appendix A. We also prove that for linear
dynamics, $\bar{\Gamma}$ does not depend on $x_{t_{i}}$.

##### Estimating the Midpoint Value.

Given the result of Theorem 1, we propose to parameterize the unknown function
$\bar{\Gamma}_{\phi}(\cdot)$ using a neural network with parameters $\phi$.
For notational simplicity, we use $\Gamma_{\phi}(x,\Delta t)$ to denote the
value of the right hand side of (3) when $\bar{\Gamma}(\cdot)$ is approximated
by $\bar{\Gamma}_{\phi}(\cdot)$. Given the predicted midpoint value
$\Gamma_{\phi}(x,\Delta t)$, we estimate the remainder term of the $p^{th}$
order Taylor-Lagrange expansion $TL^{\Delta t}_{\theta,\phi}(x_{t_{i}})$ as
$\mathcal{R}_{\phi}(f_{\theta},x_{t_{i}},\Delta t)\approx\Delta
t^{p}f_{\theta}^{[p]}(\Gamma_{\phi}(x_{t_{i}},\Delta t))$.

##### The Proposed TL-NODE Evaluation Algorithm.

Algorithm 1 summarizes the proposed approach for the numerical evaluation of
neural ODEs. In lines $1$ and $2$, the prediction interval $[t_{0},T]$ is
broken into $H$ sub-intervals. The for loop in lines $3-6$ iterates over these
sub-intervals, and uses the midpoint prediction network $\Gamma_{\phi}(\cdot)$
to estimate the midpoint value (line $4$), before using this estimate to
approximate the state value $\hat{x}_{t_{i+1}}$ at the end of the sub-interval
(line 5).

Algorithm 1 Evaluating $\mathrm{TLNODE}_{\theta,\phi}(x_{0},t_{0},T)$

Input: $x_{0}$, $t_{0}$, $T$
Parameter: $\theta$, $\phi$, $p$, $H$
Output: Model prediction $\hat{x}_{T}$.

1: $\hat{x}_{t_{0}}\leftarrow x_{0}$; $\Delta t\leftarrow\frac{T-t_{0}}{H}$

2: for $i=0,1,\ldots,H$ do{ $t_{i}\leftarrow t_{0}+i\Delta t$}

3: for $i=0,1,\ldots,H-1$ do

4: $\Gamma\leftarrow\hat{x}_{t_{i}}+\bar{\Gamma}_{\phi}(\hat{x}_{t_{i}},\Delta
t)\odot f_{\theta}(\hat{x}_{i})$

5: $\hat{x}_{t_{i+1}}\leftarrow\hat{x}_{t_{i}}+\sum_{l=1}^{p-1}\Delta
t^{l}f_{\theta}^{[l]}(\hat{x}_{i})+\Delta t^{p}f_{\theta}^{[p]}(\Gamma)$

6: end for

7: return $\hat{x}_{H}$

##### Bounding the Error of the TL-NODE Evaluation Algorithm.

Given a fixed dynamics function $f_{\theta}(\cdot)$, we seek to bound the
error on a $p^{th}$ order Taylor-Lagrange expansion which uses a learned
midpoint value predictor $\Gamma_{\phi}(\cdot)$ to estimate the expansion’s
remainder $\mathcal{R}_{\phi}(\cdot)$. Such an error bounds straightforwardly
depends on how well $\Gamma_{\phi}(\cdot)$ approximates the true midpoint
$\Gamma$, as described in Theorem 2. A proof of Theorem 2 is provided in
Appendix C.

######  Theorem 2 (Integration Accuracy).

If the learned midpoint predictor $\Gamma_{\phi}(\cdot)$ is a
$\mathcal{O}(\eta)$ approximator to the midpoint $\Gamma$ of the Taylor-
Lagrange expansion of $TL^{\Delta t}_{\theta,\phi}(x_{t_{i}})$, then
$\|x_{t_{i+1}}-TL^{\Delta t}_{\theta,\phi}(x_{t_{i}})\|\leq c\eta\Delta t^{p}$
for some $c>0$ that depends on $f_{\theta}$.

##### A Note on the Evaluating the Taylor Coefficients.

The Taylor coefficients
$f_{\theta}^{[1]}(\cdot),\ldots,f_{\theta}^{[p]}(\cdot)$ can in principle be
evaluated using repeated application of forward-mode automatic differentiation
to iteratively compute the Jacobian-vector products $[\frac{\partial
f^{[l]}_{\theta}}{\partial x}f_{\theta}](x)$. However, doing so would incur a
time cost of $\mathcal{O}(exp(p))$. We instead use Taylor mode automatic
differentiation, which computes the first $p$ Taylor coefficients
$f^{[1]}_{\theta}(\cdot),\ldots,f^{[p]}_{\theta}(\cdot)$ in a single pass,
with a time cost of only $\mathcal{O}(p^{2})$ or of $\mathcal{O}(p\log p)$,
depending on the underlying operations involved Griewank and Walther (2008);
Bettencourt et al. (2019); Kelly et al. (2020).

###  3.2 Training Taylor-Lagrange Neural Ordinary Differential Equations

Given a training dataset $\mathcal{D}$, we wish to train both components of
the TL-NODE: the dynamics network $f_{\theta}(\cdot)$ and the midpoint
prediction network $\Gamma_{\phi}(\cdot)$. To do so, we propose an algorithm
that alternates between training each of the components via stochastic
gradient descent while keeping the parameters of the other component fixed. We
assume that each datapoint within the dataset
$\mathcal{D}=\\{(x_{0}^{j},t_{0}^{j},T^{j},y^{j})\\}_{j=1}^{|\mathcal{D}|}$ is
comprised of an initial state $x^{j}$, an initial time $t_{0}^{j}$, a
prediction time $T^{j}$, and a labeled output value $y^{j}$.

##### Training the Dynamics Network $\boldsymbol{f_{\theta}(\cdot)}$.

We begin by holding the parameter vector $\phi$ of the midpoint prediction
network to some fixed value $\hat{\phi}$, and training the dynamics network
$f_{\theta}(\cdot)$ by solving the optimization problem (4) via stochastic
gradient descent.

$\begin{split}\min_{\theta}\sum_{(x_{0}^{j},t_{0}^{j},T^{j},y^{j})\in\mathcal{D}}&\Bigl{[}\mathcal{L}(\mathrm{TLNODE}_{\theta,\hat{\phi}}(x_{0}^{j},t_{0}^{j},T^{j}),y^{j})\Bigr{.}\\\
&\Bigl{.}+\lambda\sum_{i=0}^{H-1}||\Delta
t^{p}f_{\theta}^{[p]}(\Gamma_{\hat{\phi}}(\hat{x}_{t_{i}},\Delta
t))||^{2}\Bigr{]}\end{split}$ (4)

Here, $\mathcal{L}(\cdot)$ is any differentiable loss function and
$\hat{x}_{t_{i}}$ denotes integrator-estimated intermediate state at time
$t_{i}$.

The summation in the second line of (4) measures the magnitude of the
remainder terms $\mathcal{R}_{\phi}(\cdot)$ of the truncated Taylor-Lagrange
expansions used for numerical integration. We may interpret this penalty term
as having two purposes. Firstly, it acts as a regularizer that penalizes the
higher-order derivatives of $f_{\theta}(\cdot)$ during training. Intuitively,
by penalizing these higher order derivatives we encourage solutions that fit
the data while also remaining as simple as possible. Secondly, the penalty
term prevents the TL-NODE from using $\mathcal{R}_{\phi}(\cdot)$ to overfit
the training data. By ensuring that the remainder term of the Taylor-Lagrange
expansion remains small during training, we learn a dynamics function
$f_{\theta}(\cdot)$ whose truncated expansions fit the training data as well
as possible, while using $\mathcal{R}_{\phi}(\cdot)$ only for small
corrections.

##### Training the Midpoint Prediction Network
$\boldsymbol{\Gamma_{\phi}(\cdot)}$.

Recall that the midpoint prediction network $\Gamma_{\phi}(\cdot)$ plays a
crucial role in accurately integrating the dynamics specified by
$f_{\theta}(\cdot)$. So, as the parameters of the dynamics network
$f_{\theta}(\cdot)$ are updated throughout training, our estimates of
$\Gamma$, the midpoint of the Taylor-Lagrange expansion of
$f_{\theta}(\cdot)$, should be updated accordingly. We thus propose to
occasionally freeze the parameters of the dynamics network $\hat{\theta}$ in
order to train $\Gamma_{\phi}(\cdot)$.

After fixing $\hat{\theta}$, we begin by generating a small dataset
$\mathcal{D}_{\hat{\theta}}$. The datapoints of $\mathcal{D}_{\hat{\theta}}$
correspond to solutions of the ODE encoded by the fixed dynamics network
$f_{\hat{\theta}}(\cdot)$. That is, for each
$(x_{0},t_{0},T,y)\in\mathcal{D}_{\hat{\theta}}$ the output label $y$ is given
by $\mathrm{ODESolve}(f_{\hat{\theta}},x_{0},t_{0},T)$, where
$\mathrm{ODESolve}(\cdot)$ is a highly accurate adaptive-step ODE solver. Once
the dataset $\mathcal{D}_{\hat{\theta}}$ has been generated, we train
$\Gamma_{\phi}(\cdot)$ by using stochastic gradient descent to solve the
optimization problem (5).

$\begin{split}\min_{\phi}\sum_{(x_{0}^{j},t_{0}^{j},T^{j},y^{j})\in\mathcal{D}_{\hat{\theta}}}&||\mathrm{TLNODE}_{\hat{\theta},\phi}(x_{0}^{j},t_{0}^{j},T^{j})-y^{j}||^{2}\end{split}$
(5)

##### The Proposed TL-NODE Training Algorithm.

Algorithm 2 details the proposed training procedure. Throughout training we
alternate between the following two subroutines: (lines $3$-$4$) fix $\phi$
and take $N_{\theta}$ stochastic gradient descent steps to train the dynamics
network $f_{\theta}(\cdot)$ according to (4), (lines $5$-$7$) fix $\theta$ and
take $N_{\phi}$ stochastic gradient descent steps to train the midpoint
prediction network $\Gamma_{\phi}(\cdot)$ according to (5).

Algorithm 2 Training the TL-NODE

Input: Training dataset $\mathcal{D}$
Parameter: $N_{\theta}$, $N_{\phi}$, $N_{train}$,
$N_{|\mathcal{D}_{\hat{\theta}}|}$
Output: Model parameters $\theta$, $\phi$

1: Initialize parameters $\theta$, $\phi$

2: for $N_{train}$ steps do

3: Fix $\hat{\phi}\leftarrow\phi$

4: for $N_{\theta}$ steps do{ $\theta\leftarrow\mathrm{sgdStep}(\text{Eq.
}\eqref{eq:loss_opt_problem},\theta,\hat{\phi},\mathcal{D})$}

5: Fix $\hat{\theta}\leftarrow\theta$;
$\\{(x_{0}^{j},t_{0}^{j},T^{j})\\}_{j}\leftarrow\mathrm{Sample}(\mathcal{D},N_{|\mathcal{D}_{\hat{\theta}}|})$

6:
$\mathcal{D}_{\hat{\theta}}\leftarrow\mathrm{ODESolve}(f_{\hat{\theta}},\\{(x_{0}^{j},t_{0}^{j},T^{j})\\}_{j})$

7: for $N_{\phi}$ steps do{ $\phi\leftarrow\mathrm{sgdStep}(\text{Eq.
}\eqref{eq:midpoint_loss_opt_problem},\hat{\theta},\phi,\mathcal{D}_{\hat{\theta}})$}

8: end for

9: return $\theta,\phi$

##  4 Numerical Experiments

Method | Train Accuracy (%) | Test Accuracy (%) | Train Time (min) | Eval. Time (ms) | NFE
---|---|---|---|---|---
TL-NODE (ours) | 99.96 | 98.23 | 2.55 | 1.04 | 62
Vanilla NODE | 99.33 | 97.87 | 42.7 | 16 | 110.6
TayNODE | 99.29 | 98.02 | 94.3 | 11 | 80.00
RNODE | 98.72 | 97.74 | 10.2 | 2.05 | 98.0
SRNODE* | 100.0 | 98.08 | 98.1 | - | 259.0
STEER* | 100.0 | 97.94 | 103 | - | 265.0
Table 1: MNIST image classification results.

We demonstrate the effectiveness of TL-NODE through several numerical
experiments: the numerical integration of known dynamics, the learning of
unknown dynamics, a supervised classification task, and a density estimation
task. As an initial illustrative example we apply TL-NODE to linear dynamics.
However, we note that the latter classification and density estimation tasks
involve non-linear, time-dependent, and high-dimensional dynamics. Additional
experimental details – including hyperparameter selection – are included in
Appendix D.

0.00.050.10.150.20.250.3$10^{-4}$$10^{-3}$Eval. Time [s]HypersolverTaylor,
$p=1$RK4TL (ours), $p=1$TL (ours),
$p=2$Dopri50.00.050.10.150.20.250.3$10^{-12}$$10^{-8}$$10^{-4}$$10^{0}$Prediction
Time $T-t_{0}$Integration Error Figure 3: Numerical integration of known stiff
dynamics. All algorithms other than Dopri5 use $H=1$. Top: average time for
numerical integration, as a function of the size of the prediction time
interval. Bottom: average normalized integration error. Averages are taken
with respect to 250 randomly sampled initial states.

###  4.1 Modeling a Dynamical System

We begin by applying TL-NODE to the task of modeling a stiff dynamical system.
More specifically, we use the proposed Taylor-Lagrange approach to learn, and
to integrate, the ODE $\dot{x}=Ax$, where $x\in\mathbb{R}^{2}$ and
$A\in\mathbb{R}^{2\times 2}$ has eigenvalues $\lambda_{1}=-1$ and
$\lambda_{2}=-1000$.

####  4.1.1 Integration of Known Stiff Dynamics.

To examine the accuracy and robustness of the midpoint prediction network
$\Gamma_{\phi}(\cdot)$, we begin by assuming the dynamics function $f(x)=Ax$
is known, and we use the proposed Taylor-Lagrange numerical integration method
to predict future system states. We note that because we assume $f(\cdot)$ is
known, there is no need to parameterize the system dynamics using a neural
network $f_{\theta}(\cdot)$. However, we may still apply the method outlined
in §3.2 to train $\Gamma_{\phi}(\cdot)$ to predict the approximation error of
the Taylor-Lagrange expansions.

Baselines. We apply both $1^{st}$ and $2^{nd}$ order Taylor expansions for
numerical integration. For comparison, we include the results of a fixed-step
RK4 method, an adaptive-step method (Dopri5), and the Hypersolver method Poli
et al. (2020). The tolerance parameters $rtol$ and $atol$ of the adaptive-step
Dopri5 integrator are both set to $1.4e^{-12}$. We also plot the result of
using a Taylor expansion for integration, without including the learned
approximation error term.

Results. Figure 3 illustrates the numerical integration results. For brevity,
in the figure we refer to the proposed Taylor-Lagrange method for integration
as TL. We observe that TL-NODE enjoys lower integration error than all of the
baseline methods except for Dopri5. However, Dopri5 requires computation times
that are more than an order of magnitude higher than that of our method. We
additionally observe that while the Hypsersolver method requires similar
computation time to TL-NODE, the error of its numerical integration results
are several orders of magnitude higher. Furthermore, we note that for any
prediction time intervals $T-t_{0}$ larger than $0.05(s)$, the fixed-step RK4,
Truncated Taylor expansion method, and Hypersolver method all have normalized
prediction errors values of $1.0$ (the highest possible value). By contrast,
our TL-NODE approach achieves an average error value on the order of
$10^{-4}$, even when $T-t_{0}=0.3(s)$. This demonstrataes the robustness of
the proposed approach to the size of the prediction interval.

Finally, we note that the integration error of the truncated Taylor expansion
method (yellow) is several orders of magnitude larger than that of TLN. The
only difference between these methods is TLN’s inclusion of the proposed
correction term that learns the approximation error, demonstrating the gain in
accuracy that this learned correction term provides.

####  4.1.2 Learning Unknown Dynamics

We now assume that the system dynamics $f(x)=Ax$ are unknown and train a TL-
NODE to model the dynamical system.

Baselines. We compare to Vanilla NODE, which uses adaptive-step Dopri5 for
numerical integration, to a NODE trained using fixed-step RK4, and to T-NODE –
a version of our approach that also uses Taylor expansions for integration,
but does not estimate their remainder term.

Results. Figure 4 illustrates the NODE’s average prediction error as a
function of the number of elapsed training steps. TL-NODE achieves smiliar
prediction error values to the Vanilla NODE throughout training, while the
prediction errors of the other two baseline methods are twice as large. The
wall-clock training time for TL-NODE is 31.9s, for the Vanilla NODE it is
609.8s, for the RK4 NODE it is 35.8s, and for T-NODE it is 21.2s. Algorithm 2
effectively balances the training of TL-NODE’s two components:
$f_{\theta}(\cdot)$ and $\Gamma_{\phi}(\cdot)$. The result is a dynamics model
that is as accurate as the Vanilla NODE trained using Dopri5, but whose
training and evaluation times are much faster.

$0$$2$$4$$6$$8$$10$$12$$14$$16$$18$$20$1.05.010Training Steps $(\cdot
10^{3})$Prediction Error $(\cdot 10^{-5})$TL-NODE (ours), $p=1$RK4T-NODE,
$p=2$Vanilla NODE Figure 4:  Predicting unknown dynamics over a prediction
time step of $T-t_{0}=0.01(s)$. We plot the average mean square error of the
predicted state as a function of the elapsed training steps.

###  4.2 Supervised Classification

We train a TL-NODE model to perform classification on the MNIST dataset Deng
(2012). Our model follows the architecture of the neural ODE-based MNIST
classifier presented in the work of Kelly et al. (2020) and further used in
Pal et al. (2021) for benchmarking. Specifically, the model uses a two-layered
neural network of size $100$ and $728$ (size of the images) with sigmoid-based
non linearities to parameterize the dynamics function $f_{\theta}(\cdot)$. The
NODE outputs propagate through a linear classifier to estimate of the image
labels.

Baselines. We compare the proposed Taylor-Lagrange networks with state-of-the-
art NODE algorithms. More specifically, we compare to RNODE Finlay et al.
(2020) and TayNode Kelly et al. (2020) using the source code provided by the
respective authors. We implicitly also compare our results with other
regularization techniques such STEER Ghosh et al. (2020) and SRNODE Pal et al.
(2021) thanks to the thorough experiments provided in Pal et al. (2021) for a
similar model of the MNIST classification problem.

Results. Table 1 lists the experimental results. TL-NODE achieves evaluation
and training times that are more than an order of magnitude faster than the
baseline approaches, while also achieving the highest accuracy on the test
dataset. TL-NODE also learns a dynamics network $f_{\theta}(\cdot)$ that
requires the smallest number of function evaluations (NFE) when it is being
numerically integrated using an adaptive-step integrator. The low NFE score of
TL-NODE indicates that the regularization term in (4) is effective at
producing learned dynamics networks $f_{\theta}(\cdot)$ that are easy to
numerically integrate.

###  4.3 Density Estimation

We apply TL-NODEs to train continuous-normalizing-flow-based generative models
Chen et al. (2018); Grathwohl et al. (2019) to approximate the distribution of
the MiniBooNE dataset Roe et al. (2005); Papamakarios et al. (2017).

Method | Loss (nat) | Train Time (min) | NFE
---|---|---|---
TL-NODE (ours) | 9.62 | 12.3 | 167.9
Vanilla NODE | 9.74 | 59.7 | 183.8
TayNODE | 9.75 | 148.3 | 168.2
RNODE | 9.78 | 10.32 | 182.0
Table 2: Density estimation results.

Results. TL-NODE achieves the best loss score and the lowest required number
of function evaluations (NFE) in comparison with the baseline approaches.

##  5 Conclusions

We present Taylor-Lagrange Neural Ordinary Differential Equations (TL-NODEs):
a class of neural ODEs (NODEs) that use fixed-order Taylor expansions for
numerical integration during NODE training and evaluation. TL-NODEs also train
a separate neural network to predict the expansion’s remainder, which is used
as a correction term to improve the accuracy of the NODE’s outputs. We
demonstrate that TL-NODEs enjoy evaluation and training times that are an
order of magnitude faster than the current state-of-the-art, without any loss
in accuracy. Future work will aim to apply the accelerated NODE evaluation
times to the online model-based control of unknown dynamical systems.

## Acknowledgements

This work was supported in part by ARL W911NF2020132, AFOSR FA9550-19-1-0005,
and NSF 1646522.

## References

  * Bettencourt et al. [2019] J. Bettencourt, Matthew J. Johnson, and D. Duvenaud.  Taylor-mode automatic differentiation for higher-order derivatives in jax.  In Workshop on Program Transformations for ML, NeurIPS, 2019.
  * Chen et al. [2018] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud.  Neural ordinary differential equations.  Advances in Neural Information Processing Systems, 2018.
  * Cranmer et al. [2020] M. Cranmer, Sam Greydanus, Stephan Hoyer, Peter W. Battaglia, David N. Spergel, and Shirley Ho.  Lagrangian neural networks.  ArXiv, abs/2003.04630, 2020.
  * Deng [2012] Li Deng.  The mnist database of handwritten digit images for machine learning research.  IEEE Signal Processing Magazine, 2012.
  * Djeumou et al. [2021] Franck Djeumou, Abraham P Vinod, Eric Goubault, Sylvie Putot, and Ufuk Topcu.  On-the-fly control of unknown smooth systems from limited data.  In 2021 American Control Conference. IEEE, 2021.
  * Djeumou et al. [2022] Franck Djeumou, Cyrus Neary, Eric Goubault, Sylvie Putot, and Ufuk Topcu.  Neural networks with physics-informed architectures and constraints for dynamical systems modeling.  In Learning for Dynamics and Control. PMLR, 2022.
  * Finlay et al. [2020] Chris Finlay, J. Jacobsen, L. Nurbekyan, and Adam M. Oberman.  How to train your neural ode: the world of jacobian and kinetic regularization.  In International Conference on Machine Learning. PMLR, 2020.
  * Finzi et al. [2020] Marc Finzi, Ke Alexander Wang, and Andrew G Wilson.  Simplifying hamiltonian and lagrangian neural networks via explicit constraints.  Advances in Neural Information Processing Systems, 2020.
  * Ghosh et al. [2020] Arna Ghosh, Harkirat Singh Behl, Emilien Dupont, Philip H. S. Torr, and Vinay Namboodiri.  Steer : Simple temporal regularization for neural odes.  ArXiv, abs/2006.10711, 2020.
  * Grathwohl et al. [2019] Will Grathwohl, Ricky T. Q. Chen, Jesse Bettencourt, Ilya Sutskever, and David Kristjanson Duvenaud.  Ffjord: Free-form continuous dynamics for scalable reversible generative models.  ArXiv, abs/1810.01367, 2019.
  * Greydanus et al. [2019] Samuel Greydanus, Misko Dzamba, and Jason Yosinski.  Hamiltonian neural networks.  Advances in Neural Information Processing Systems, 2019.
  * Griewank and Walther [2008] Andreas Griewank and Andrea Walther.  Evaluating derivatives: principles and techniques of algorithmic differentiation.  SIAM, 2008.
  * Gupta et al. [2020] Jayesh K Gupta, Kunal Menda, Zachary Manchester, and Mykel Kochenderfer.  Structured mechanical models for robot learning and control.  In Learning for Dynamics and Control. PMLR, 2020.
  * Kelly et al. [2020] Jacob Kelly, Jesse Bettencourt, Matthew J Johnson, and David K Duvenaud.  Learning differential equations that are easy to solve.  Advances in Neural Information Processing Systems, 2020.
  * Kingma and Ba [2014] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * Mathieu and Nickel [2020] Emile Mathieu and Maximilian Nickel.  Riemannian continuous normalizing flows.  ArXiv, abs/2006.10605, 2020.
  * Menda et al. [2019] Kunal Menda, Jayesh K Gupta, Zachary Manchester, and Mykel J Kochenderfer.  Structured mechanical models for efficient reinforcement learning.  In Workshop on Structure and Priors in Reinforcement Learning, International Conference on Learning Representations, 2019.
  * Pal et al. [2021] Avik Pal, Yingbo Ma, Viral Shah, and Christopher V Rackauckas.  Opening the blackbox: Accelerating neural differential equations by regularizing internal solver heuristics.  In International Conference on Machine Learning. PMLR, 2021.
  * Papamakarios et al. [2017] George Papamakarios, Iain Murray, and Theo Pavlakou.  Masked autoregressive flow for density estimation.  ArXiv, abs/1705.07057, 2017.
  * Poli et al. [2020] Michael Poli, Stefano Massaroli, Atsushi Yamashita, Hajime Asama, and Jinkyoo Park.  Hypersolvers: Toward fast continuous-depth models.  ArXiv, abs/2007.09601, 2020.
  * Roe et al. [2005] Byron P Roe, Hai-Jun Yang, Ji Zhu, Yong Liu, Ion Stancu, and Gordon McGregor.  Boosted decision trees as an alternative to artificial neural networks for particle identification.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 2005.
  * Salman et al. [2018] Hadi Salman, Payman Yadollahpour, Tom Fletcher, and Nematollah Batmanghelich.  Deep diffeomorphic normalizing flows.  ArXiv, abs/1810.03256, 2018.
  * Zhong et al. [2021] Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty.  Benchmarking energy-conserving neural networks for learning dynamics from data.  In Learning for Dynamics and Control. PMLR, 2021.

Supplementary Material - Taylor-Lagrange Neural Ordinary Differential
Equations: Toward Fast Learning and Evaluation of Neural ODEs

##  Appendix A Proof of Theorem 1

The proof of this theorem is adapted from the results provided in Theorem $2$
of Djeumou et al. [2021].

We assume that the dynamics function $f_{\theta}$ is locally Lipschitz
continuous. Indeed, most common activation functions enforce Lipschitz
continuity of the neural network encoding $f_{\theta}$. Given that the state
$x$ lies in a domain $\mathcal{X}\subseteq\mathbb{R}^{n}$, The Lipschitz
assumption provides that there exists a constant
$L_{f_{\theta}}\in\mathbb{R}^{n}$ such that

$(L_{f_{\theta}})_{k}=\sup\\{L\in\mathbb{R}\>|\>|(f_{\theta})_{k}(x)-(f_{\theta})_{k}(y)|\leq
L\|x-y\|_{2},x,y\in\mathcal{X},x\neq y\\},\;\text{for all }k\in[1,\ldots,n],$

where $(f_{\theta})_{k}$ and $(L_{f_{\theta}})_{k}$ denote the $k$-th
component of the vectors $f_{\theta}$ and $L_{f_{\theta}}$ respectively. We
use such Lipschitz assumption to provide the proof of Theorem 1 in two steps.

First, we provide a bound on the rate of changes of solutions $x(t)$ of the
differential equation $\dot{x}(t)=f_{\theta}(x(t))$.

######  Lemma 1 (Variation of trajectories of $\dot{x}(t)=f_{\theta}(x(t))$).

Let $x:\mathbb{R}_{+}\mapsto\mathcal{X}$ be a continuous-time signal solution
of the differential equation $\dot{x}(t)=f_{\theta}(x(t))$ from the initial
state $x(t_{i})=x_{t_{i}}$. Then, for all $t\in[t_{i},t_{i+1}=t_{i}+\Delta
t]$, we have that

$\displaystyle\|x(t)-x_{t_{i}}\|_{2}\leq\|f_{\theta}(x_{t_{i}})\|_{2}\|L_{f_{\theta}}\|_{2}^{-1}(\mathrm{e}^{\|L_{f_{\theta}}\|_{2}\Delta
t}-1),,$ (6)

###### Proof.

For all $t\in[t_{i},t_{i+1]}$, we have that

$\displaystyle\|x(t)-x_{t_{i}}\|_{2}$
$\displaystyle=\|\int_{t_{i}}^{t}f_{\theta}(x(s))ds\|_{2}$ (7)
$\displaystyle\leq\int_{t_{i}}^{t}\|f_{\theta}(x(s))-f_{\theta}(x_{t_{i}})\|_{2}ds+\int_{t_{i}}^{t}\|f_{\theta}(x_{t_{i}})\|_{2}$
(8)
$\displaystyle\leq\int_{t_{i}}^{t}\|L_{f_{\theta}}\|_{2}\|x(s)-x_{t_{i}}\|_{2}ds+\int_{t_{i}}^{t}\|f_{\theta}(x_{t_{i}})\|_{2}ds$
(9)
$\displaystyle=\|f_{\theta}(x_{t_{i}})\|_{2}(t-t_{i})+\int_{t_{i}}^{t}\|L_{f_{\theta}}\|_{2}\|x(s)-x_{t_{i}}\|_{2}ds.$
(10)

We obtain (7) since $x$ is a solution of $\dot{x}(t)=f_{\theta}(x(t))$. The
passage from (7) to (8) results from applying the triangle inequality. We use
The definition of the Lipschtiz bound to obtain (9). Finally we will use the
Grönwall’s inequality to transform (10) into (6). Specifically, the Grönwall’s
inequality is given as follows. Let $\Psi,\alpha,\gamma$ be real-valued
functions on $[t_{0},T]$ with $t_{0},T\in\mathbb{R}$. Suppose $\Psi$
satisfies, for all $t\in[t_{0},T]$, the inequality

$\displaystyle\Psi(t)\leq\alpha(t)+\int_{t_{0}}^{t}\gamma(s)\Psi(s)ds,$ (11)

with $\gamma(s)\geq 0$ for all $s\in[t_{0},T]$. Then, for all $t\in[t_{0},T]$,

$\displaystyle\Psi(t)\leq\alpha(t)+\int_{t_{0}}^{t}\alpha(s)\gamma(s)\exp\Big{(}\int_{s}^{t}\gamma(r)dr\Big{)}ds.$
(12)

The inequality (10) satisfies the condition (12) of the Grönwall’s inequality
with $\Psi(t)=\|x(t)-x_{t_{i}}\|_{2}$,
$\alpha(t)=\|f_{\theta}(x_{t_{i}})\|_{2}(t-t_{i})$, and
$\gamma(t)=\|L_{f_{\theta}}\|_{2}$. Thus, we have that

$\displaystyle\|x(t)-x_{t_{i}}\|_{2}\leq\begin{aligned}
&\|f_{\theta}(x_{t_{i}})\|_{2}(t-t_{i})+\|L_{f_{\theta}}\|_{2}\|f_{\theta}(x_{t_{i}})\|_{2}\int_{t_{i}}^{t}(s-t_{i})\mathrm{e}^{\|L_{f_{\theta}}\|_{2}(t-s)}ds.\end{aligned}$
(13)

By integration by parts, we have that

$\displaystyle\int_{t_{i}}^{t}(s-t_{i})\mathrm{e}^{\|L_{f_{\theta}}\|_{2}(t-s)}ds$
$\displaystyle=-\frac{t-t_{i}}{\|L_{f_{\theta}}\|_{2}}+\frac{1}{\|L_{f_{\theta}}\|_{2}}\int_{t_{i}}^{t}\mathrm{e}^{\|L_{f_{\theta}}\|_{2}(t-s)}ds=-\frac{t-t_{i}}{\|L_{f_{\theta}}\|_{2}}+\frac{\mathrm{e}^{\|L_{f_{\theta}}\|_{2}(t-t_{i})}-1}{\|L_{f_{\theta}}\|_{2}^{2}}.$
(14)

Finally, by combining (14) and (13), we obtain (6). ∎

Next, we obtain the result in Theorem 1 by exploiting the following closed-
form expression of an a priori enclosure
$\mathcal{S}_{i+1}\subseteq\mathcal{X}$ of the system’s state $x(t_{i+1})$ at
time $t_{i+1}$. That is, $x(t_{i+1})\in\mathcal{S}_{i+1}$. Indeed, the
expression below obtained for $\mathcal{S}_{i+1}$ provide insight on how to
compute the Midpoint in a manner that incorporate information on the
underlying dynamics $f_{\theta}$.

######  Lemma 2 (A priori Enclosure Estimation).

Given the Lipschitz continuous function $f_{\theta}$, the Lipschitz bound
$L_{f_{\theta}}$, and an initial state $x(t_{i})=x_{t_{i}}$, then under the
assumption that $\Delta t\sqrt{n}\|L_{f_{\theta}}\|_{2}<1$, an a priori rough
enclosure $\mathcal{S}_{i+1}\subseteq\mathcal{X}$ of $x(t_{i}+\Delta t)$ is
given by

$\displaystyle x(t_{i}+\Delta t)\in\mathcal{S}_{i}=x_{t_{i}}+[-1,1]^{n\times
n}\frac{\Delta tf_{\theta}(x_{t_{i}})}{1-\sqrt{n}\Delta
t\|L_{f_{\theta}}\|_{2}},$ (15)

where $[-1,1]^{n\times n}$ is a set matrix of size $n\times n$ with elements
between $-1$ and $1$.

###### Proof.

First, it has been proved that a set $\mathcal{S}_{i}$ satisfying the fixed-
point equation $x_{t_{i}}+[0,\Delta
t]\mathrm{Im}(f_{\theta},\mathcal{S}_{i})\subseteq\mathcal{S}_{i}$ is such
that the solution of the differential equation at $t_{i}+\Delta t$ satisfies
$x(t_{i}+\Delta t)\in\mathcal{S}_{i}$. Here, all the operations are set-based
operations, $[0,\Delta t]$ is the interval of values between $0$ and $\Delta
t$, and $\mathrm{Im}(\cdot,\cdot)$ represents the range of a function over a
given domain. In the remainder of this proof, all sets are intervals outer
approximations and the operations between sets are propagated through interval
arithmetic. We refer the reader to the work by Djeumou et al. [2021] for more
details and references on the a priori rough enclosure.

The expression (15) is derived by scaling adequately the bound (6) from Lemma
1. Specifically, we seek for an a priori rough enclosure $\mathcal{S}_{i}$
such that

$\mathcal{S}_{i}=x_{t_{i}}+\mu\frac{\|f_{\theta}(x_{t_{i}})\|_{2}}{\|L_{f_{\theta}}\|_{2}}(\mathrm{e}^{\|L_{f_{\theta}}\|_{2}\Delta
t}-1)[-1,1]^{n},$ (16)

where $\mu>0$ is a parameter to find in order for $\mathcal{S}_{i}$ to satisfy
the fixed-point equation.

We over-approximate the set $\mathrm{Im}(f_{\theta},\mathcal{S}_{i})$ as a
function of $\Delta t$, $\|L_{f_{\theta}}\|_{2}$. Specifically, for all
$s_{i}\in\mathcal{S}_{i}$, we have that

$\displaystyle\|f_{\theta}(s_{i})-f_{\theta}(x_{t_{i}})\|\leq\|L_{f_{\theta}}\|_{2}\|s_{i}-x_{t_{i}}\|\leq\sqrt{n}\mu\|f_{\theta}(x_{t_{i}})\|_{2}(\mathrm{e}^{\|L_{f_{\theta}}\|_{2}\Delta
t}-1).$ (17)

∎

Here, the definition of $\mathcal{S}_{i}$ in (16) provides an upper bound on
$\|s_{i}-x_{i}\|$ that yields (17). Additionally, the inequality (17) implies
that

$\mathrm{Im}(f_{\theta},\mathcal{S}_{i})\subseteq
f_{\theta}(x_{t_{i}})+\sqrt{n}\mu\|f_{\theta}(x_{t_{i}})\|_{2}(\mathrm{e}^{\|L_{f_{\theta}}\|_{2}\Delta
t}-1)[-1,1]^{n}.$

Hence, $\mathcal{S}_{i}$ from (16) solves the fixed-point equation if

$[0,\Delta
t]\big{(}f_{\theta}(x_{t_{i}})+\sqrt{n}\mu\|f_{\theta}(x_{t_{i}})\|_{2}(\mathrm{e}^{\|L_{f_{\theta}}\|_{2}\Delta
t}-1)[-1,1]^{n}\big{)}\subseteq\mu\frac{\|f_{\theta}(x_{t_{i}})\|_{2}}{\|L_{f_{\theta}}\|_{2}}(\mathrm{e}^{\|L_{f_{\theta}}\|_{2}\Delta
t}-1)[-1,1]^{n}.$ (18)

For notation brevity, let
$c_{1}=\|f_{\theta}(x_{t_{i}})\|_{2}(\mathrm{e}^{\|L_{f_{\theta}}\|_{2}\Delta
t}-1)$ and $\beta=\|L_{f_{\theta}}\|_{2}$. Observe that with interval
arithmetic, $[0,\Delta t][a,b]=[\min(0,a),\max(0,b)]$. We use the observation
to find $\mu>0$ such that the inclusion (18) holds. That is, the inequalities

$\displaystyle\Delta t\Big{(}{(f_{\theta}(x_{t_{i}}))}_{k}+\sqrt{n}\mu
c_{1}\Big{)}\leq\frac{\mu c_{1}}{\beta}\Longleftrightarrow(\frac{1}{\Delta
t\beta}-\sqrt{n})\mu\geq\frac{1}{c_{1}}{(f_{\theta}(x_{t_{i}}))}_{k}$

and

$\displaystyle\Delta t\Big{(}{(f_{\theta}(x_{t_{i}}))}_{k}-\sqrt{n}\mu
c_{1}\Big{)}\geq-\frac{\mu
c_{1}}{\beta_{i}}\Longleftrightarrow(\sqrt{n}-\frac{1}{\Delta
t\beta_{i}})\mu\geq\frac{1}{c_{1}}{(f_{\theta}(x_{t_{i}}))}_{k}$

hold for all $k[1,\ldots,n]$. Therefore, for a step size $\Delta t$ satisfying
$\Delta t\sqrt{n}\|L_{f_{\theta}}\|_{2}<1$, $\mu[-1,1]^{n}$ given by

$\displaystyle\mu[-1,1]^{n}=\frac{[-1,1]^{n\times
n}f_{\theta}(x_{t_{i}})}{c_{1}(\frac{1}{\Delta t\beta_{i}}-\sqrt{n})},$

is such that $\mu$ satisfies the above inequalities. Thus, the inclusion (18)
holds and the set $\mathcal{S}_{i}$ is solution of the fixed-point equation.
By replacing $\mu$ in $\mathcal{S}_{i}$ given by (16), we obtain (15).

######  Theorem 1 (Simplified Midpoint Expression).

If $f_{\theta}$ is a Lipschitz-continuous function, then there exists a
function
$\bar{\Gamma}:\mathbb{R}^{n}\times\mathbb{R}_{+}\to\mathbb{R}^{n\times n}$
such that the midpoint value $\Gamma$ of the Taylor-Lagrange expansion
$TL^{\Delta t}_{\theta,\phi}(x_{t_{i}})$ is related to $x_{t_{i}}$,
$f_{\theta}(x_{t_{i}})$, and $\bar{\Gamma}(x_{t_{i}},\Delta t)$ through

$\Gamma=x_{t_{i}}+\bar{\Gamma}(x_{t_{i}},\Delta t)\odot
f_{\theta}(x_{t_{i}}),$ (19)

where $\odot$ denotes matrix-vector multiplication.

###### Proof.

This is a direct application of the result proved in Lemma 2. Specifically,
the midpoint $\Gamma=x_{\epsilon}$ is a point that lies in the rough a priori
enclosure $\mathcal{S}_{i}$, i.e., $\Gamma\ \in\mathcal{S}_{i}$, for some
$\epsilon\in[t_{i},t_{i}+\Delta t]$. By using Lemma 2 and the expression (15),
we can parameterize any point inside $\mathcal{S}_{i}$ as
$x_{t_{i}}+\bar{\Gamma}(x_{t_{i}},\Delta t)\odot f_{\theta}(x_{t_{i}})$. Hence
the results provided by the theorem. ∎

##  Appendix B Simplified Midpoint Expression for Linear Systems

In this section, we provide proof that, for linear systems, the function
$\bar{\Gamma}$ in the simplified midpoint expression (3) does not depend on
the state $x$. We consider linear systems in the form

$\displaystyle\dot{x}(t)=Ax(t),$ (20)

where $x(t)\in\mathbb{R}^{n}$ and $A\in\mathbb{R}^{n\times n}$ is a time-
independent matrix. First, by performing a zero-order Taylor-Lagrange
expansion on the solution $x(t)$ of (20) from an initial point $x(t_{0})$, we
have that

$\displaystyle x(t)$ $\displaystyle=x(t_{0})+(t-t_{0})A\Gamma$ (21)
$\displaystyle=x(t_{0})+(t-t_{0})A\Big{(}x(t_{0})+\bar{\Gamma}(x_{t_{0}},t-t_{0})Ax(t_{0})\Big{)}$
(22)
$\displaystyle=x(t_{0})+(t-t_{0})Ax(t_{0})+(t-t_{0})\Big{(}A\bar{\Gamma}(x_{t_{0}},t-t_{0})A\Big{)}x(t_{0}),$
(23)

where (21) comes from the Taylor expansion and (22) results from the
simplified midpoint expression provided in (3). Next, we consider
$\bar{\Gamma}(x_{t_{0}},t-t_{0})$ given by the state-independent expression

$\displaystyle\bar{\Gamma}(x_{t_{0}},t-t_{0})=\bar{\Gamma}(t-t_{0})=\sum_{i=1}^{\infty}\frac{(t-t_{0})^{i}}{(i+1)!}A^{i-1}.$
(24)

By substituting (24) in (23), we obtain that

$\displaystyle x(t)$
$\displaystyle=x(t_{0})+(t-t_{0})Ax(t_{0})+(t-t_{0})A\Big{(}\sum_{i=1}^{\infty}\frac{(t-t_{0})^{i}}{(i+1)!}\Big{)}Ax(t_{0})$
(25)
$\displaystyle=x(t_{0})+(t-t_{0})Ax(t_{0})+\sum_{i=2}^{\infty}\frac{(t-t_{0})^{i}}{i!}A^{i}x(t_{0})$
(26) $\displaystyle=\exp^{A(t-t_{0})}x(t_{0}),$ (27)

where the right-hand side in the last line equation is traditionally known as
the solution of the linear differential equation (20). Thus, for the linear
differential equation (20), the function $\bar{\Gamma}$ in the simplified
midpoint expression (3) does not depend on $x$

##  Appendix C Proof of Theorem 2

In this section, we prove Theorem 2 on the integration accuracy of the
proposed Taylor-Lagrange expression when the remainder is approximated through
deep neural networks.

######  Theorem 2 (Integration Accuracy).

If the learned midpoint predictor $\Gamma_{\phi}(\cdot)$ is a
$\mathcal{O}(\eta)$ approximator to the midpoint $\Gamma$ of the Taylor-
Lagrange expansion of $TL^{\Delta t}_{\theta,\phi}(x_{t_{i}})$, then
$\|x_{t_{i+1}}-TL^{\Delta t}_{\theta,\phi}(x_{t_{i}})\|\leq c\eta\Delta t^{p}$
for some $c>0$ that depends on $f_{\theta}$.

###### Proof.

The proof of the Theorem follows the classical proof for obtaining truncation
error of Taylor expansion. Specifically, if $\Delta
t^{p}f_{\theta}^{[p]}(\Gamma(x_{t_{i}},\Delta t))$ denotes the unknown
midpoint function we seek to approximate, we have that

$\displaystyle\|x_{t_{i+1}}-TL^{\Delta t}_{\theta,\phi}(x_{t_{i}})\|$
$\displaystyle=\|\Delta t^{p}f_{\theta}^{[p]}(\Gamma_{\phi}(x_{t_{i}},\Delta
t))-\Delta t^{p}f_{\theta}^{[p]}(\Gamma(x_{t_{i}},\Delta t))\|$ (28)
$\displaystyle\leq\Delta
t^{p}L_{f_{\theta}^{[p]}}\|\Gamma_{\phi}(x_{t_{i}},\Delta
t)-\Gamma(x_{t_{i}},\Delta t)\|$ (29) $\displaystyle\leq\Delta
t^{p}L_{f_{\theta}^{[p]}}\eta,$ (30)

where $L_{f_{\theta}^{[p]}}$ is the Lipschitz constant of the $p$-th Taylor
coefficient and $\eta$ is such that $\|\Gamma_{\phi}(x_{t_{i}},\Delta
t)-\Gamma(x_{t_{i}},\Delta t)\|\leq\eta$ by the assumption that
$\Gamma_{\phi}(\cdot)$ is a $\mathcal{O}(\eta)$ approximator to the midpoint
$\Gamma$. ∎

##  Appendix D Additional Experimental Details

##### Code.

All the implementations are written and tested in Python $3.8$, and we will
release the full code upon the paper is accepted. We attach with the initial
submission the latest version of the code with instructions on how to
reproduce each of the results in the paper.

##### Datasets.

We use the MNIST dataset Deng [2012] for the supervised classsification
experiments and the MiniBoONe dataset Roe et al. [2005] for the density
estimation task.

###  D.1 Experiments on a Stiff Dynamical System

In this section, we provide details on the numerical integration experiments
for the stiff dynamics when the vector field is assumed to be known. Then, we
provide additional details for the case where the dynamics are unknown and
must be learned from data.

####  D.1.1 Integrating the Known Stiff Dynamics

We conside a $2$ dimensional linear system with dynamics given by
$\dot{x}=Ax$, where $x\in\mathbb{R}^{2}$ and $A\in\mathbb{R}^{2\times 2}$ has
eigenvalues $\lambda_{1}=-1$ and $\lambda_{2}=-1000$. The matrix $A$ was
chosen to have such a large gap $\lambda_{1}-\lambda_{2}$ in order to obtain a
system with stiff dynamics.

##### Training and Testing Dataset.

The dataset used to train $\bar{\Gamma}_{\phi}$ consists of $100$ separate
trajectories beginning from initial states that are randomly sampled from the
subset $[-0.5,0.5]\times[-0.5,0.5]$. Each trajectory is integrated for $T=10$
seconds, for varying sizes of the integration time step $\Delta t$. The
testing dataset consists of $10$ trajectories over the same duration $T$. By
saying that the trajectories are obtained from integration of the ODE, we mean
that we use the analytic form of the solution to linear dynamics in order to
obtain the testing and training trajectories.

##### TL-NODE Parameterization and Training.

We parameterize the unknown term $\bar{\Gamma}_{\phi}(\cdot)$ in (4) as a
fully connected MLP with relu activation functions and one hidden layer of
size $16$. The network is optimized using ADAM Kingma and Ba [2014] with a
learning rate of $1e^{-3}$ and an exponential decay of $1e^{-4}$ per
iteration. We train the neural network representing the midpoint for $1000$
epochs with a mini batch of size $512$ at each iteration.

##### Hypersolver Parameterization and Training.

We closely follow the formulation in Poli et al. [2020] to implement the
comparison with Hypersolver. Specifically, in this example, we use the
HyperEuler formulation, where the solution of the ODE is given by a first
order Euler approximation and the error of the approximation is given by a
neural network. We parameterize the neural network to learn the approximation
error of the Euler formulation as a fully connected MLP with relu activation
functions and one hidden layer of size $32$. Note that this provides an
advantage to TL-NODE as TL-NODE only uses $16$ nodes in the hidden layer. We
refer the reader to Poli et al. [2020] for more details.

In our experiments, we observe that Hypersolver is more challenging to tune.
For example, the choice of the initial values of the weights significantly
impacts the algorithm convergence rate and the solution at which it converges.
We did our best to tune Hypersolver such that it attains an accuracy close to
the adaptive time step solver. We train the corresponding neural network for
$1000$ epochs using a mini batch of size $512$.

####  D.1.2 Learning the Unknown Stiff Dynamics

We parameterize the unknown stiff dynamics by a two-layered neural network of
size $64$ and $2$, respectively, with no activation functions. That is, the
output of the neural network is a polynomial function of the states. We first
train the vanilla neural ODE approach for $150$ epochs using Adam optimizer
with a mini batch of size $512$. We use a learning rate of $1e-2$ and an
exponential decay of $1e^{-4}$. We use the hyperparameters for the vanilla
NODE to train RK4 and the truncated Taylor method.

We tune our proposed approach (TL-NODE) such that it matches the performance
achieved by the vanilla neural ODE. Specifically, in addition to the dynamics
network $f_{\theta}(\cdot)$, we parameterize the unknown term
$\bar{\Gamma}_{\phi}(\cdot)$ in (3) as a fully connected MLP with relu
activation functions and one hidden layer of size $16$. Then, we use the same
hyperparameters as for the vanilla ODE to update the dynamics network
$f_{\theta}(\cdot)$. We use Adam optimizer with a learning rate of $1e^{-4}$
and exponential decay of $1e^{-4}$ to simultaneously train the dynamics
network $f_{\theta}(\cdot)$ and the midpoint prediction network
$\Gamma_{\phi}(\cdot)$. We choose the parameter $N_{\phi}=200$, which is the
period at which the algorithm switches from updating the dynamics network
parameters $\theta$ to updating the midpoint network parameters $\phi$ (see
Algorithm 2). $N_{\phi}=200$ was picked to optimize the accuracy and
performance of TL-NODE.

###  D.2 Supervised Classification

We train a Neural ODE and a Linear Classifier to map flattened MNIST Images
Deng [2012] to their corresponding labels. Our model uses a two layered neural
network $f_{\theta}$, as the ODE dynamics, followed by a linear classifier
$g_{\theta}$, identical to the architecture used in Kelly et al. [2020]. In
fact, we directly use the code provided by Kelly et al. [2020] to provide the
comparisons shown in this paper.

##### TL-NODE Parameterization and Training

We train our Neural ODE using a batch size of $512$ and for $1000$ epochs
without any early stopping methods. We use Adam optimizer with a linear decay
learning rate from $10^{-3}$ to $10^{-5}$. Additionally, we parameterize the
midpoint value by a single-layer dense neural network with $24$ node. We use
the penalty term value $\lambda=2\cdot 10^{2}$. These values were obtained by
checking how the model converges during the first $10$ epochs.

###  D.3 Density Estimation

For the models trained on the MiniBooNE tabular dataset, we used the same
architecture as in Table 4 in the appendix of Grathwohl et al. [2019]. Our
implementation was also modified from the implementation provided by Kelly et
al. [2020]. The comparisons with TayNODE and RNODE are obtained from the same
architecture and hyper-parameters described in Kelly et al. [2020]. The data
was obtained as made available from Papamakarios et al. [2017], which was
already processed and split into train/validation/test. In particular, the
training set has $29556$ examples, the validation set has $3284$ examples, and
the test set has 3648 examples, which consist of 43 features.

For TL-NODE, we choose the number of epochs and learning rate for both
training the ODE and the midpoint such that we achieve the performance
provided in Kelly et al. [2020]. Specifically, we parameterize the unknown
term $\bar{\Gamma}_{\phi}(\cdot)$ in (3) as a fully connected MLP with relu
activation functions and one hidden layer of size $32$. We train the neural
network encoding the ODE using Adam optimizer with a training batch size of
$1000$ and a number of epochs of $400$. We pick the learning rate to be $1e-3$
for the first $300$ epochs and reduce the learning rate to $1e-5$ for the last
$100$ epochs.

Then, we use Adam optimizer with a learning rate of $1e^{-4}$ and exponential
decay of $1e^{-4}$ to simultaneously train the midpoint value. We choose the
parameter $N_{\phi}=50$, which is the period at which the algorithm switch
from updating the ODE parameters to the midpoint parameters (see Algorithm 2).
Note that in this experiment, we need a more frequent update of the midpoint
value as it helps correcting with respect to the solution of the ODE using an
adaptive time step solver.


# Convergence rate analysis of randomized and cyclic coordinate descent for
convex optimization through semidefinite programming

Hadi Abbaszadehpeivasti1, Etienne de Klerk1,∗, Moslem Zamani1

1Tilburg University, Department of Econometrics and Operations Research,
Tilburg, The Netherlands

Abstract. In this paper, we study randomized and cyclic coordinate descent for
convex unconstrained optimization problems. We improve the known convergence
rates in some cases by using the numerical semidefinite programming
performance estimation method. As a spin-off we provide a method to analyse
the worst-case performance of the Gauss-Seidel iterative method for linear
systems where the coefficient matrix is positive semidefinite with a positive
diagonal.

Keywords. cyclic and randomized coordinate descent; semidefinite programming;
Gauss-Seidel method.

AMS subject classification. 90C25, 90C22

††footnotetext:  ∗Corresponding author. E-mail addresses:
<EMAIL_ADDRESS>(H. Abbaszadehpeivasti),
<EMAIL_ADDRESS>(E. de Klerk),
<EMAIL_ADDRESS>(M. Zamani). Received …; Accepted February ….

##  1\. Introduction

We consider the unconstrained optimization problem

$\displaystyle f^{\star}=\min_{x\in\mathbb{R}^{n}}f(x),$ (1.1)

where $f:\mathbb{R}^{n}\to\mathbb{R}$ is convex. We assume that $f$ attains
its minimum and $f^{\star}$ denotes the optimal value. In addition, we assume
that $f$ is an $L$-smooth function, that is,

$\|\nabla f(y)-\nabla f(x)\|\leq L\|y-x\|,\ \ \ \forall y,x\in\mathbb{R}^{n}.$

Moreover, we denote the component Lipschitz constants by $\ell_{i}$
($i\in\\{1,\ldots,n\\}$), i.e.,

$\displaystyle|[\nabla f(x+te_{i})]_{i}-[\nabla f(x)]_{i}|\leq\ell_{i}|t|,\ \
\ \forall x\in\mathbb{R}^{n},t\in\mathbb{R},$ (1.2)

where $e_{i}$ is the $i$th standard unit vector. Let $\ell_{\max}:=\max_{1\leq
i\leq n}\ell_{i}$, and note that $1\leq\tfrac{L}{\ell_{\max}}\leq n$.

Due to the simplicity and small per-iteration cost, coordinate descent methods
have been employed extensively for large-scale optimization problems [1, 2].

The generic coordinate descent method is shown in Algorithm 1.

Algorithm 1 Generic coordinate descent

Set $N$ and $\\{t_{k}\\}_{k=0}^{N-1}$ (step lengths) and pick
$x^{0}\in\mathbb{R}^{n}$.

For $k=0,1,\ldots,N-1$ perform the following step:

  1. (1)

Choose an index $i_{k}$ from $\\{1,2,...,n\\}$.

  2. (2)

$x^{k+1}=x^{k}-t_{k}[\nabla f(x^{k})]_{i_{k}}e_{i_{k}}$.

In this paper, we revisit the worst-case convergence rate analysis for
Algorithm 1 for two of the best known variants, namely randomized coordinate
descent, and cyclic coordinate descent. In the former, the index $i_{k}$ is
chosen uniformly at random from $\\{1,2,...,n\\}$, and in the latter, the
cyclic ordering is used.

We will improve the best-known convergence rates from the literature for some
specific values of the parameters $n,L,N,t_{k}$ for $k\in\\{0,1,\ldots,N-1\\}$
and $\ell_{i}$ for $i\in\\{1,\ldots,n\\}$. Finally, the Gauss-Seidel iterative
method for positive semidefinite linear systems is a special case cyclic
coordinate descent for convex quadratic functions, and we will investigate the
implications of our analysis for this classical method as well.

Recently, Kamri et al. [3] studied the convergence of the coordinate descent
algorithm using the semidefinite programming (SDP) performance estimation
method, that was introduced by Drori and Teboulle [4]. We will also use SDP
performance estimation in our analysis, and our main contribution may be seen
as the extension and refinement of the approach by Kamri et al. [3]. SDP
performance estimation has been applied to the analysis of many iterative
methods (other than coordinate descent); the interested reader may consult [5,
6, 7, 8, 9] and the references therein. For general background information on
SDP, see e.g. [10].

### Notation and background results

We use $\langle\cdot,\cdot\rangle$ and $\|\cdot\|$ to denote the Euclidean
inner product and norm, respectively, unless indicated otherwise. The column
vector $e_{i}$ represents the $i$-th standard unit vector and $I$ stands for
the identity matrix. For a matrix $A$, $a_{ij}$ denotes its $(i,j)$-th entry,
and $A^{\top}$ represents the transpose of $A$. The function
$f:\mathbb{R}^{n}\to\mathbb{R}$ is called $\mu$-strongly convex function if
the function $x\mapsto f(x)-\tfrac{\mu}{2}\|x\|^{2}$ is convex. Clearly, any
convex function is $0$-strongly convex. We denote the set of real-valued
convex functions which are $L$-smooth and $\mu$-strongly convex by
$\mathcal{F}_{\mu,L}(\mathbb{R}^{n})$.

Let $\mathcal{I}$ be a finite index set and let
$\\{x^{i};g^{i};f^{i}\\}_{i\in\mathcal{I}}\subseteq\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}$.
A set $\\{x^{i};g^{i};f^{i}\\}_{i\in\mathcal{I}}$ is called
$\mathcal{F}_{\mu,L}$-interpolable if there exists
$f\in\mathcal{F}_{\mu,L}(\mathbb{R}^{n})$ with

$f(x^{i})=f^{i},\ g^{i}\in\partial f(x^{i})\ \ i\in\mathcal{I}.$

The next theorem gives necessary and sufficient conditions for
$\mathcal{F}_{\mu,L}$-interpolablity.

######  Theorem 1.1.

[5] Let $L\in(0,\infty)$ and $\mu\in[0,\infty)$ and
$f\in\mathcal{F}_{\mu,L}(\mathbb{R}^{n})$. For any $x,y\in\mathbb{R}^{n}$, we
have

$\displaystyle\tfrac{1}{2(1-\tfrac{\mu}{L})}\left(\tfrac{1}{L}\left\|\nabla
f(x)-\nabla
f(y)\right\|^{2}+\mu\left\|x-y\right\|^{2}-\tfrac{2\mu}{L}\left\langle\nabla
f(y)-\nabla f(x),y-x\right\rangle\right)\leq f(x)-f(y)-\left\langle\nabla
f(y),x-y\right\rangle.$ (1.3)

Conversely, if $\mathcal{I}$ is a finite index set and
$\\{x^{i};g^{i};f^{i}\\}_{i\in\mathcal{I}}\subseteq\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}$
are given data, then the data set is $\mathcal{F}_{\mu,L}$-interpolable if it
satisfies (1.3) in the sense that, for each pair $i,j\in\mathcal{I}$:

$\tfrac{1}{2(1-\tfrac{\mu}{L})}\left(\tfrac{1}{L}\left\|g^{i}-g^{j}\right\|^{2}+\mu\left\|x^{i}-x^{j}\right\|^{2}-\tfrac{2\mu}{L}\left\langle
g^{j}-g^{i},x^{j}-x^{i}\right\rangle\right)\leq f^{i}-f^{j}-\left\langle
g^{j},x^{i}-x^{j}\right\rangle.$

##  2\. Convergence rate of randomized coordinate descent

The randomized coordinate descent method is shown in Algorithm 2 for easy
reference.

Algorithm 2 Randomized coordinate descent

Set $N$ and $\\{t_{k}\\}_{k=0}^{N-1}$ (step lengths) and pick
$x^{0}\in\mathbb{R}^{n}$.

For $k=0,1,\ldots,N-1$ perform the following step:

  1. (1)

Choose index $i_{k}$ with uniform probability from $\\{1,2,...,n\\}$.

  2. (2)

$x^{k+1}=x^{k}-t_{k}[\nabla f(x^{k})]_{i_{k}}e_{i_{k}}$.

We proceed to revisit its worst-case convergence rate for three classes for
function, namely convex $L$-smooth functions, convex quadratic functions, and
strongly convex, $L$-smooth functions.

###  2.1. The case of $L$-smooth functions

Regarding the convergence of Algorithm 2 for $L$-smooth convex functions, the
following is known. (We state the result as in the survey [2, Theorem 1], but
it is originally due to Nesterov [1]).

######  Theorem 2.1.

[2, Theorem 1] Let $f:\mathbb{R}^{n}\to\mathbb{R}$ is an $L$-smooth convex
function for some $L>0$. If $t_{k}=\tfrac{1}{\ell_{\max}}$ for all $k$, then,
for each $k>0$,

$\displaystyle\mathbb{E}\left(f(x^{k})\right)-f^{\star}\leq\left(\tfrac{2n\ell_{\max}}{k}\right)R_{0}^{2},$
(2.1)

where $R_{0}$ satisfies
$\max_{x^{\star}\in\mathbb{S}}\max_{x}\\{\|x-x^{\star}\|:f(x)\leq
f(x^{0})\\}\leq R_{0}$ and $\mathbb{S}$ denotes the optimal solution set.

In this section, we study the behaviour of randomized coordinate descent
method for $L$-smooth convex functions. The worst-case convergence rate of the
Algorithm 2 can be formulated as follows.

$\displaystyle\max$ $\displaystyle\ \mathbb{E}[f(x^{N})]-f(x^{\star})$
$\displaystyle\operatorname*{s.t.}$ $\displaystyle\
f\in\mathcal{F}_{0,L}(\mathbb{R}^{n})$ $f$ satisfies (1.2) for every
$x\in\mathbb{R}^{n}$ for some $\ell_{i}$, $i\in\\{1,\dots,n\\}$ (2.2)
$\displaystyle\|x^{0}-x^{\star}\|^{2}\leq\Delta$ $x^{k}$ $k\in\\{1,2,...,N\\}$
are generated by Algorithm 2 with respect to $x^{0}$ and step length $t_{k}$
$\displaystyle x^{0}\in\mathbb{R}^{n},\ \nabla f(x^{\star})=0,$

where $f,x^{k},x^{\star}$ are decision variables and $t,L,n$ and
$\ell_{i},i\in\\{1,\dots,n\\}$, are the given parameters. Problem (2.1) in
general is intractable. Moreover, note that $x^{k}$ depends on the index
$i_{k}$ which is chosen uniformly at random from the set $\\{1,\dots,n\\}$
therefore (2.1) is a stochastic programming problem. To deal with this we
introduce a random variable $d^{k}$ which depends on the index $i_{k}$ and is
defined by $d^{k}:=[\nabla f(x^{k})]_{i_{k}}e_{i_{k}}$. Note that $d^{k}$ has
the following properties:

$\displaystyle\operatorname{\mathbb{E}}\left[\|d^{k}\|^{2}\right]=\tfrac{1}{n}\operatorname{\mathbb{E}}\left[\|\nabla
f(x^{k})\|^{2}\right]$
$\displaystyle\operatorname{\mathbb{E}}\left[\left\langle d^{k},\nabla
f(x^{k})\right\rangle\right]=\tfrac{1}{n}\operatorname{\mathbb{E}}\left[\|\nabla
f(x^{k})\|^{2}\right]$ (2.3)
$\displaystyle\operatorname{\mathbb{E}}\left[\left\langle
d^{k},x^{k}\right\rangle\right]=\tfrac{1}{n}\operatorname{\mathbb{E}}\left[\left\langle\nabla
f(x^{k}),x^{k}\right\rangle\right],$

where the expectation again refers to the joint distribution of all the random
variables $d^{k}$ for $k\in\\{0,1,\ldots,N\\}$ and $x^{k},\nabla
f(x^{k}),f(x^{k})$ for $k\in\\{0,1,\ldots,N\\}$. By Taylor’s theorem and
(1.2), we have

$\displaystyle f(x^{k+1})\leq f(x^{k})+\left\langle\nabla
f(x^{k}),x^{k+1}-x^{k}\right\rangle+\tfrac{\ell_{\max}}{2}\|x^{k+1}-x^{k}\|^{2}$
(2.4) $\displaystyle f(x^{k})\leq f(x^{k+1})+\left\langle\nabla
f(x^{k+1}),x^{k}-x^{k+1}\right\rangle+\tfrac{\ell_{\max}}{2}\|x^{k}-x^{k+1}\|^{2},$

where, $\ell_{\max}=\max_{i\in\\{1,\dots,n\\}}\ell_{i}$ as before. Therefore,
the relaxation of problem (2.1) is given by

$\displaystyle\max$ $\displaystyle\ \mathbb{E}[f(x^{N})]-f(x^{\star})$
$\displaystyle\operatorname*{s.t.}$ $\displaystyle\ \\{(x^{k};\nabla
f(x^{k});f(x^{k}))\\}\ \textrm{satisfy \eqref{inter.1} for
$k\in\\{0,1,\cdots,N,\star\\}$ w.r.t.}\ \mu=0,L$ $\displaystyle\
\\{(x^{k};\nabla f(x^{k});f(x^{k}))\\}\ \textrm{satisfy \eqref{re.co.lip} for
$k\in\\{0,1,\cdots,N\\}$ w.r.t.}\ \ell_{\max}$ $\displaystyle\ \\{x^{k};\nabla
f(x^{k});d^{k}\\}\ \textrm{satisfies \eqref{d.P}}\;\;\;(k\in\\{0,\cdots,N\\})$
(2.5) $\displaystyle\|x^{0}-x^{\star}\|^{2}\leq\Delta$ $\displaystyle\
x^{k+1}=x^{k}-t_{k}d^{k}$ $\displaystyle x^{0}\in\mathbb{R}^{n},\ \nabla
f(x^{\star})=0,$

where $f(x^{k}),x^{k},x^{\star},\nabla f(x^{k})$ and $d^{k}$ are decision
variables. Note that because the problem (1.1) is invariant under translation,
without loss of generality we may assume that $x^{\star}$ is the zero vector.
Since $x^{k+1}=x^{k}-t_{k}d^{k}$ is a recursive relation, $x^{k}$ can be
written as linear combination of $x^{0}$ and $d^{i}$s. In this way, all the
unknowns appear as entries in the following matrix:

$\displaystyle
G=\begin{pmatrix}\operatorname{\mathbb{E}}[\|x^{0}\|^{2}]&\operatorname{\mathbb{E}}[\langle
x^{0},\nabla f(x^{0})\rangle]&\cdots&\operatorname{\mathbb{E}}[\langle
x^{0},\nabla f(x^{N})\rangle]&\operatorname{\mathbb{E}}[\langle
x^{0},d^{0}\rangle]&\cdots&\operatorname{\mathbb{E}}[\langle
x^{0},d^{N}\rangle]\\\ \operatorname{\mathbb{E}}[\langle\nabla
f(x^{0}),x^{0}\rangle]&\operatorname{\mathbb{E}}[\|\nabla
f(x^{0})\|^{2}]&\cdots&\operatorname{\mathbb{E}}[\langle\nabla f(x^{0}),\nabla
f(x^{N})\rangle]&\operatorname{\mathbb{E}}[\langle\nabla
f(x^{0}),d^{0}\rangle]&\cdots&\operatorname{\mathbb{E}}\langle\nabla
f(x^{0}),d^{N}\rangle\\\ \vdots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
\operatorname{\mathbb{E}}[\langle
d^{N},x^{0}\rangle]&\operatorname{\mathbb{E}}[\langle d^{N},\nabla
f(x^{0})\rangle]&\cdots&\operatorname{\mathbb{E}}[\langle d^{N},\nabla
f(x^{N})\rangle]&\operatorname{\mathbb{E}}[\langle
d^{N},d^{0}\rangle]&\cdots&\operatorname{\mathbb{E}}[\|d^{N}\|^{2}]\\\
\end{pmatrix}.$

Note that $G$ is the expectation of a random Gram matrix. Since every
realization of this random matrix is positive semidefinite, and the
expectation preserves positive semidefiniteness, it follows that $G$ is
positive semidefinite as well. Therefore, problem (2.1) can be written as an
SDP problem, where the variables are $G$ and
$\operatorname{\mathbb{E}}[f(x^{i})]$.

In what follows we compare the convergence rate derived by solving the problem
(2.1) and the bound by Wright (2.1) for some specific values of the parameters
$n,L,N,t_{k}$ for $k\in\\{0,1,\cdots,N\\}$ and $\ell_{i}$ for
$i\in\\{1,\cdots,n\\}$. All the figures in this paper were obtained by solving
the SDP problems with the solver Mosek [11], using the Yalmip [12] Matlab
interface.

$5$$10$$15$$20$$25$$30$$0$$2$$4$Number of iterations
$(N)$$\mathbb{E}[f(x^{N})]-f(x^{\star})$Bound (2.1) ($n=2$)Bound (2.1)
($n=4$)PEP bound (2.1) ($n=2$)PEP bound (2.1) ($n=4$) Figure 1. Convergence
rate for Algorithm 2 computed by performance estimation problem (2.1) (dashed
lines) and the bound given by (2.1) (thick lines) for
$L=2,l_{\max}=1,t=1,\Delta=1$ and different $n$
$5$$10$$15$$20$$25$$30$$0$$2$$4$Number of iterations
$(N)$$\mathbb{E}[f(x^{N})]-f(x^{\star})$Bound (2.1) ($n=2$)Bound (2.1)
($n=4$)PEP bound (2.1) ($n=2$)PEP bound (2.1) ($n=4$) Figure 2. Convergence
rate for Algorithm 2 computed by performance estimation problem (2.1) (dashed
lines) and the bound given by (2.1) (thick lines) for
$L=4,l_{\max}=2,t=0.5,\Delta=1$ and different $n$

Note that the convergence rate provided by solving performance estimation is
strictly better than the bound given by Wright. In other words, the bound
(2.1) is not tight for the values of the parameters that we considered.
Moreover, the bound given by performance estimation can also be calculated for
different step lengths than the fixed step lengths $1/\ell_{\max}$ in the
bound (2.1).

###  2.2. The case of convex quadratic functions

In this section we study the convergence of rate of the randomized coordinate
descent method in case that the objective function is quadratic function of
the form

$\displaystyle\min_{x\in\mathbb{R}^{n}}f(x):=\tfrac{1}{2}x^{\top}Ax-b^{\top}x,$
(2.6)

where, $A$ is a symmetric positive semidefinite matrix. To study this case we
need to add additional constraint to restrict our model to quadratic
functions. The following necessary condition for $f$ to be a quadratic
function can be verified easily, and has been used in SDP performance analysis
by Drori et al [7]:

$\displaystyle\tfrac{1}{2}\langle\nabla f(x)-\nabla
f(y),x-y\rangle=f(x)-f(y)-\langle\nabla f(y),x-y\rangle.$

Since this constraint holds for every point in the domain we just consider the
relaxed constraint that only holds for the point generated by the method in
addition to the initial point and the optimal point. In this case we add the
following constraint to the problem (2.1).

$\displaystyle\tfrac{1}{2}\langle\nabla f(x^{i})-\nabla
f(x^{j}),x^{i}-x^{j}\rangle=f(x^{i})-f(x^{j})-\langle\nabla
f(x^{j}),x^{i}-x^{j}\rangle\ \ \ \forall i,j\in\\{0,\ldots,N,*\\}.$ (2.7)

In what follows we compare the convergence rate of the randomized coordinate
descent method for the general problem (2.1) to the convergence rate for the
quadratic problems.

$5$$10$$15$$20$$25$$30$$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$Number
of iterations $(N)$$\operatorname{\mathbb{E}}[f(x^{N})]-f(x^{\star})$PEP bound
(2.1)PEP bound for quadratic functions Figure 3. Convergence rate for
Algorithm 2 computed by performance estimation problem for quadratic functions
(red line) and the bound given by (2.1) (blue line) for
$n=10,L=2,l_{\max}=1,t=1,\Delta=1$.
$5$$10$$15$$20$$25$$30$$0$$0.1$$0.2$$0.3$$0.4$$0.5$Number of iterations
$(N)$$\operatorname{\mathbb{E}}[f(x^{N})]-f(x^{\star})$PEP bound (2.1)PEP
bound for quadratic functions Figure 4. Convergence rate for Algorithm 2
computed by performance estimation problem for quadratic functions (red line)
and the bound given by (2.1) (blue line) for
$n=2,L=2,l_{\max}=1,t=1,\Delta=1$.

Note that the convergence rate for the quadratic problem is sightly better
than that of the general case.

###  2.3. The case of $\mu$-strongly convex $L$-smooth functions

In this section, we study the convergence rate of the $\mu$-strongly convex
$L$-smooth functions. If $\mu>0$ the optimal value of problem (2.1) for one
iteration of Algorithm 2, i.e. $N=1$, appears to be the same as the following
bound (2.8) given by Wright [2].

######  Theorem 2.2.

[2, Theorem 1] Let $f\in\mathcal{F}_{\mu,L}(\mathbb{R}^{n})$. If
$t_{k}=\tfrac{1}{\ell_{\max}}$ for each $k$, and $\mu>0$, then, for all $N>0$,

$\displaystyle\operatorname{\mathbb{E}}\left[f(x^{N})\right]-f^{\star}\leq\left(1-\tfrac{\mu}{n\ell_{\max}}\right)^{N}\left(f(x^{0})-f^{\star}\right).$
(2.8)

The fact that the two bounds seem to coincide does not suggest that (2.8) is
tight, since the SDP bound (2.1) is a relaxation, and not exact.

##  3\. Cyclic coordinate descent

Cyclic coordinate descent is one of the most important coordinate descent
algorithms due to its simplicity. The convergence rate of cyclic coordinate
descent method for the class of $L$-smooth convex functions is studied by
Kamri et al using the performance estimation method [3]. This method is
described in Algorithm 3.

Algorithm 3 Cyclic coordinate descent

Set number of cycles $K$, $\\{t_{k}\\}_{k=0}^{N-1}$ (step lengths), pick
$x^{0}\in\mathbb{R}^{n}$ and set $N=nK$.

For $k=0,1,2,\ldots,N-1$ perform the following step:

  1. (1)

Set $i=k\pmod{n}+1$

  2. (2)

$x^{k+1}=x^{k}-t_{k}[\nabla f(x^{k})]_{i}e_{i}$.

In each iteration the method updates the current point over one of the
coordinates in cyclic order.

The following result is known about the rate of convergence. We present it as
in [2], but it is originally due to Beck and Tetruashvili [13].

######  Theorem 3.1.

[2, Theorem 3] Let $f:\mathbb{R}^{n}\to\mathbb{R}$ is an $L$-smooth convex
function for some $L>0$. If $t_{k}=\tfrac{1}{\ell_{\max}}$ for all $k$, then,
for $k=n,2n,3n,\ldots$,

$\displaystyle\left(f(x^{k})\right)-f^{\star}\leq\left(\frac{4nR_{0}^{2}\ell_{\max}(1+nL^{2}/\ell^{2}_{\max})}{k+8}\right),$
(3.1)

where $R_{0}$ satisfies
$\max_{x^{\star}\in\mathbb{S}}\max_{x}\\{\|x-x^{\star}\|:f(x)\leq
f(x^{0})\\}\leq R_{0}$ and $\mathbb{S}$ denotes the optimal solution set. If
$f$ is also strongly convex with parameter $\mu>0$, then one has, for
$k=n,2n,3n,\ldots$,

$\left(f(x^{k})\right)-f^{\star}\leq\left(1-\frac{\mu}{2\ell_{\max}(1+nL^{2}/\ell^{2}_{\max})}\right)^{k/n}\left(f(x^{0})-f^{\star}\right).$

For easy reference, we recall the interpolation conditions from Theorem 1.1 in
the case that $\mu=0$: The set $\\{x^{i},\nabla f(x^{i}),f(x^{i})\\}$ for
$i\in\\{0,1,\cdots,N,\star\\}$ is $\mathcal{F}_{0,L}-$ interpolable if and
only if

$\displaystyle f(x^{i})\geq f(x^{j})+\langle\nabla
f(x^{j}),x^{i}-x^{j}\rangle+\tfrac{1}{2L}\left\|\nabla f(x^{i})-\nabla
f(x^{j})\right\|^{2},\ \ \ \forall i,j\in\\{0,1,\cdots,N,\star\\}.$ (3.2)

Using these conditions, we may formulate the worst-case convergence rate as
performance estimation problem.

$\displaystyle\max$ $\displaystyle\ f(x^{N})-f(x^{\star})$
$\displaystyle\operatorname*{s.t.}$ $\displaystyle\ \\{(x^{i};\nabla
f(x^{i});f(x^{i}))\\}\ \textrm{satisfy \eqref{int_glin} for
$i\in\\{0,1,\cdots,N,\star\\}$ w.r.t.}\ L$
$\displaystyle\|x^{0}-x^{\star}\|^{2}\leq\Delta$ (3.3) $\displaystyle\ x^{k}\
\;(k\in\\{1,2,\cdots,N\\})\text{ are generated using Algorithm \ref{Alg-cyc}}$
$\displaystyle x^{0}\in\mathbb{R}^{n},\ \nabla f(x^{\star})=0.$

Problem (3) can be formulated as a semidefinite programming problem, and this
is precisely what was done by Kamri et al. [3].

Since the univariate function $t\mapsto f(x^{k}+te_{i})$ is convex and
$\ell_{i}$-smooth, it follows from (3.2) that, for every two consecutive
points $x^{k}$ and $x^{k+1}$ generated by Algorithm 3, the following
inequalities hold if $i=k\pmod{n}+1$:

$\displaystyle f(x^{k})\geq f(x^{k+1})+\nabla
f(x^{k+1})_{i}(x^{k}_{i}-x^{k+1}_{i})+\tfrac{1}{2\ell_{i}}(\nabla
f(x^{k})_{i}-\nabla f(x^{k+1})_{i})^{2}$ (3.4) $\displaystyle f(x^{k+1})\geq
f(x^{k})+\nabla
f(x^{k})_{i}(x^{k+1}_{i}-x^{k}_{i})+\tfrac{1}{2\ell_{i}}(\nabla
f(x^{k+1})_{i}-\nabla f(x^{k})_{i})^{2}.$

By adding the above inequalities to (3) one can get a better upper bound for
the worst-case convergence rate, i.e.

$\displaystyle\max$ $\displaystyle\ f(x^{N})-f(x^{\star})$
$\displaystyle\operatorname*{s.t.}$ $\displaystyle\ \\{(x^{i};\nabla
f(x^{i});f(x^{i}))\\}\ \textrm{satisfy \eqref{int_glin} for
$i\in\\{1,\cdots,N,\star\\}$ w.r.t.}\ L$ $\displaystyle\ \\{(x^{i};\nabla
f(x^{i});f(x^{i}))\\}\ \textrm{satisfy \eqref{l_ismooth} for
$i\in\\{1,\cdots,N,\star\\}$ w.r.t.}\ \\{\ell_{1},\cdots,\ell_{n}\\}$
$\displaystyle\|x^{0}-x^{\star}\|^{2}\leq\Delta$ (3.5) $\displaystyle\ x^{k}\
\;(k\in\\{1,\cdots,N\\})\text{ are generated using Algorithm \ref{Alg-cyc}}$
$\displaystyle x^{0}\in\mathbb{R}^{n},\ \nabla f(x^{\star})=0.$

In order to obtain an SDP relaxation to (3), we proceed in the same way as
Kamri et al. [3]. We view (3) as a quadratically constrained quadratic program
(QCQP) in variables corresponding to the unknowns

$x^{k}_{i},\;\frac{\partial f(x^{k}_{i})}{\partial x_{i}},\;f(x^{k})\quad
k\in\\{0,1,2,\cdots,N\\},\;i\in\\{1,\ldots,n\\},$

Next we use the following relations to eliminate variables:

$x_{i}^{k+1}-x_{i}^{k}=\left\\{\begin{array}[]{ll}-t_{k}\frac{\partial
f(x^{k}_{i})}{\partial x_{i}}&\mbox{if $i=k\pmod{n}+1$}\\\
0&else\end{array}\right.$

which hold for all $k\in\\{0,1,2,\cdots,N-1\\}$, and $i\in\\{1,\ldots,n\\}$.
Subsequently we form the standard Shor SDP relaxation (see e.g. [14]) of the
resulting QCQP. Note that this is different to the approach we followed for
randomized coordinate descent. In particular, the size of the SDP relaxation
now depends on $n$, which was not the case before. This also limits the
parameter values for which we may solve the SDP relaxations.

In Figure 5 we compare the SDP upper bounds from (3) and (3) for various
parameter values.

The figure shows that the bound can be improved slightly by adding the set of
constraints (3.4) to the model provided by Kamri et al. [3]. Moreover, we add
the constraint which correspond to the quadratic functions (2.7) to the model
(3) which provides us with a better bound for quadratic functions. The
computed values are much better that the theoretical bound (3.1), to the
extent that we do not include this bound in the plot. Indeed, Kamri et al. [3]
already mentioned in their paper that the computed values for their model are
much better than the theoretical bound (3.1).

$5$$10$$15$$20$$25$$30$$0$$0.1$$0.2$$0.3$Number of cycles
$(N/n)$$f(x^{N})-f(x^{\star})$PEP bound (3)PEP bound (3)PEP bound for
quadratic functions Figure 5. Convergence rate Algorithm 3 computed by
performance estimation problem (3) (blue), the bound given by (3) (red) and
the bound for quadratic functions (green) for
$n=2,L=2,\ell_{1}=1,\ell_{2}=1,t=0.5,\Delta=1$.
$5$$10$$15$$0$$0.1$$0.2$$0.3$Number of cycles
$(N/n)$$f(x^{N})-f(x^{\star})$PEP bound of quadratic functions for $n=2$PEP
bound of quadratic functions for $n=3$PEP bound of quadratic functions for
$n=4$ Figure 6. Convergence rate for Algorithm 3 computed by performance
estimation problem of quadratic functions for $\ell_{i}=1\
i\in\\{1,\cdots,n\\},L=\sum_{i=1}^{n}\ell_{i},t=1,\Delta=1$ and different $n$.

Note that our discussion for coordinate-wise cyclic coordinate descent in this
section could be extended to block-wise cyclic coordinate descent in a similar
way as was done by Kamri et al. [3].

###  3.1. Relation to the Gauss-Seidel method

The minimization of the convex quadratic function in (2.6) is equivalent to
the solution of the linear system $Ax=b$. Here, we may assume w.l.o.g. that
$A$ has a positive diagonal. Cyclic coordinate decent for problem (2.6) is
closely related to the iterative Gauss-Seidel method for solving this linear
system. For this reason, cyclic coordinate descent is sometimes also referred
to as nonlinear Gauss-Seidel. It is therefore an interesting question whether
the SDP performance estimation framework yields any new insights on the
performance of the Gauss-Seidel method.

Denoting $A=(a_{ij})$, the iterative Gauss-Seidel method may be described as
follows.

Algorithm 4 Gauss-Seidel method

Set $N$ and pick $x^{0}\in\mathbb{R}^{n}$.

For $k=0,1,\ldots,N-1$ perform the following:

  1. $x_{i}^{k+1}=\frac{1}{a_{ii}}\left(b_{i}-\sum_{j=1}^{i-1}a_{ij}x^{k+1}_{j}-\sum_{j=i+1}^{n}a_{ij}x^{k}_{j}\right)$ ($i=1,\ldots,n$).

This is exactly cyclic coordinate descent with unit step lengths if the
gradient at a point $x$ is replaced by $D^{-1}\nabla f(x)$, where
$f(x)=\frac{1}{2}x^{\top}Ax-b^{\top}x$ as before, and $D$ is the diagonal
matrix with the same diagonal entries as $A$. To see this, recall that the
Fréchet derivative of a differentiable function
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ at a point $x\in\mathbb{R}^{n}$ is the
unique linear operator, say $D_{f}(x):\mathbb{R}^{n}\rightarrow\mathbb{R}$,
such that

$\lim_{\|{h}\|\rightarrow 0}\frac{f({x}+{h})-f({x})-D_{f}(x){h}}{\|{h}\|}=0.$

Once an inner product on $\mathbb{R}^{n}$ is fixed, say
$\langle\cdot,\cdot\rangle$, one may, by the Riesz representation theorem,
express $D_{f}(x){h}=\langle g(x),h\rangle$, where $g(x)$ is called the
gradient vector of $f$ at $x$ with respect to $\langle\cdot,\cdot\rangle$. In
particular, if $\langle\cdot,\cdot\rangle$ is the Euclidean dot product, then
$g(x)=\nabla f(x)$. If one changes to the inner product
$\langle\cdot,\cdot\rangle_{D}$ defined by

$\langle
u,v\rangle_{D}=\sum_{i=1}^{n}{a_{ii}}u_{i}v_{i}\quad\quad(u,v\in\mathbb{R}^{n}),$
(3.6)

then the gradient vector at $x$ becomes $D^{-1}\nabla f(x)$, by the uniqueness
of the Fréchet derivative.

It was shown in [6] that the interpolation condition in Theorem 1.1 holds for
any reference inner product $\langle\cdot,\cdot\rangle$, provided that the
gradient vector is interpreted accordingly.

In other words, the following SDP performance estimation problem gives a bound
on the worst-case performance of the Gauss-Sidel method after $N$ iterations,
when $A$ is a symmetric positive semidefinite matrix with a positive diagonal.

$\displaystyle\max$ $\displaystyle\ f(x^{N})-f(x^{\star})$
$\displaystyle\operatorname*{s.t.}$ $\displaystyle\ \\{(x^{i};\nabla
f(x^{i});f(x^{i}))\\}\ \textrm{satisfy \eqref{int_glin} and \eqref{qua_cons}
for $i\in\\{0,1,\cdots,N,\star\\}$ w.r.t.\ $L=\lambda_{\max}(D^{-1}A)$}$
$\displaystyle\ \\{(x^{i};\nabla f(x^{i});f(x^{i}))\\}\ \textrm{satisfy
\eqref{l_ismooth} for $i\in\\{0,1,\cdots,N,\star\\}$ w.r.t.\
$\ell_{1}=\ldots=\ell_{n}=1$}$
$\displaystyle\|x^{0}-x^{\star}\|^{2}\leq\Delta$ (3.7) $\displaystyle\ x^{k}\
(k\in\\{1,2,\cdots,N\\})\text{ is generated using Algorithm \ref{Alg-cyc}}$
$\displaystyle x^{0}\in\mathbb{R}^{n},\ \nabla f(x^{\star})=0,$

where the inner product is now understood to be the one in (3.6), and the norm
the induced norm for this inner product, and $\lambda_{\max}(D^{-1}A)$ denotes
the largest eigenvalue of $D^{-1}A$. (Note that the eigenvalues of $D^{-1}A$
are real.) Importantly, the reference inner product is not visible in the SDP
performance estimation problem reformulation of (3.1), since only a Gram
matrix for this inner product appears. It is therefore equally valid, for any
inner product, provided that the inner product and norm are interpreted
accordingly. Of course, the Lipschitz constants like (1.2) depend on the norm
as well. It is easy to verify that, for the inner product (3.6), and
$f(x)=\frac{1}{2}x^{\top}Ax-b^{\top}x$, one has $\ell_{1}=\ldots=\ell_{n}=1$
and $L=\lambda_{\max}(D^{-1}A)$ as is used in (3.1).

In summary, we have shown the following.

######  Theorem 3.2.

Consider a solvable system of linear equations $Ax=b$ where $A$ is a symmetric
positive semidefinite matrix with positive diagonal, and let $x^{\star}$
denote a solution. Letting $f(x)=\frac{1}{2}x^{\top}Ax-b^{\top}x$, after $N$
iterations of the Gauss-Seidel method, an upper bound on
$f(x^{N})-f(x^{\star})$ is given by the optimal value of the SDP problem
(3.1), provided that the starting point $x^{0}$ satisfies
$\|x^{0}-x^{\star}\|\leq\Delta$ for a given $\Delta$, where the norm is the
induced norm of the inner product (3.6).

The Gauss-Seidel method is known to be convergent when $A$ is symmetric
positive-definite, e.g. [15, Theorem 10.1.2], or strictly or irreducibly
diagonally dominant, e.g. [16]. The case when $A$ is only positive
semidefinite (with positive diagonal) seems to be less well-understood, and
our approach sheds more light on this case. In particular, numerical results
of the type shown in Figure 6 apply here.

##  4\. Conclusion

We have studied SDP performance estimation approaches to analyse randomized
and cyclic coordinate descent, thereby complementing recent results in [3].
For randomized coordinate descent, we have given the first known SDP
performance estimation bound. For cyclic coordinate descent, we were able to
improve slightly on the numerical values given in [3], and we also discussed
the link with the Gauss-Seidel method in the case of convex quadratic
functions. Of course, to obtain new rates of convergence in general, it is
necessary to solve the SDP performance estimation problems analytically, as
opposed to numerically, but we have been unable to obtain analytic solutions
for the SDP problems presented in this paper.

Acknowledgments

This work was supported by the Dutch Scientific Council (NWO) grant
OCENW.GROOT.2019.015, _Optimization for and with Machine Learning (OPTIMAL)_.

## References

  * [1] Yu Nesterov.  Efficiency of coordinate descent methods on huge-scale optimization problems.  SIAM Journal on Optimization, 22(2):341–362, 2012.
  * [2] Stephen J Wright.  Coordinate descent algorithms.  Mathematical Programming, 151(1):3–34, 2015.
  * [3] Yassine Kamri, Julien M Hendrickx, and François Glineur.  On the worst-case analysis of cyclic coordinate-wise algorithms on smooth convex functions.  arXiv preprint arXiv:2211.17018, 2022.
  * [4] Yoel Drori and Marc Teboulle.  Performance of first-order methods for smooth convex minimization: a novel approach.  Mathematical Programming, 145(1):451–482, 2014.
  * [5] Adrien B Taylor, Julien M Hendrickx, and François Glineur.  Smooth strongly convex interpolation and exact worst-case performance of first-order methods.  Mathematical Programming, 161(1):307–345, 2017.
  * [6] Etienne De Klerk, François Glineur, and Adrien B Taylor.  Worst-case convergence analysis of inexact gradient and Newton methods through semidefinite programming performance estimation.  SIAM Journal on Optimization, 30(3):2053–2082, 2020.
  * [7] Yoel Drori and Ohad Shamir.  The complexity of finding stationary points with stochastic gradient descent.  In International Conference on Machine Learning, pages 2658–2667. PMLR, 2020.
  * [8] Hadi Abbaszadehpeivasti, Etienne de Klerk, and Moslem Zamani.  On the rate of convergence of the difference-of-convex algorithm (DCA).  arXiv preprint arXiv:2109.13566, 2021.
  * [9] Teodor Rotaru, François Glineur, and Patrinos Panagiotis.  Tight convergence rates of the gradient method on hypoconvex functions.  arXiv preprint arXiv:2203.00775, 2022.
  * [10] Henry Wolkowicz, Romesh Saigal, and Lieven Vandenberghe.  Handbook of semidefinite programming: theory, algorithms, and applications, volume 27.  Springer Science & Business Media, 2012.
  * [11] MOSEK ApS.  The MOSEK optimization toolbox for MATLAB manual. Version 9.0., 2019.
  * [12] J. Löfberg.  YALMIP : A toolbox for modeling and optimization in MATLAB.  In In Proceedings of the CACSD Conference, Taipei, Taiwan, 2004\.
  * [13] Amir Beck and Luba Tetruashvili.  On the convergence of block coordinate descent type methods.  SIAM Journal on Optimization, 23(4):2037–2060, 2013.
  * [14] Alex L. Wang and Fatma Kilinc-Karzan.  A geometric view of SDP exactness in QCQPs and its applications.  arXiv preprint arXiv:2011.07155, 2020.
  * [15] G.H. Golub and C.F. Van Loan.  Matrix computations.  The John Hopkins University Press, Baltimore, 2nd edition, 1989.
  * [16] Roberto Bagnara.  A unified proof for the convergence of Jacobi and Gauss–Seidel methods.  SIAM Review, 37(1):93–97, 1995.