# Efficient and robust high-dimensional sparse logistic regression via
nonlinear primal-dual hybrid gradient algorithms

Jérôme Darbon and Gabriel P. Langlois

###### Abstract

Logistic regression is a widely used statistical model to describe the
relationship between a binary response variable and predictor variables in
data sets. It is often used in machine learning to identify important
predictor variables. This task, variable selection, typically amounts to
fitting a logistic regression model regularized by a convex combination of
$\ell_{1}$ and $\ell_{2}^{2}$ penalties. Since modern big data sets can
contain hundreds of thousands to billions of predictor variables, variable
selection methods depend on efficient and robust optimization algorithms to
perform well. State-of-the-art algorithms for variable selection, however,
were not traditionally designed to handle big data sets; they either scale
poorly in size or are prone to produce unreliable numerical results. It
therefore remains challenging to perform variable selection on big data sets
without access to adequate and costly computational resources. In this paper,
we propose a nonlinear primal-dual algorithm that addresses these
shortcomings. Specifically, we propose an iterative algorithm that provably
computes a solution to a logistic regression problem regularized by an elastic
net penalty in $O(T(m,n)\log(1/\epsilon))$ operations, where
$\epsilon\in(0,1)$ denotes the tolerance and $T(m,n)$ denotes the number of
arithmetic operations required to perform matrix-vector multiplication on a
data set with $m$ samples each comprising $n$ features. This result improves
on the known complexity bound of $O(\min(m^{2}n,mn^{2})\log(1/\epsilon))$ for
first-order optimization methods such as the classic primal-dual hybrid
gradient or forward-backward splitting methods.

#### Significance statement

Logistic regression is a widely used statistical model to describe the
relationship between a binary response variable and predictor variables in
data sets. With the trends in big data, logistic regression is now commonly
applied to data sets whose predictor variables range from hundreds of
thousands to billions. State-of-the-art algorithms for fitting logistic
regression models, however, were not traditionally designed to handle big data
sets; they either scale poorly in size or are prone to produce unreliable
numerical results. This paper proposes a nonlinear primal-dual algorithm that
provably computes a solution to a logistic regression problem regularized by
an elastic net penalty in $O(T(m,n)\log(1/\epsilon))$ operations, where
$\epsilon\in(0,1)$ denotes the tolerance and $T(m,n)$ denotes the number of
arithmetic operations required to perform matrix-vector multiplication on a
data set with $m$ samples each comprising $n$ features. This result improves
on the known complexity bound of $O(\min(m^{2}n,mn^{2})\log(1/\epsilon))$ for
first-order optimization methods such as the classic primal-dual hybrid
gradient or forward-backward splitting methods.

##  1 Introduction

Logistic regression is a widely used statistical model to describe the
relationship between a binary response variable and predictor variables in
data sets [36]. It is often used in machine learning to identify important
predictor variables [22, 76]. This task, variable selection, typically amounts
to fitting a logistic regression model regularized by a convex combination of
$\ell_{1}$ and $\ell_{2}^{2}$ penalties. Variable selection is frequently
applied to problems in medicine [2, 9, 30, 52, 56, 71, 77], natural language
processing [6, 48, 29, 55, 62], economics [46, 64, 74, 75], and social science
[1, 39, 50], among others.

Since modern big data sets can contain up to billions of predictor variables,
variable selection methods require efficient and robust optimization
algorithms to perform well [47]. State-of-the-art algorithms for variable
selection methods, however, were not traditionally designed to handle big data
sets; they either scale poorly in size [14] or are prone to produce unreliable
numerical results [8, 45, 72, 73]. These shortcomings in terms of efficiency
and robustness make variable selection methods on big data sets essentially
impossible without access to adequate and costly computational resources [18,
57]. Further exacerbating this problem is that machine learning applications
to big data increasingly rely on computing power to make progress [19, 37, 47,
42]. Without efficient and robust algorithms to minimize monetary and energy
costs, these shortcomings prevent scientific discoveries. Indeed, it is
expected that progress will rapidly become economically and environmentally
unsustainable as computational requirements become a severe constraint [65].

This paper proposes a novel optimization algorithm that addresses the
shortcomings of state-of-the-art algorithms used for variable selection. Our
proposed algorithm is an accelerated nonlinear variant of the classic primal-
dual hybrid gradient (PDHG) algorithm, a first-order optimization method
initially developed to solve imaging problems [23, 54, 78, 12, 35, 13]. Our
proposed accelerated nonlinear PDHG algorithm, which is based on the work the
authors recently provided in [16], uses the Kullback–Leibler divergence to
efficiently fit a logistic regression model regularized by a convex
combination of $\ell_{1}$ and $\ell_{2}^{2}$ penalties. Specifically, our
algorithm provably computes a solution to a logistic regression problem
regularized by an elastic net penalty in $O(T(m,n)\log(1/\epsilon))$
operations, where $\epsilon\in(0,1)$ denotes the tolerance and $T(m,n)$
denotes the number of arithmetic operations required to perform matrix-vector
multiplication on a data set with $m$ samples each comprising $n$ features.
This result improves on the known complexity bound of
$O(\min(m^{2}n,mn^{2})\log(1/\epsilon))$ for first-order optimization methods
such as the classic primal-dual hybrid gradient or forward-backward splitting
methods.

### Organization of this preprint

In Section 2, we describe how variable selection works with logistic
regression regularized by the elastic net penalty, why this problem is
challenging, what the state-of-the-art algorithms are, and what their
limitations are. In Section 3, we describe our approach for solving this
problem using the Kullback–Leibler divergence, we derive an explicit algorithm
for solving this problem, and we explain why our algorithm overcomes the
limitations of current state-of-the-art algorithms. We also describe how our
approach can be adapted to solve a broad class logistic regression problems
regularized by an appropriate convex penalty, including, for example, the
Lasso penalty. Finally, Section 4 provides a detailed derivation of the
explicit algorithms described in Section 3.

##  2 Preliminaries

### Description of the problem

Suppose we receive $m$ independent samples
$\\{(\boldsymbol{x}_{i},y_{i})\\}_{i=1}^{m}$, each comprising an
$n$-dimensional vector of predictor variables
$\boldsymbol{x}_{i}\in\mathbb{R}^{n}$ and a binary response variable
$y_{i}\in\\{0,1\\}$. The predictor variables are encoded in an $m\times n$
matrix $\bm{A}$ whose rows are the vectors
$\boldsymbol{x}_{i}=(x_{i1},\dots,x_{in})$, and the binary response variables
are encoded in an $m$-dimensional vector $\boldsymbol{y}$. The goal of
variable selection is to identify which of the $n$ predictor variables best
describe the $m$ response variables. A common approach to do so is to fit a
logistic regression model regularized by a convex combination of $\ell_{1}$
and $\ell_{2}^{2}$ penalties:

$\inf_{\boldsymbol{\theta}\in\mathbb{R}^{n}}f(\boldsymbol{\theta};\alpha,\lambda)=\inf_{\boldsymbol{\theta}\in\mathbb{R}^{n}}\left\\{\frac{1}{m}\sum_{i=1}^{m}\log\left(1+\exp{((\bm{A}\boldsymbol{\theta})_{i})}\right)-\frac{1}{m}\left\langle\boldsymbol{y},\bm{A}\boldsymbol{\theta}\right\rangle+\lambda\left(\alpha\left\|{\boldsymbol{\theta}}\right\|_{1}+\frac{1-\alpha}{2}\left\|{\boldsymbol{\theta}}\right\|_{2}^{2}\right)\right\\},$
(1)

where $\lambda>0$ is a tuning parameter and $\alpha\in(0,1)$ is a fixed
hyperparameter. The function
$\boldsymbol{\theta}\mapsto\lambda(\alpha\left\|{\boldsymbol{\theta}}\right\|_{1}+(1-\alpha)\left\|{\boldsymbol{\theta}}\right\|_{2}^{2}/2)$
is called the elastic net penalty [79]. It is a compromise between the ridge
penalty ($\alpha=0$) [34] and the lasso penalty ($\alpha=1$) [66]. The choice
of $\alpha$ depends on the desired prediction model; for variable selection
its value is often chosen to be close to but not equal to one [63].

The elastic net regularizes the logistic regression model in three ways.
First, it ensures that the logistic regression problem (1) has a unique
solution (global minimum) [21, Chapter II, Proposition 1.2]. Second, the
$\ell_{2}^{2}$ penalty shrinks the coefficients of correlated predictor
variables toward each other (and zero), which alleviates negative correlation
effects (e.g., high variance) between highly correlated predictor variables.
Third, the $\ell_{1}$ penalty promotes sparsity in the solution of (1); that
is, the global minimum of (1) has a number of entries that are identically
zero [26, 22, 76]. We note that other penalties are sometimes used in practice
to promote sparsity, including, for example, the group lasso penalty [49]. In
any case, the non-zero entries are identified as the important predictor
variables, and the zero entries are discarded. The number of non-zero entries
itself depends on the value of the fixed hyperparameter $\alpha$ and the
tuning parameter $\lambda$.

In most applications, the desired value of $\lambda$ proves challenging to
estimate. To determine an appropriate value for it, variable selection methods
first compute a sequence of minimums $\boldsymbol{\theta}^{*}(\lambda)$ of
problem (1) from a chosen sequence of values of the parameter $\lambda$ and
then choose the parameter that gives the preferred minimum [8, 28]. Variable
selection methods differ in how they choose the sequence of parameters
$\lambda$ and how they repeatedly compute global minimums of problem (1), but
the procedure is generally the same. The sequence of parameters thus computed
is called a regularization path [28].

Unfortunately, computing a regularization path to problem (1) can be
prohibitively expensive for big data sets. To see why, fix $\alpha\in(0,1)$
and $\lambda>0$, and let
$\boldsymbol{\theta}_{\epsilon}(\alpha,\lambda)\in\mathbb{R}^{n}$ with
$\epsilon>0$ denote an $\epsilon$-approximate solution to the true global
minimum $\boldsymbol{\theta}^{*}(\alpha,\lambda)$ in (1), i.e.,

$f(\boldsymbol{\theta}_{\epsilon}(\alpha,\lambda);\alpha,\lambda)-f(\boldsymbol{\theta}^{*}(\alpha,\lambda);\alpha,\lambda)<\epsilon.$

Then the best achievable rate of convergence for computing
$\boldsymbol{\theta}_{\epsilon}(\alpha,\lambda)$ in the Nesterov class of
optimal first-order methods is linear, that is, $O(\log(1/\epsilon))$ in the
number of iterations [51]. While optimal, this rate of convergence is
difficult to achieve in practice because it requires a precise estimate of the
largest singular value of the matrix $\bm{A}$, a quantity essentially
impossible to compute for large matrices due to its prohibitive computational
cost of $O(\min{(m^{2}n,mn^{2})})$ operations [31]. This issue generally makes
solving problem (1) difficult and laborious. As computing a regularization
path entails repeatedly solving problem (1) for different values of $\lambda$,
this process can become particularly time consuming and resource intensive for
big data sets.

In summary, variable selection methods work by repeatedly solving an
optimization problem that can be prohibitively computationally expensive for
big data sets. This issue has driven much research in the development of
robust and efficient algorithms to minimize costs and maximize performance.

### Algorithms for variable selection methods and their shortcomings

The state of the art for computing regularization paths to problem (1) is
based on coordinate descent algorithms [27, 28, 32, 60, 61, 67, 70, 73]. These
algorithms are implemented, for example, in the popular glmnet software
package [32], which is available in the Python, MATLAB, and R programming
languages. Other widely used variable selection methods include those based on
the least angle regression algorithm and its variants [20, 33, 41, 68, 79],
and those based on the forward-backward splitting algorithm and its variants
[5, 11, 17, 58, 59]. Here, we focus on these algorithms, but before doing so
we wish to stress that many more algorithms have been developed to compute
minimums of (1); see [7, 22, 43, 69, 76] for recent surveys and comparisons of
different methods and models.

Coordinate descent algorithms are considered the state of the art because they
are scalable, with steps in the algorithms generally having an asymptotic
space complexity of at most $O(mn)$ operations. Some coordinate descent
algorithms, such as those implemented in the glmnet software [32], also offer
options for parallel computing. Despite these advantages, coordinate descent
algorithms generally lack robustness and good convergence properties. For
example, the glmnet implementation depends on the sparsity of the matrix
$\bm{A}$ to converge fast [79], and it is known to be slowed down when the
predictor variables are highly correlated [27]. This situation often occurs in
practice, and it would be desirable to have a fast algorithm for this case.
Another issue is that the glmnet implementation approximates the logarithm
term in problem (1) with a quadratic in order to solve the problem
efficiently. Without costly step-size optimization, which glmnet avoids to
improve performance, the glmnet implementation may not converge [28, 41]. Case
in point, Yuan et al. [72] provides two numerical experiments in which glmnet
does not converge. Although some coordinate descent algorithms recently
proposed in [10] and in [24] can provably solve the logistic regression
problem (1) (with parameter $\alpha=1$), in the first case, the convergence
rate is strictly less than the achievable rate, and in the second case, the
method fails to construct meaningful regularization paths to problem (1), in
addition to having large memory requirements.

The least angle regression algorithm is another popular tool for computing
regularization paths to problem (1). This algorithm, however, scales poorly
with the size of data sets because the entire sequence of steps for computing
regularization paths has an asymptotic space complexity of at most
$O(\min{(m^{2}n+m^{3},mn^{2}+n^{3})})$ operations [20]. It also lacks
robustness because, under certain conditions, it fails to compute meaningful
regularization paths to problem (1) [8, 45]. Case in point, Bringmann et al.
[8] provides an example for which the least angle regression algorithm fails
to converge.

The forward-backward splitting algorithm and its variants are widely used
because they are robust and can provably compute $\epsilon$-approximate
solutions of (1) in at most $O(\log(1/\epsilon))$ iterations. To achieve this
convergence rate, the step size parameter in the algorithm needs to be fine-
tuned using a precise estimate of the largest singular value of the matrix
$\bm{A}$. As mentioned before, however, computing this estimate is essentially
impossible for large matrices due to its prohibitive computational cost, which
has an asymptotic computational complexity of at most
$O(\min{(m^{2}n,mn^{2})})$ operations. Line search methods and other
heuristics are often employed to bypass this problem, but they come at the
cost of slowing down the convergence of the forward-backward splitting
algorithm. Another approach is to compute a crude estimate of the largest
singular value of the matrix $\bm{A}$, but doing so dramatically reduces the
speed of convergence of the algorithm. This problem makes regularization path
construction methods based on the forward-backward splitting algorithm and its
variants generally inefficient and impractical for big data sets.

In summary, state-of-the-art and other widely used variable selection methods
for computing regularization paths to problem (1) either scale poorly in size
or are prone to produce unreliable numerical results. These shortcomings in
terms of efficiency and robustness make it challenging to perform variable
selection on big data sets without access to adequate and costly computational
resources. This paper proposes an efficient and robust optimization algorithm
for solving problem (1) that addresses these shortcomings.

##  3 Methodology

We consider the problem of solving the logistic regression problem (1) with
$\alpha\in(0,1)$. Our approach is to reformulate problem (1) as a saddle-point
problem and solve the latter using an appropriate primal-dual algorithm. Based
on work the authors recently provided in [16], we propose to use a nonlinear
PDHG algorithm with Bregman divergence terms tailored to the logistic
regression model and the elastic net penalty in (1). Specifically, we propose
to use the Bregman divergence generated from the negative sum of $m$ binary
entropy functions. This divergence is the function
$D_{H}\colon\mathbb{R}^{m}\times\mathbb{R}^{m}\to[0,+\infty]$ given by

$D_{H}(\boldsymbol{s},\boldsymbol{s}^{\prime})=\begin{dcases}&\sum_{i=1}^{m}s_{i}\log\left(\frac{s_{i}}{s_{i}^{\prime}}\right)+(1-s_{i})\log\left(\frac{1-s_{i}}{1-s_{i}^{\prime}}\right)\quad\mathrm{if}\,\boldsymbol{s},\boldsymbol{s}^{\prime}\in[0,1]^{m},\\\
&+\infty,\quad\mathrm{otherwise}.\end{dcases}$ (2)

We also show how to adapt our approach for solving the logistic regression
problem (1) with the lasso penalty ($\alpha=0$) and a broad class of convex
penalties, including for example the group lasso.

### Numerical optimization algorithm

The starting point of our approach is to express the logistic regression
problem (1) in saddle-point form. To do so, we use the convex conjugate
formula of the sum of logarithms that appears in (1), namely

$\psi(\boldsymbol{s})=\sup_{\boldsymbol{u}\in\mathbb{R}^{n}}\left\\{\left\langle\boldsymbol{s},\boldsymbol{u}\right\rangle-\sum_{i=1}^{m}\log(1+\exp{(u_{i})})\right\\}=\begin{dcases}&\sum_{i=1}^{m}s_{i}\log(s_{i})+(1-s_{i})\log(1-s_{i})\quad\mathrm{if}\,\boldsymbol{s}\in[0,1]^{m},\\\
&+\infty,\quad\mathrm{otherwise}.\end{dcases}$ (3)

Hence we have the representation

$\sum_{i=1}^{m}\log(1+\exp{((\bm{A}\boldsymbol{\theta})_{i})})=\sup_{\boldsymbol{s}\in[0,1]^{m}}\left\\{\left\langle\boldsymbol{s},\bm{A}\boldsymbol{\theta}\right\rangle-\psi(\boldsymbol{s})\right\\},$

and from it we can express problem (1) in saddle-point form as

$\inf_{\boldsymbol{\theta}\in\mathbb{R}^{n}}\sup_{\boldsymbol{s}\in[0,1]^{m}}\left\\{-\frac{1}{m}\psi(\boldsymbol{s})-\frac{1}{m}\left\langle\boldsymbol{y}-\boldsymbol{s},\bm{A}\boldsymbol{\theta}\right\rangle+\lambda\left(\alpha\left\|{\boldsymbol{\theta}}\right\|_{1}+\frac{1-\alpha}{2}\left\|{\boldsymbol{\theta}}\right\|_{2}^{2}\right)\right\\}.$
(4)

A solution to the convex-concave saddle-point problem (4) is called a saddle
point. For $\alpha\in(0,1)$, the saddle-point problem (1) has a unique saddle
point $(\boldsymbol{\theta}^{*},\boldsymbol{s}^{*})$, where the element
$\boldsymbol{\theta}^{*}$ itself is the unique global minimum of the original
problem (1) [21, Proposition 3.1, page 57]. Hence for our purpose it suffices
to compute a solution to the saddle-point problem (4), and to do so we can
take advantage of the fact that the saddle point
$(\boldsymbol{\theta}^{*},\boldsymbol{s}^{*})$ satisfies the following
optimality conditions:

$\frac{1}{m}\bm{A}^{T}(\boldsymbol{y}-\boldsymbol{s}^{*})-\lambda(1-\alpha)\boldsymbol{\theta}^{*}\in\lambda\alpha\partial\left\|{\boldsymbol{\theta}^{*}}\right\|_{1}\quad\mathrm{and}\quad
s_{i}^{*}=\frac{1}{1+\exp{(-(\bm{A}\theta^{*})_{i})}}\quad\mathrm{for}\,i\in\\{1,\dots,m\\}.$
(5)

The next step of our approach is to split the infimum and supremum problems in
(4) with an appropriate primal-dual scheme. We propose to alternate between a
nonlinear proximal ascent step using the Kullback–Leibler divergence (2) and a
proximal descent step using a quadratic function:

$\displaystyle\boldsymbol{s}^{(k+1)}=\operatorname*{arg\,max}_{\boldsymbol{s}\in(0,1)^{m}}\left\\{-\psi(\boldsymbol{s})+\left\langle\boldsymbol{s},\bm{A}(\boldsymbol{\theta}^{(k)}+\rho(\boldsymbol{\theta}^{(k)}-\boldsymbol{\theta}^{(k-1)}))\right\rangle-\frac{1}{\sigma}D_{H}(\boldsymbol{s},\boldsymbol{s}^{(k)})\right\\},$
(6)
$\displaystyle\boldsymbol{\theta}^{(k+1)}=\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\mathbb{R}^{n}}\left\\{\left(\lambda_{1}\left\|{\boldsymbol{\theta}}\right\|_{1}+\frac{\lambda_{2}}{2}\left\|{\boldsymbol{\theta}}\right\|_{2}^{2}\right)+\left\langle\boldsymbol{s}^{(k+1)}-\boldsymbol{y},\bm{A}\boldsymbol{\theta}\right\rangle+\frac{1}{2\tau}\left\|{\boldsymbol{\theta}-\boldsymbol{\theta}^{(k)}}\right\|_{2}^{2}\right\\},$

where $\lambda_{1}=m\lambda\alpha$, $\lambda_{2}=m\lambda(1-\alpha)$, and
$\rho,\sigma,\tau>0$ are parameters to be specified in the next step of our
approach. The scheme starts from initial values
$\boldsymbol{s}^{(0)}\in(0,1)^{m}$ and
$\boldsymbol{\theta}^{(-1)}=\boldsymbol{\theta}^{(0)}\in\mathbb{R}^{n}$.

The key element in this primal-dual scheme is the choice of the
Kullback–Leibler divergence (2) in the first line of (LABEL:eq:kl-nPDHG-alg).
Its choice is motivated by two facts. First, because it is generated from the
sum of $m$ binary entropy that appears _explicitly_ in the saddle-point
problem (4) as the function $\psi$ defined in (3), i.e.,

$D_{H}(\boldsymbol{s},\boldsymbol{s}^{\prime})=\psi(\boldsymbol{s})-\psi(\boldsymbol{s}^{\prime})-\left\langle\boldsymbol{s}-\boldsymbol{s}^{\prime},\nabla\psi(\boldsymbol{s}^{\prime})\right\rangle.$
(7)

This fact will make the maximization step in (LABEL:eq:kl-nPDHG-alg) easy to
evaluate. Second, because it is strongly convex with respect to the $\ell
1$-norm in that

$D_{H}(\boldsymbol{s},\boldsymbol{s}^{\prime})\geqslant\frac{1}{2}\left\|{\boldsymbol{s}-\boldsymbol{s}^{\prime}}\right\|_{1}^{2}$

for every $\boldsymbol{s},\boldsymbol{s}^{\prime}\in[0,1]^{m}$, which is a
direct consequence of a fundamental result in information theory known as
Pinsker’s inequality [4, 15, 38, 40, 53].

The latter fact, notably, implies that the primal-dual scheme (LABEL:eq:kl-
nPDHG-alg) alternates between solving a 1-strongly concave problem over the
space $(\mathbb{R}^{m},\left\|{\cdot}\right\|_{1})$ and a
$\lambda_{2}$-strongly convex problem over the space
$(\mathbb{R}^{n},\left\|{\cdot}\right\|_{2})$. The choice of these spaces is
significant, for it induces the matrix norm

$\left\|{\bm{A}}\right\|_{op}=\sup_{\left\|{\boldsymbol{s}}\right\|_{1}=1}\left\|{\bm{A}^{T}\boldsymbol{s}}\right\|_{2}=\max_{i\in\\{1,\dots,m\\}}\sqrt{\sum_{j=1}^{n}{A_{ij}^{2}}}=\max_{i\in\\{1,\dots,m\\}}\left\|{\boldsymbol{x}_{i}}\right\|_{2},$
(8)

which can be computed in _optimal_ $\Theta(mn)$ time. This is unlike most
first-order optimization methods, such as the forward-backward splitting
algorithm, where instead the matrix norm is the largest singular value of the
matrix $\bm{A}$, which takes $O(\min{(m^{2}n,mn^{2})})$ operations to compute.
This point is _crucial_ : the smaller computational cost makes it easy and
efficient to estimate all the parameters of the nonlinear PDHG algorithm,
which is needed to achieve an optimal rate of convergence.

The last step of our approach is to choose the parameters $\rho$, $\sigma$,
and $\tau$ so that the iterations in the primal-dual scheme (LABEL:eq:kl-
nPDHG-alg) converge. Based on the analysis of accelerated nonlinear PDHG
algorithms the authors recently provided in [16, Section 5.4], the choice of
parameters

$\rho=1-\frac{\lambda_{2}}{2\left\|{\bm{A}}\right\|_{op}^{2}}\left(\sqrt{1+\frac{4\left\|{\bm{A}}\right\|_{op}^{2}}{\lambda_{2}}}-1\right),\quad\sigma=\frac{1-\rho}{\rho},\quad\mathrm{and}\quad\tau=\frac{(1-\rho)}{\lambda_{2}\rho},$

ensure that the iterations converge to the unique saddle point
$(\boldsymbol{\theta}^{*},\boldsymbol{s}^{*})$ of problem (4). In particular,
the rate of convergence is linear in the number of iterations, with

$\frac{1}{2}\left\|{\boldsymbol{\theta}^{*}-\boldsymbol{\theta}^{(k)}}\right\|_{2}^{2}\leqslant\rho^{k}\left(\frac{1}{2}\left\|{\boldsymbol{\theta}^{*}-\boldsymbol{\theta}^{(0)}}\right\|_{2}^{2}+\frac{1}{\lambda_{2}}D_{H}(\boldsymbol{s}^{*},\boldsymbol{s}^{(0)})\right).$
(9)

This convergence rate is optimal: it is the best achievable rate of
convergence in the Nesterov class of optimal first-order methods [51].

An important feature of our proposed algorithm is that the minimization steps
in (LABEL:eq:kl-nPDHG-alg) can be computed exactly. Specifically, with the
auxiliary variables $\boldsymbol{u}^{(k)}=\bm{A}\boldsymbol{\theta}^{(k)}$ and
$v_{i}^{(k)}=\log\left(s_{i}^{(k)}/(1-s_{i}^{(k)})\right)$ for $i\in\\{1,\dots
m\\}$, the steps in algorithm (LABEL:eq:kl-nPDHG-alg) can be expressed
explicitly as follows:

$\begin{dcases}\boldsymbol{v}^{(k+1)}&=\frac{1}{1+\sigma}\left(\sigma\boldsymbol{u}^{(k)}+\sigma\rho\left(\boldsymbol{u}^{(k)}-\boldsymbol{u}^{(k-1)}\right)+\boldsymbol{v}^{(k)}\right)\\\
s^{(k+1)}_{i}&=\frac{1}{1+\exp{\left(-v_{i}^{(k+1)}\right)}}\quad\mathrm{for}\,i\in\left\\{1,\dots,m\right\\},\\\
\hat{\boldsymbol{\theta}}^{(k+1)}&=\boldsymbol{\theta}^{(k)}-\tau\bm{A}^{T}\left(\boldsymbol{s}^{(k+1)}-\boldsymbol{y}\right)\\\
\theta^{(k+1)}_{j}&=\mathrm{sign~{}}{\hat{\theta}^{(k+1)}_{j}}\max{\left(0,\frac{\left|\hat{\theta}^{(k+1)}_{j}\right|-\lambda_{1}\tau}{1+\lambda_{2}\tau}\right)}\quad\mathrm{for}\,j\in\\{1,\dots,n\\}\\\
\boldsymbol{u}^{(k+1)}&=\bm{A}\boldsymbol{\theta}^{(k+1)}.\end{dcases}$ (10)

In addition, from the auxiliary variables and the optimality condition on the
right in (5), we have the limit

$\lim_{k\to+\infty}\left\|{\boldsymbol{u}^{(k)}-\boldsymbol{v}^{(k)}}\right\|_{2}=0,$

which can serve as a convergence criterion. We refer to Material and Methods
for the derivation of algorithm (10) from the iterations in (LABEL:eq:kl-
nPDHG-alg).

Our proposed explicit nonlinear PDHG algorithm (10) offers many advantages in
terms of efficiency and robustness. First, the computational bottlenecks in
algorithm (10) consist of matrix-vector multiplications and the estimation of
the induced matrix $\left\|{\bm{A}}\right\|_{op}$ given by (8). If $T(m,n)$
denotes the number of arithmetic operations required to perform matrix-vector
multiplication with the matrix $\bm{A}$, then the asymptotic space complexity
for computing the iterations in algorithm (10) as well as the induced matrix
$\left\|{\bm{A}}\right\|_{op}$ is at most $O(T(m,n))$ operations. As mentioned
before, this is unlike most first-order optimization methods, such as the
forward-backward splitting algorithm, where instead the matrix norm is the
largest singular value of the matrix $\bm{A}$, which takes
$O(\min{(m^{2}n,mn^{2})})$ operations to compute. This fact is crucial because
the smaller computational cost makes it easy and efficient to estimate all the
parameters of the nonlinear PDHG algorithm, which is needed to achieve an
optimal rate of convergence.

Another advantage of our algorithm is that it exhibits scalable parallelism
because the matrix-vector multiplication operations can be implemented via
parallel algorithms. This makes it possible to implement our proposed
algorithms in a way that takes advantage of emerging hardware, such as field-
programmable gate arrays architectures.

Finally, our algorithm also provably computes an $\epsilon$-approximate
solution of (1) in $O(\log(1/\epsilon))$ operations [16, Section 5.4]. The
size of the parameter $\rho$ dictates this linear rate of convergence; it
depends on the matrix $\bm{A}$, the tuning parameter $\lambda$, and the
hyperparameter $\alpha$. Hence the overall complexity required to compute a
global minimum of the elastic net regularized logistic regression problem (1)
with tolerance $\epsilon\in(0,1)$ is on the order of
$O(T(m,n)\log(1/\epsilon))$ operations.

With these advantages, algorithm (10) overcomes the limitations of the state-
of-the-art and other widely-used algorithms for solving the logistic
regression problem (1). We are unaware of any other algorithm that offers
these advantages in terms of efficiency and robustness simultaneously.

In general, the nonlinear PDHG algorithm (10) can be adapted to any
regularized logistic regression problem for which the penalty is strongly
convex on the space $(\mathbb{R}^{n},\left\|{\cdot}\right\|_{2})$. To do so,
substitute this penalty for the elastic net penalty in the minimization
problem of the scheme (LABEL:eq:kl-nPDHG-alg) and use its solution in place of
the third and fourth lines in the explicit algorithm (10).

### Special case: Logistic regression with the lasso penalty

In some situations, it may be desirable to fit the regularized logistic
regression model (1) without the $\ell{2}^{2}$ penalty ($\alpha=1$). In this
case, algorithm (10) does not apply since it depends on the strong convexity
of $\ell{2}^{2}$ penalty. We present here an algorithm for fitting a logistic
regression model regularized by an $\ell 1$ penalty, or in principle, any
convex penalty that is not strongly convex, such as the group lasso.

The $\ell{1}$-regularized logistic regression problem is

$\inf_{\boldsymbol{\theta}\in\mathbb{R}^{n}}\left\\{\frac{1}{m}\sum_{i=1}^{m}\log\left(1+\exp{((\bm{A}\boldsymbol{\theta})_{i})}\right)-\frac{1}{m}\left\langle\boldsymbol{y},\bm{A}\boldsymbol{\theta}\right\rangle+\lambda\left\|{\boldsymbol{\theta}}\right\|_{1}\right\\},$
(11)

and its associated saddle-point problem is

$\inf_{\boldsymbol{\theta}\in\mathbb{R}^{n}}\sup_{\boldsymbol{s}\in[0,1]^{m}}\left\\{-\frac{1}{m}\psi(\boldsymbol{s})-\frac{1}{m}\left\langle\boldsymbol{y}-\boldsymbol{s},\bm{A}\boldsymbol{\theta}\right\rangle+\lambda\left\|{\boldsymbol{\theta}}\right\|_{1}\right\\}.$
(12)

The $\ell 1$ penalty in (11) guarantees that problem (11) has at least one
solution. Accordingly, the saddle-point problem (12) also has at least one
saddle point. As before, we split the infimum and supremum problems in (12) by
alternating between a nonlinear proximal ascent step using the
Kullback–Leibler divergence (2) and a proximal descent step using a quadratic
function, but this time we also update the stepsize parameters at each
iteration:

$\displaystyle\boldsymbol{s}^{(k+1)}=\operatorname*{arg\,max}_{\boldsymbol{s}\in(0,1)^{m}}\left\\{-\psi(\boldsymbol{s})+\left\langle\boldsymbol{s},\bm{A}(\boldsymbol{\theta}^{(k)}+\rho^{(k)}(\boldsymbol{\theta}^{(k)}-\boldsymbol{\theta}^{(k-1)}))\right\rangle-\frac{1}{\sigma^{(k)}}D_{H}(\boldsymbol{s},\boldsymbol{s}^{(k)})\right\\},$
(13)
$\displaystyle\boldsymbol{\theta}^{(k+1)}=\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\mathbb{R}^{n}}\left\\{m\lambda\left\|{\boldsymbol{\theta}}\right\|_{1}+\left\langle\boldsymbol{s}^{(k+1)}-\boldsymbol{y},\bm{A}\boldsymbol{\theta}\right\rangle+\frac{1}{2\tau^{(k)}}\left\|{\boldsymbol{\theta}-\boldsymbol{\theta}^{(k)}}\right\|_{2}^{2}\right\\},$
$\displaystyle\rho^{(k+1)}=1/\sqrt{1+\sigma^{(k)}},\quad\sigma^{(k+1)}=\rho^{(k+1)}\sigma^{(k)},\quad\tau^{(k+1)}=\tau^{(k)}/\rho^{(k+1)}.$

The scheme starts from initial stepsize parameters $\rho^{(0)}\in(0,1)$,
$\tau^{(0)}>0$ and
$\sigma^{(0)}=1/(\tau^{(0)}\left\|{\bm{A}}\right\|_{op}^{2})$, and initial
values $\boldsymbol{s}^{(0)}\in(0,1)^{m}$ and
$\boldsymbol{\theta}^{(-1)}=\boldsymbol{\theta}^{(0)}\in\mathbb{R}^{n}$.

The following accelerated nonlinear PDHG algorithm computes a global minimum
of (11):

$\begin{dcases}\boldsymbol{v}^{(k+1)}&=\frac{1}{1+\sigma^{(k)}}\left(\sigma^{(k)}\boldsymbol{u}^{(k)}+\sigma^{(k)}\rho^{(k)}\left(\boldsymbol{u}^{(k)}-\boldsymbol{u}^{(k-1)}\right)+\boldsymbol{v}^{(k)}\right)\\\
s^{(k+1)}_{i}&=\frac{1}{1+\exp{\left(-v_{i}^{(k+1)}\right)}}\quad\mathrm{for}\,i\in\left\\{1,\dots,m\right\\},\\\
\hat{\boldsymbol{\theta}}^{(k+1)}&=\boldsymbol{\theta}^{(k)}-\tau^{(k)}\bm{A}^{T}\left(\boldsymbol{s}^{(k+1)}-\boldsymbol{y}\right)\\\
\theta^{(k+1)}_{j}&=\mathrm{sign~{}}{\hat{\theta}^{(k+1)}_{j}}\max{\left(0,\left|\hat{\theta}^{(k+1)}_{j}\right|-m\lambda\tau^{(k)}\right)}\quad\mathrm{for}\,j\in\\{1,\dots,n\\},\\\
\boldsymbol{u}^{(k+1)}&=\bm{A}\boldsymbol{\theta}^{(k+1)},\\\
\rho^{(k+1)}&=1/\sqrt{1+\sigma^{(k)}},\quad\sigma^{(k+1)}=\rho^{(k+1)}\sigma^{(k)},\quad\tau^{(k+1)}=\tau^{(k)}/\rho^{(k+1)}.\end{dcases}$
(14)

In addition, from the auxiliary variables and the optimality condition on the
right in (5), we have the limit

$\lim_{k\to+\infty}\left\|{\boldsymbol{u}^{(k)}-\boldsymbol{v}^{(k)}}\right\|_{2}=0,$

which can serve as a convergence criterion. The derivation of algorithm (14)
from the iterations in (LABEL:eq:kl-lasso-alg) follows from the derivation of
algorithm (10) from the iterations in (LABEL:eq:kl-nPDHG-alg) described in
Material and Methods by setting $\alpha=1$. According to results provided by
the authors in [16, Proposition 5.2], the sequence of iterates
$\\{(\boldsymbol{\theta}^{(k)},\boldsymbol{s}^{(k)})\\}_{k=1}^{+\infty}$
converges to a saddle point $(\boldsymbol{\theta}^{*},\boldsymbol{s}^{*})$ of
(12) at a sublinear rate $O(1/k^{2})$ in the iterations. Moreover, this
sublinear rate satisfies the lower bound

$\frac{2\tau^{(0)}\left\|{\bm{A}}\right\|^{2}_{op}}{1+2\tau^{(0)}\left\|{\bm{A}}\right\|^{2}_{op}}k+\frac{2\tau^{(0)}}{(1+2\tau^{(0)}\left\|{\bm{A}}\right\|^{2}_{op})^{2}}k^{2}.$

In particular, the constant term multiplying $k^{2}$ is maximized when
$\tau^{(0)}=1/(2\left\|{\bm{A}}\right\|_{op}^{2})$. This suggests a practical
choice for the free parameter $\tau^{(0)}$.

In general, the nonlinear PDHG algorithm (14) can be adapted to any
regularized logistic regression problem for which the penalty is proper, lower
semicontinuous and convex, and for which a solution exists. To do so,
substitute the $\ell 1$ penalty in the minimization problem of the scheme
(LABEL:eq:kl-nPDHG-alg) and use its solution in place of the third and fourth
lines in the explicit algorithm (14).

##  4 Material and Methods

### Derivation of the explicit algorithm (10)

We derive here the explicit algorithm (LABEL:eq:kl-nPDHG-alg) from the
iterations in (10). Consider the first line of (LABEL:eq:kl-nPDHG-alg). This
maximization problem has a unique maximum inside the interval $(0,1)^{m}$ [3,
Proposition 3.21-3.23, Theorem 3.24, Corollary 3.25], and the objective
function is differentiable. Thus it suffices to compute the gradient with
respect to $\boldsymbol{s}$ and solve for $\boldsymbol{s}$ to compute its
global maximum. To do so, it helps to first rearrange the objective function.
Substitute $\boldsymbol{s}^{(k)}$ for $\boldsymbol{s}^{\prime}$ in equation
(2), use equation (7), and rearrange to obtain the objective function

$\displaystyle-\psi(\boldsymbol{s})$
$\displaystyle+\left\langle\boldsymbol{s},\bm{A}(\boldsymbol{\theta}^{(k)}+\rho(\boldsymbol{\theta}^{(k)}-\boldsymbol{\theta}^{(k-1)}))\right\rangle-\frac{1}{\sigma}D_{H}(\boldsymbol{s},\boldsymbol{s}^{(k)})$
$\displaystyle=-\psi(\boldsymbol{s})+\left\langle\boldsymbol{s},\bm{A}(\boldsymbol{\theta}^{(k)}+\rho(\boldsymbol{\theta}^{(k)}-\boldsymbol{\theta}^{(k-1)}))\right\rangle-\frac{1}{\sigma}\left(\psi(\boldsymbol{s})-\psi(\boldsymbol{s}^{(k)})-\left\langle\boldsymbol{s}-\boldsymbol{s}^{(k)},\nabla\psi(\boldsymbol{s}^{(k)})\right\rangle)\right)$
$\displaystyle=-\left(1+\frac{1}{\sigma}\right)\psi(\boldsymbol{s})+\left\langle\boldsymbol{s},\bm{A}(\boldsymbol{\theta}^{(k)}+\rho(\boldsymbol{\theta}^{(k)}-\boldsymbol{\theta}^{(k-1)}))\right\rangle+\frac{1}{\sigma}\left\langle\boldsymbol{s}-\boldsymbol{s}^{(k)},\nabla\psi(\boldsymbol{s}^{(k)})\right\rangle+\psi(\boldsymbol{s}^{(k)}).$

The optimality condition is then

$\nabla\psi(\boldsymbol{s}^{(k+1)})=\frac{\sigma}{1+\sigma}\left(\bm{A}(\boldsymbol{\theta}^{(k)}+\rho(\boldsymbol{\theta}^{(k)}-\boldsymbol{\theta}^{(k-1)}))\right)+\frac{1}{1+\sigma}\nabla\psi(\boldsymbol{s}^{(k)}),$

where
$(\nabla\psi(\boldsymbol{s}))_{i}=\log\left(s_{i}^{(k+1)}/(1-s_{i}^{(k+1)})\right)$
for $i\in\\{1,\dots,m\\}$ and $\boldsymbol{s}\in(0,1)^{m}$. With the auxiliary
variables $\boldsymbol{u}^{(k)}=\bm{A}\boldsymbol{\theta}^{(k)}$ and
$v_{i}^{(k)}=\log\left(s_{i}^{(k)}/(1-s_{i}^{(k)})\right)$ for $i\in\\{1,\dots
m\\}$, the optimality condition can be written as

$\boldsymbol{v}^{(k+1)}=\frac{1}{1+\sigma^{(k)}}\left(\sigma^{(k)}\boldsymbol{u}^{(k)}+\sigma^{(k)}\rho^{(k)}\left(\boldsymbol{u}^{(k)}-\boldsymbol{u}^{(k-1)}\right)+\boldsymbol{v}^{(k)}\right).$

This gives the first line in (10). The second line follows upon solving for
$\boldsymbol{s}^{(k+1)}$ in terms of $\boldsymbol{v}^{(k+1)}$. The fifth line
follows from the definition of the auxiliary variable $\boldsymbol{u}^{(k)}$.

Now, consider the second line of (LABEL:eq:kl-nPDHG-alg). Complete the square
and multiply by $\tau/(1+\lambda_{2}\tau)$ to get the equivalent minimization
problem

$\boldsymbol{\theta}^{(k+1)}=\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\mathbb{R}^{n}}\left\\{\frac{\lambda_{1}\tau}{1+\lambda_{2}\tau}\left\|{\boldsymbol{\theta}}\right\|_{1}+\frac{1}{2}\left\|{\boldsymbol{\theta}-\left(\boldsymbol{\theta}^{(k)}-\tau\bm{A}^{T}(\boldsymbol{s}^{(k+1)}-\boldsymbol{y})\right)/(1+\lambda_{2}\tau)}\right\|_{2}^{2}\right\\}.$

The unique minimum is computed using the soft thresholding operator [17, 25,
44]. With the notation

$\hat{\boldsymbol{\theta}}^{(k+1)}=\boldsymbol{\theta}^{(k)}-\tau\bm{A}^{T}\left(\boldsymbol{s}^{(k+1)}-\boldsymbol{y}\right),$

the soft thresholding operator is defined component-wise by

$\theta^{(k+1)}_{i}=\mathrm{sign~{}}{\hat{\theta}^{(k+1)}_{j}}\max{\left(0,\frac{\left|\hat{\theta}^{(k+1)}_{j}\right|-\lambda_{1}\tau}{1+\lambda_{2}\tau}\right)}\quad\mathrm{for}\,j\in\\{1,\dots,n\\}.$

The third and fourth lines of (10) are precisely these two equations.

## References

  * Achia et al. [2010] Thomas NO Achia, Anne Wangombe, and Nancy Khadioli.  A logistic regression model to identify key determinants of poverty using demographic and health survey data.  _European Journal of Social Sciences_ , 13(1), 2010.
  * Bagley et al. [2001] Steven C Bagley, Halbert White, and Beatrice A Golomb.  Logistic regression in the medical literature:: Standards for use and reporting, with particular attention to one medical domain.  _Journal of Clinical Epidemiology_ , 54(10):979–985, 2001.  ISSN 0895-4356.  doi: https://doi.org/10.1016/S0895-4356(01)00372-9.
  * Bauschke et al. [2003] Heinz H Bauschke, Jonathan M Borwein, and Patrick L Combettes.  Bregman monotone optimization algorithms.  _SIAM Journal on control and optimization_ , 42(2):596–636, 2003.
  * Beck and Teboulle [2003] Amir Beck and Marc Teboulle.  Mirror descent and nonlinear projected subgradient methods for convex optimization.  _Operations Research Letters_ , 31(3):167–175, 2003.
  * Beck and Teboulle [2009] Amir Beck and Marc Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  _SIAM journal on imaging sciences_ , 2(1):183–202, 2009.
  * Berger et al. [1996] Adam Berger, Stephen A Della Pietra, and Vincent J Della Pietra.  A maximum entropy approach to natural language processing.  _Computational linguistics_ , 22(1):39–71, 1996\.
  * Bertsimas et al. [2019] Dimitris Bertsimas, Jean Pauphilet, and Bart Van Parys.  Sparse regression: Scalable algorithms and empirical performance.  _arXiv preprint arXiv:1902.06547_ , 2019.
  * Bringmann et al. [2018] Björn Bringmann, Daniel Cremers, Felix Krahmer, and Michael Möller.  The homotopy method revisited: Computing solution paths of $\ell_{1}$-regularized problems.  _Mathematics of Computation_ , 87(313):2343–2364, 2018.  URL https://doi.org/10.1090/mcom/3287.
  * Bursac et al. [2008] Zoran Bursac, C Heath Gauss, David Keith Williams, and David W Hosmer.  Purposeful selection of variables in logistic regression.  _Source code for biology and medicine_ , 3(1):1–8, 2008.
  * Catalina et al. [2018] Alejandro Catalina, Carlos M Alaíz, and José R Dorronsoro.  scho.  In _2018 International Joint Conference on Neural Networks (IJCNN)_ , pages 1–8. IEEE, 2018.
  * Chambolle and Pock [2016a] A. Chambolle and T. Pock.  An introduction to continuous optimization for imaging.  _Acta Numer._ , 25:161–319, 2016a.
  * Chambolle and Pock [2011] Antonin Chambolle and Thomas Pock.  A first-order primal-dual algorithm for convex problems with applications to imaging.  _Journal of mathematical imaging and vision_ , 40(1):120–145, 2011.
  * Chambolle and Pock [2016b] Antonin Chambolle and Thomas Pock.  On the ergodic convergence rates of a first-order primal–dual algorithm.  _Mathematical Programming_ , 159(1-2):253–287, 2016b.
  * Chu et al. [2007] Cheng Chu, Sang Kyun Kim, Yian Lin, YuanYuan Yu, Gary Bradski, Andrew Y Ng, and Kunle Olukotun.  Map-reduce for machine learning on multicore.  _Advances in neural information processing systems_ , 19:281, 2007.
  * Csiszár [1967] Imre Csiszár.  Information-type measures of difference of probability distributions and indirect observation.  _studia scientiarum Mathematicarum Hungarica_ , 2:229–318, 1967.
  * Darbon and Langlois [2021] Jérôme Darbon and Gabriel Provencher Langlois.  Accelerated nonlinear primal-dual hybrid gradient algorithms with applications to machine learning, 2021.
  * Daubechies et al. [2004] Ingrid Daubechies, Michel Defrise, and Christine De Mol.  An iterative thresholding algorithm for linear inverse problems with a sparsity constraint.  _Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences_ , 57(11):1413–1457, 2004.
  * Demchenko et al. [2013] Yuri Demchenko, Paola Grosso, Cees De Laat, and Peter Membrey.  Addressing big data issues in scientific data infrastructure.  In _2013 International Conference on Collaboration Technologies and Systems (CTS)_ , pages 48–55. IEEE, 2013.
  * Dhar [2020] Payal Dhar.  The carbon impact of artificial intelligence.  _Nat Mach Intell_ , 2:423–5, 2020.
  * Efron et al. [2004] Bradley Efron, Trevor Hastie, Iain Johnstone, Robert Tibshirani, et al.  Least angle regression.  _The Annals of statistics_ , 32(2):407–499, 2004\.
  * Ekeland and Temam [1999] Ivar Ekeland and Roger Temam.  _Convex analysis and variational problems_.  SIAM, 1999.
  * El Guide et al. [2020] M El Guide, K Jbilou, C Koukouvinos, and A Lappa.  Comparative study of l 1 regularized logistic regression methods for variable selection.  _Communications in Statistics-Simulation and Computation_ , pages 1–16, 2020.
  * Esser et al. [2010] Ernie Esser, Xiaoqun Zhang, and Tony F Chan.  A general framework for a class of first order primal-dual algorithms for convex optimization in imaging science.  _SIAM Journal on Imaging Sciences_ , 3(4):1015–1046, 2010.
  * Fercoq and Richtárik [2016] Olivier Fercoq and Peter Richtárik.  Optimization in high dimensions via accelerated, parallel, and proximal coordinate descent.  _SIAM Review_ , 58(4):739–771, 2016.
  * Figueiredo and Nowak [2001] Mário AT Figueiredo and Robert D Nowak.  Wavelet-based image estimation: an empirical Bayes approach using Jeffrey’s noninformative prior.  _IEEE Transactions on Image Processing_ , 10(9):1322–1331, 2001.
  * Foucart and Rauhut [2013] Simon Foucart and Holger Rauhut.  _Sparse Solutions of Underdetermined Systems_ , pages 41–59.  Springer New York, New York, NY, 2013.  doi: 10.1007/978-0-8176-4948-7˙2.
  * Friedman et al. [2007] Jerome Friedman, Trevor Hastie, Holger Höfling, Robert Tibshirani, et al.  Pathwise coordinate optimization.  _The annals of applied statistics_ , 1(2):302–332, 2007.
  * Friedman et al. [2010] Jerome Friedman, Trevor Hastie, and Rob Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  _Journal of statistical software_ , 33(1):1, 2010\.
  * Genkin et al. [2007] Alexander Genkin, David D Lewis, and David Madigan.  Large-scale bayesian logistic regression for text categorization.  _Technometrics_ , 49(3):291–304, 2007.
  * Greene et al. [2014] Casey S Greene, Jie Tan, Matthew Ung, Jason H Moore, and Chao Cheng.  Big data bioinformatics.  _Journal of cellular physiology_ , 229(12):1896–1900, 2014.
  * Hastie et al. [2009] Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  _The elements of statistical learning_.  Springer Series in Statistics. Springer, New York, second edition, 2009\.  doi: 10.1007/978-0-387-84858-7.  Data mining, inference, and prediction.
  * Hastie et al. [2021] Trevor Hastie, Junyang Qian, and Kenneth Tay.  An introduction to glmnet (2021).  Available at ”https://glmnet.stanford.edu/articles/glmnet.html”. Accessed on 17 February 2021., 2021.
  * Hesterberg et al. [2008] Tim Hesterberg, Nam Hee Choi, Lukas Meier, Chris Fraley, et al.  Least angle and $l^{1}$ penalized regression: A review.  _Statistics Surveys_ , 2:61–93, 2008.
  * Hoerl and Kennard [1970] Arthur E Hoerl and Robert W Kennard.  Ridge regression: Biased estimation for nonorthogonal problems.  _Technometrics_ , 12(1):55–67, 1970.
  * Hohage and Homann [2014] Thorsten Hohage and Carolin Homann.  A generalization of the chambolle-pock algorithm to banach spaces with applications to inverse problems.  _arXiv preprint arXiv:1412.0126_ , 2014.
  * Hosmer Jr et al. [2013] David W Hosmer Jr, Stanley Lemeshow, and Rodney X Sturdivant.  _Applied logistic regression_ , volume 398.  John Wiley & Sons, 2013.
  * Kambatla et al. [2014] Karthik Kambatla, Giorgos Kollias, Vipin Kumar, and Ananth Grama.  Trends in big data analytics.  _Journal of parallel and distributed computing_ , 74(7):2561–2573, 2014.
  * Kemperman [1969] Johannes HB Kemperman.  On the optimum rate of transmitting information.  In _Probability and information theory_ , pages 126–169. Springer, 1969.
  * King and Zeng [2001] Gary King and Langche Zeng.  Logistic regression in rare events data.  _Political analysis_ , 9(2):137–163, 2001.
  * Kullback [1967] Solomon Kullback.  A lower bound for discrimination information in terms of variation (corresp.).  _IEEE transactions on Information Theory_ , 13(1):126–127, 1967.
  * Lee et al. [2006] Su-In Lee, Honglak Lee, Pieter Abbeel, and Andrew Y Ng.  Efficient l~ 1 regularized logistic regression.  In _Aaai_ , volume 6, pages 401–408, 2006.
  * Leiserson et al. [2020] Charles E Leiserson, Neil C Thompson, Joel S Emer, Bradley C Kuszmaul, Butler W Lampson, Daniel Sanchez, and Tao B Schardl.  There’s plenty of room at the top: What will drive computer performance after moore’s law?  _Science_ , 368(6495), 2020.
  * Li et al. [2020] Xiaoping Li, Yadi Wang, and Rubén Ruiz.  A survey on sparse learning models for feature selection.  _IEEE Transactions on Cybernetics_ , 2020.
  * Lions and Mercier [1979] Pierre-Louis Lions and Bertrand Mercier.  Splitting algorithms for the sum of two nonlinear operators.  _SIAM Journal on Numerical Analysis_ , 16(6):964–979, 1979.
  * Loris [2008] Ignace Loris.  L1packv2: A mathematica package for minimizing an l1-penalized functional.  _Computer physics communications_ , 179(12):895–902, 2008.
  * Lowe and Parvar [2004] David J Lowe and Jamshid Parvar.  A logistic regression approach to modelling the contractor’s decision to bid.  _Construction Management and Economics_ , 22(6):643–653, 2004.
  * L’heureux et al. [2017] Alexandra L’heureux, Katarina Grolinger, Hany F Elyamany, and Miriam AM Capretz.  Machine learning with big data: Challenges and approaches.  _Ieee Access_ , 5:7776–7797, 2017.
  * Manning and Klein [2003] Christopher Manning and Dan Klein.  Optimization, maxent models, and conditional estimation without magic.  In _Proceedings of the 2003 Conference of the North American Chapter of the Association for Computational Linguistics on Human Language Technology: Tutorials-Volume 5_ , pages 8–8, 2003.
  * Meier et al. [2008] Lukas Meier, Sara Van De Geer, and Peter Bühlmann.  The group lasso for logistic regression.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 70(1):53–71, 2008.
  * Muchlinski et al. [2016] David Muchlinski, David Siroky, Jingrui He, and Matthew Kocher.  Comparing random forest with logistic regression for predicting class-imbalanced civil war onset data.  _Political Analysis_ , pages 87–103, 2016.
  * Nesterov [2018] Yurii Nesterov.  _Lectures on Convex Optimization_.  Springer International Publishing, 2018.
  * Pereira et al. [2009] Francisco Pereira, Tom Mitchell, and Matthew Botvinick.  Machine learning classifiers and fmri: A tutorial overview.  _NeuroImage_ , 45(1, Supplement 1):S199–S209, 2009.  ISSN 1053-8119.  Mathematics in Brain Imaging.
  * Pinsker [1964] Mark S Pinsker.  _Information and information stability of random variables and processes_.  Holden-Day, 1964.
  * Pock et al. [2009] Thomas Pock, Daniel Cremers, Horst Bischof, and Antonin Chambolle.  An algorithm for minimizing the mumford-shah functional.  In _2009 IEEE 12th International Conference on Computer Vision_ , pages 1133–1140. IEEE, 2009.
  * Pranckevičius and Marcinkevičius [2017] Tomas Pranckevičius and Virginijus Marcinkevičius.  Comparison of naive bayes, random forest, decision tree, support vector machines, and logistic regression classifiers for text reviews classification.  _Baltic Journal of Modern Computing_ , 5(2):221, 2017.
  * Privé et al. [2018] Florian Privé, Hugues Aschard, Andrey Ziyatdinov, and Michael GB Blum.  Efficient analysis of large-scale genome-wide data with two r packages: bigstatsr and bigsnpr.  _Bioinformatics_ , 34(16):2781–2787, 2018.
  * Sculley et al. [2014] D. Sculley, Gary Holt, Daniel Golovin, Eugene Davydov, Todd Phillips, Dietmar Ebner, Vinay Chaudhary, and Michael Young.  Machine learning: The high interest credit card of technical debt.  In _SE4ML: Software Engineering for Machine Learning (NIPS 2014 Workshop)_ , 2014.
  * Shi et al. [2010] Jianing Shi, Wotao Yin, Stanley Osher, and Paul Sajda.  A fast hybrid algorithm for large-scale l1-regularized logistic regression.  _The Journal of Machine Learning Research_ , 11:713–741, 2010.
  * Shi et al. [2013] Jianing V Shi, Wotao Yin, and Stanley J Osher.  Linearized bregman for l1-regularized logistic regression.  In _Proceedings of the 30th international conference on machine learning (ICML)_. Citeseer, 2013.
  * Simon et al. [2011] Noah Simon, Jerome Friedman, Trevor Hastie, and Rob Tibshirani.  Regularization paths for cox’s proportional hazards model via coordinate descent.  _Journal of statistical software_ , 39(5):1, 2011\.
  * Simon et al. [2013] Noah Simon, Jerome Friedman, and Trevor Hastie.  A blockwise descent algorithm for group-penalized multiresponse and multinomial regression.  _arXiv preprint arXiv:1311.6529_ , 2013.
  * Taddy [2013] Matt Taddy.  Multinomial inverse regression for text analysis.  _Journal of the American Statistical Association_ , 108(503):755–770, 2013.
  * Tay et al. [2021] J Kenneth Tay, Balasubramanian Narasimhan, and Trevor Hastie.  Elastic net regularization paths for all generalized linear models.  _arXiv preprint arXiv:2103.03475_ , 2021.
  * Theodossiou [1998] Ioannis Theodossiou.  The effects of low-pay and unemployment on psychological well-being: a logistic regression approach.  _Journal of health economics_ , 17(1):85–104, 1998.
  * Thompson et al. [2021] Neil C Thompson, Kristjan Greenewald, Keeheon Lee, and Gabriel F Manso.  Deep learning’s diminishing returns: The cost of improvement is becoming unsustainable.  _IEEE Spectrum_ , 58(10):50–55, 2021.
  * Tibshirani [1996] Robert Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 58(1):267–288, 1996.
  * Tibshirani et al. [2012] Robert Tibshirani, Jacob Bien, Jerome Friedman, Trevor Hastie, Noah Simon, Jonathan Taylor, and Ryan J Tibshirani.  Strong rules for discarding predictors in lasso-type problems.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 74(2):245–266, 2012.
  * Tibshirani et al. [2013] Ryan J Tibshirani et al.  The lasso problem and uniqueness.  _Electronic Journal of statistics_ , 7:1456–1490, 2013\.
  * Vidaurre et al. [2013] Diego Vidaurre, Concha Bielza, and Pedro Larrañaga.  A survey of l1 regression.  _International Statistical Review_ , 81(3):361–387, 2013.
  * Wu et al. [2008] Tong Tong Wu, Kenneth Lange, et al.  Coordinate descent algorithms for lasso penalized regression.  _The Annals of Applied Statistics_ , 2(1):224–244, 2008.
  * Wu et al. [2009] Tong Tong Wu, Yi Fang Chen, Trevor Hastie, Eric Sobel, and Kenneth Lange.  Genome-wide association analysis by lasso penalized logistic regression.  _Bioinformatics_ , 25(6):714–721, 2009.
  * Yuan et al. [2010] Guo-Xun Yuan, Kai-Wei Chang, Cho-Jui Hsieh, and Chih-Jen Lin.  A comparison of optimization methods and software for large-scale l1-regularized linear classification.  _The Journal of Machine Learning Research_ , 11:3183–3234, 2010.
  * Yuan et al. [2012] Guo-Xun Yuan, Chia-Hua Ho, and Chih-Jen Lin.  An improved glmnet for l1-regularized logistic regression.  _The Journal of Machine Learning Research_ , 13:1999–2030, 2012.
  * Zaghdoudi [2013] Taha Zaghdoudi.  Bank failure prediction with logistic regression.  _International Journal of Economics and Financial Issues_ , 3(2):537, 2013.
  * Zaidi and Amirat [2016] M Zaidi and A Amirat.  Forecasting stock market trends by logistic regression and neural networks: Evidence from ksa stock market.  _Int. J. Econ. Commer. Manag_ , 4:4–7, 2016.
  * Zanon et al. [2020] Mattia Zanon, Giuliano Zambonin, Gian Antonio Susto, and Seán McLoone.  Sparse logistic regression: Comparison of regularization and bayesian implementations.  _Algorithms_ , 13(6):137, 2020.
  * Zhang et al. [2018] Zhongheng Zhang, Victor Trevino, Sayed Shahabuddin Hoseini, Smaranda Belciug, Arumugam Manivanna Boopathi, Ping Zhang, Florin Gorunescu, Velappan Subha, and Songshi Dai.  Variable selection in logistic regression model with genetic algorithm.  _Annals of translational medicine_ , 6(3), 2018.
  * Zhu and Chan [2008] Mingqiang Zhu and Tony Chan.  An efficient primal-dual hybrid gradient algorithm for total variation image restoration.  _UCLA CAM Report_ , 34, 2008.
  * Zou and Hastie [2005] Hui Zou and Trevor Hastie.  Regularization and variable selection via the elastic net.  _Journal of the royal statistical society: series B (statistical methodology)_ , 67(2):301–320, 2005.


# Phase Space Analysis of Cardiac Spectra

Onder Pekcan
Department of Molecular Biology and Genetics
Kadir Has University
Istanbul, Turkey
<EMAIL_ADDRESS>
Taner Arsan
Department of Computer Engineering
Kadir Has University
Istanbul, Turkey
<EMAIL_ADDRESS>
Corresponding Author

###### Abstract

Cardiac diseases are one of the main reasons of mortality in modern,
industrialized societies, and they cause high expenses in public health
systems. Therefore, it is important to develop analytical methods to improve
cardiac diagnostics. Electric activity of heart was first modeled by using a
set of nonlinear differential equations. Latter, variations of cardiac spectra
originated from deterministic dynamics are investigated. Analyzing the power
spectra of a normal human heart presents His-Purkinje network, possessing a
fractal like structure. Phase space trajectories are extracted from the time
series graph of ECG. Lower values of fractal dimension, $D$ indicate dynamics
that are more coherent. If $D$ has non-integer values greater than two when
the system becomes chaotic or strange attractor. Recently, the development of
a fast and robust method, which can be applied to multichannel physiologic
signals, was reported. This manuscript investigates two different ECG systems
produced from normal and abnormal human hearts to introduce an auxiliary phase
space method in conjunction with ECG signals for diagnoses of heart diseases.
Here, the data for each person includes two signals based on $V_{4}$ and
modified lead III (MLIII) respectively. Fractal analysis method is employed on
the trajectories constructed in phase space, from which the fractal dimension
$D$ is obtained using the box counting method. It is observed that, MLIII
signals have larger $D$ values than the first signals ($V_{4}$), predicting
more randomness yet more information. The lowest value of $D$ (1.708)
indicates the perfect oscillation of the normal heart and the highest value of
$D$ ( 1.863) presents the randomness of the abnormal heart. Our significant
finding is that the phase space picture presents the distribution of the peak
heights from the ECG spectra, giving valuable information about heart
activities in conjunction with ECG.

_Keywords_ Electrocardiography $\cdot$ Analysis $\cdot$ Diagnostic method
$\cdot$ Cardiac electrophysiology $\cdot$ Computer-based model

##  1 Introduction

It is well known that human heart has an electric activity, which can be
detected by measuring the potential difference from various points on the
surface of the body. The measured electric potential versus time is called
electrocardiogram (ECG), which possesses three separate parts. $P$-wave
presents the excitation of the atria, $QRS$ complex shows the ventricles (His-
Purkinje network) and $T$-wave is associated with the recovery of initial
electrical state of ventricles (see Figure 1(a)). Although ECG presents
periodic behavior, some irregularities can be seen in details of the record.
In fact, these irregularities belong to the intrinsic part of heart activity
and/or to the random noise that can be found in such systems. These activities
in ECG spectra are highly important to understand the cardiac dynamics.
Cardiac oscillations are sometimes perturbed by unpredictable contributions
which are part of the cardiac dynamics and therefore physiologically
important. These finding predict that heart is not a perfect oscillator and/or
cardiac muscles do not always vibrate harmonically.

It has been also well established that the transformation of a sequence of
values in time to a geometrical object in space is a highly considered topic
by Liebovitch (1998). This procedure replaces an analysis in time with an
analysis in space. Here, the space is called phase space and the procedure of
transforming the time series into the object in space is called an embedding.
Afterwards, topological properties of the object are determined based on its
fractal dimension. The fractal dimension characterizes the properties of the
phase space set, not the original time series. The measured dimension in the
phase space set determines whether a data set is generated by random or
deterministic processes. A large value of the fractal dimension indicates
random generation of the time series. This means that, the number of variables
and equations are so large that there are no ways to predict the future values
from the early parameters. In other words, the multiplicity of interacting
factors precludes the possibility of understanding how the underlying
mechanisms work. On the other hand, a low fractal dimension value indicates
that the data is generated by deterministic mechanisms, based on a small
number of independent variables, which helps to understand how the values in
the past can be used to predict the values in the future. If the mechanism
generating the data is deterministic, but the time series behaves like
generated by random processes, then the system is considered as chaotic. A
chaotic system is deterministic but not predictable in the long range. In a
deterministic system which is not chaotic, the value of a variable at a given
time can be used to generate the value of that variable at all times in
future.

Box counting technique is generally used to determine the fractal dimension of
the phase space, which is generated from the time series data. The values of
fractal dimension give the number of independent values to construct the time
series from which the phase space is generated. The box counting algorithm
analyzes the phase space set generated from the points with coordinates
$x(t)$, $x(t+\Delta t)$, . . . , $x(t+(n-1)\Delta t)$, where $x(t)$ are the
time series values and $\Delta t$ is the lag.

Figure 1: $PQRST$ in ECG and phase space. (a) ECG (b) Phase space.

At early times, electric activity of heart was modeled by using a set of
nonlinear differential equations van der Pol Jun Docts. Sc. and van der Mark
(1928); Katholi et al. (1977); West et al. (1985). Afterwards, variations of
cardiac spectra originated from deterministic dynamics called "chaotic
dynamics" which is highly sensitive to heart’s initial conditions are analyzed
by Babloyantz et al. (1985); Babloyantz and Destexhe (1986). Analyzing the
power spectra of a normal human heart shows that His-Purkinje network
possesses a fractal like structure in L.Goldberger et al. (1985) and the
existence of chaotic behavior of heart can be found out from the phase space
picture by evaluating a fractal dimension, $D$, where phase space trajectories
are extracted from the time series graph of ECG (Figure 1(b)) Babloyantz and
Destexhe (1988). Lower values of $D$ indicate more coherent dynamics. If
$D=1$, the oscillation is periodic and the phase space picture shows limiting
cycle. However, $D$ becomes larger than one when a limiting cycle is perturbed
by random noise Bergé et al. (1984). $D$ has non integer values greater than
two when the system becomes chaotic or strange attractor Bergé et al. (1984).
In this case, although trajectories in time do not converge towards a limiting
cycle, they stay in a bounded region in the phase space. Instead of $D$,
Babloyantz and Destexhe (1988) evaluates the correlation dimension $D_{2}$
from a time series of finite length using the existing algorithms, Grassberger
and Procaccia (1983a, b). The produced correlation dimensions are obtained
from a total of 36 ECG leads taken from 4 normal resting persons. Within the
range of computational errors, they find values of $D_{2}$ ranging from $3.6$
to $5.2$. These values suggest that the normal cardiac oscillations follow a
deterministic dynamics of chaotic nature.

It is also shown that a short-term heart rate variability analysis yields a
prognostic value in risk stratification, independent of clinical and
functional variables, Rovere et al. (2003). However, the detailed description
and classification of dynamical changes using time and frequency measures are
often not sufficient, especially in dynamical diseases as characterized by
Mackey and Glass (1977, 1979).

Wessel et al. (2009) try to answer the question: is the normal heart rate
chaotic due to respiration? In their work, they give an example of the
influence of respiration on heart beat dynamics, showing that observed
fluctuations can mostly be explained by respiratory modulations of heart rate
and blood pressure. Recently, the development of a fast and robust method
which can be applied to multichannel physiologic signals was reported by
Wilson and Haueisen (2017). This method elaborates either removing a selected
interfering signal or separating signals that arise from temporally correlated
and spatially distributed signals such as maternal or fetal ECG spectra.
Convolutional neural networks (CNNs) method was also applied to patient
specific ECG classification for real-time heart monitoring, Kiranyaz et al.
(2016).

Nowadays, it is well understood that cardiac diseases are one of the main
reasons of mortality in modern, industrialized societies, and they cause high
expenses in public health systems. Therefore, it is important to develop
analytical methods to improve cardiac diagnostics. In this work, we
investigate two different ECG systems taken from normal and abnormal human
hearts, Goldberger et al. (2000). Our aim is to introduce auxiliary phase
space method in conjunction with ECG signals to diagnose heart diseases. We
apply fractal analysis to the given data through trajectories produced in
phase space, from where fractal dimension $D$ is obtained with the use of the
box counting method, Liebovitch and Toth (1989).

##  2 Methods

###  2.1 Data

The data are taken from European $ST$-$T$ Database intended to be used for the
evaluation of algorithms for analysis of $ST$ and $T$-wave changes, Goldberger
et al. (2000). We have selected three different ECG records of two persons.
Person 1, considered as having a normal heart, is a man aged 51 with resting
angina and normal coronary arteries. Person 2, considered as having an
abnormal heart, is a man aged 58 with resting angina, anterior myocardial
infarction, 1-vessel disease (LAD) and aortic valvular regurgitation. The ECG
records e0118, e0121 and e0122 of person 1 (Figure 2) and the ECG records
e0123, e0125 and e0126 of person 2 (Figure 3) are examined. Each record has
two signals registered based on lead $V_{4}$ and modified lead III (MLIII),
respectively. For each signal, 200,000 samples are used.

Figure 2: First and Second ECG signal of normal heart taken at two different
records of the same person. Figure 3: First and Second ECG signal of abnormal
heart taken at two different records of the same person.

###  2.2 Algorithms

####  2.2.1 Phase space

We construct the phase space based on heart voltage values over time (i.e.,
$V(t)$) and their first derivative (i.e., $dV(t)/dt$). Figure 4, and Figure 5
show the phase spaces, $V(t)$ versus $dV(t)/dt$.

To obtain the first derivative of the function $V(t)$, we use third-order
forward difference Taylor Series derivative approximation of $f^{\prime}(x)$
Burden and Faires (1993); Khan and Ohba (1999); Ronco et al. (1999). The
simple approximation of the first derivative of a function $f$ at a point $x$
is defined as the limit of a difference quotient as follows:

$f^{\prime}(x)=\lim_{h\rightarrow 0}\frac{f(x+h)-f(x)}{h}$ (1)

Figure 4: Phase space of normal heart for the first and the second signal.
Figure 5: Phase space of abnormal heart for the first and the second signal.

If $h>0$, meaning that $h$ is a finite positive number, then

$f^{\prime}(x)=\frac{f(x+h)-f(x)}{h}$ (2)

is called the first-order forward difference approximation of $f^{\prime}(x)$.
The approximation of $f^{\prime}(x)$ can be obtained by combining Taylor
series expansions. If the values ($x_{(i+1)}$,$f_{(i+1)}$),
($x_{(i+2)}$,$f_{(i+2)}$) and ($x_{(i+3)}$,$f_{(i+3)}$) are known, the first
order derivative of $f_{i}(x)$ can be calculated as shown by
$f^{\prime}_{i}(x)$.

$f_{(i+1)}=f_{i}+\frac{f^{\prime}_{i}}{1!}h+\frac{f^{\prime\prime}_{i}}{2!}h^{2}+\frac{f^{\prime\prime\prime}_{i}}{3!}h^{3}+\ldots+\frac{f_{i}^{n}}{n!}h^{n}+R_{n}$
(3) $(x_{i+1},f_{i+1})\rightarrow
f_{i+1}=f_{i}+\frac{f^{\prime}_{i}}{1!}h+\frac{f^{\prime\prime}_{i}}{2!}h^{2}\rightarrow
f_{i+1}=f_{i}+hf^{\prime}_{i}+\frac{h^{2}}{2}f^{\prime\prime}_{i}$ (4)
$(x_{i+2},f_{i+2})\rightarrow
f_{i+2}=f_{i}+\frac{f^{\prime}_{i}}{1!}(2h)+\frac{f^{\prime\prime}_{i}}{2!}(2h)^{2}\rightarrow
f_{i+2}=f_{i}+2hf^{\prime}_{i}+2h^{2}f^{\prime\prime}_{i}$ (5)
$(x_{i+3},f_{i+3})\rightarrow
f_{i+3}=f_{i}+\frac{f^{\prime}_{i}}{1!}(3h)+\frac{f^{\prime\prime}_{i}}{2!}(3h)^{2}\rightarrow
f_{i+3}=f_{i}+3hf^{\prime}_{i}+\frac{9}{2}h^{2}f^{\prime\prime}_{i}$ (6)

The second derivative $f^{\prime\prime}_{i}$ is canceled by adding Equations
4, 5 and 6 after multiplying them by $18$, $-9$, and $2$, respectively.
Therefore, we obtain:

$18f_{i+1}-9f_{i+2}+2f_{i+3}=11f_{i}+6hf^{\prime}_{i}$ (7)

Finally, third-order forward difference Taylor series approximation of
$f^{\prime}(x)$ is obtained as follows:

$f^{\prime}_{i}=\frac{1}{6h}(-11f_{i}+18f_{i+1}-9f_{i+2}+2f_{i+3})$ (8)

####  2.2.2 Box counting

The pseudocode of the box counting algorithm is shown in Figure 6. The samples
of the signal constitute the $x$ coordinates of the time series. A $y$
coordinate is calculated based on four consecutive $x$ values as explained in
the previous section (line 3). $x$ and $y$ values are then normalized (lines
6-9). A rectangle is indicated based on the ($x$,$y$) coordinates of its left
top corner, width and height (e.g., line 10). The for loop in lines 11-35
gives the number of boxes containing points in the phase space. The algorithm
starts with the smallest rectangle containing all points in the space (line
10). In each iteration, rectangles in the set Rectangle are divided into four
rectangles (lines 12-24). In lines 25-33, we count the number of boxes
containing at least one point in the space. The number of boxes per iteration
are given at the end (line 34).

Figure 6: Box Counting Algorithm.

##  3 Results and discussions

Figure 4, and 5 present the phase space pictures of ECG (time series) produced
from Figure 2, and Figure 3 for normal and abnormal hearts, respectively.
Phase space pictures of normal hearts for the first signal in Figure 4 show
deterministic behavior, meaning that this signal possesses a perfect
oscillatory character.

When the phase space pictures in Figure 4 are compared with the phase space
pictures of the first signal of the abnormal heart in Figure 5, random
behaviors are observed. Phase space pictures in Figure 5 imply that some
perturbations start to develop on the top of the heart oscillations at the
same time, which makes it difficult to understand the individual mechanisms
that occur while producing the first signal of the abnormal heart. Figure 4
(d) to (f) present the second signal of the normal heart, which shows slight
deviation from the oscillatory behavior, but still obeys the deterministic
character. The broadening of this signal is most probably due to the resting
angina which effects the heart oscillation. A similar broadening may be
interpreted for the abnormal heart. However, as seen in Figure 5, all the
phase space pictures have a strong random character, implying that the
abnormal heart has serious heart failure due to anterior myocardial
infarction, 1-vessel disease (LAD) and Aortic valvular regurgitation. In other
words, in phase space terminology, numerous simultaneous factors are behind
this random behavior of the abnormal heart. Here, in short, these phase space
pictures can imply that the abnormal heart is unable to perform normal
oscillatory action due to its cardiac deficiencies.

Figure 4, and Figure 5 can also be interpreted with the help of the peak
height, V(t) patterns in ECG in Figure 2,and Figure 3. In phase space
pictures, $R$ peaks always have larger values than $P$, $Q$, $S$ and $T$
peaks. On the other hand, $S$ peaks have the smallest values only in the first
signal of both normal and abnormal hearts, while in the second signals $Q$
peaks have the smallest values compared to the other peaks. Moreover, the
phase space picture in Figure 7(b) provides well localized $R$, $P$, $T$, $S$
and $Q$ values, implying that peak heights of the ECG signals in Figure 2(b)
have almost the same values and are distinguished from each other by their
location in the phase space picture. Spreading of $R$, $P$, $T$, $S$ and $Q$
values are more pronounced for the second signal’s phase space pictures in
Figure 4(e) due to a single defect, i.e. resting angina in normal heart. On
the other hand, as seen in Figure 5(a-b-c), especially the spreading of the R
values reflects the randomness of peak heights of $R$ in the ECG pattern,
implying a serious heart failure due to anterior myocardial infarction,
1-vessel disease (LAD) and Aortic valvular regurgitation of the abnormal
heart. This randomness of the peak heights of $R$, $P$, $T$, $S$ and $Q$
tremendously increases for the second signal of the abnormal heart, as seen in
Figure 5(d-e-f), where $R$ and $P$, $T$, $Q$, $S$ values spread in all
directions and become indistinguishable from each other in phase space,
showing a strong random character. In that sense, one of the significant
outcomes of our work is that the phase space analysis can be useful for the
diagnosis of heart diseases in conjunction with ECG patterns, because the
broadening of $R$, $P$, $T$, $S$ and $Q$ values can easily imply the
irregularities in ECG patterns and provides information for peak height
distribution.

The data in Figure 4, and Figure 5 are elaborated using box counting method in
Figure 6 , where the following equation (9) is used

$D=\frac{\log N}{\log\displaystyle\frac{1}{r}}$ (9)

to produce fractal dimension, $D$. In this equation, $N$ is the minimum number
of boxes needed to cover the set of points and $r$ is the box size. The plots
of $\log N$ versus $\log r$ are given in Figure 7 and Figure 8, from where $D$
values are calculated.

Figure 7: $\log N$ versus $\log r$ plots and best fit for normal heart data.
Figure 8: $\log N$ versus $\log r$ plots and best fit for abnormal heart data.

The results are shown in Table I together with the correlation coefficients,
$R^{2}$ which are found out to be in a reasonable range.

Fractal dimensions, $D$, calculated from the phase space pictures in Figure 4,
and Figure 5 by using the box counting method are listed in Table I. $D$
values calculated from the normal heart data (i.e., 1.787, 1.749, 1.708 for
the first signal and 1.804, 1.816, 1.821 for the second signal) are smaller
than the abnormal heart dimensions (i.e., 1.816, 1.814, 1.816 for the first
signal and 1.863, 1.861, 1.860 for the second signal), supporting the existing
deterministic behavior of the normal heart compared to the random behavior of
the abnormal heart. In other words, the former is a good oscillator, while the
latter has some serious difficulties to make perfect harmonic oscillations.

Table 1: Fractal Dimensions Produced from Phase Space in Figure 4, and Figure 5 and Fitting Procedures in Figure 7, and Figure 8. |  | Fr.Dim.$D$ |  | Fitting $R^{2}$ |
---|---|---|---|---|---
| Patient | Signal1 | Signal2 | Signal1 | Signal2
Normal | e0118 | 1.787 | 1.804 | 0.9981 | 0.9974
| e0121 | 1.749 | 1.816 | 0.9992 | 0.9980
| e0122 | 1.708 | 1.821 | 0.9992 | 0.9987
Abnormal | e0123 | 1.816 | 1.863 | 0.9978 | 0.9983
| e0125 | 1.814 | 1.861 | 0.9982 | 0.9980
| e0126 | 1.816 | 1.860 | 0.9982 | 0.9977

##  4 Conclusion

The most crucial observation in this study is the behavior of the second
signals (MLIII) which gives more information than the first signals ($V_{4}$)
for the action of normal and abnormal hearts. The second signals have larger
fractal dimension $D$ values than the first signals, predicting more
randomness yet more information about MLIII measurements. The lowest value of
$D$ (i.e., 1.708) indicates a perfect oscillation of the heart and the highest
value of $D$ (i.e., $1.863$) presents randomness of the heart. In fact, phase
space picture presents the distribution of the peak heights in the ECG
spectra, giving valuable information about heart activities in conjunction
with ECG itself. In future work, we plan to apply, both fractal dimension and
peak height distribution analysis in phase space for various abnormal human
hearts to improve the novel diagnoses method for heart diseases.

#####

Acknowledgment

We would like to thank Dr. Eliya Buyukkaya from Wageningen University for
visualizing the ECG data as well as for her fruitful discussions with us.

#####

Conflict of interest

The author declares no conflict of interest.

#####

Author’s contribution

The authors contributed equally to this work, and read and approved the final
manuscript.

#####

Availability of supporting data

All data used in this paper are obtained from European ST-T Database. It is
available at https://physionet.org/physiobank/database/edb/.

#####

Ethical approval and consent to participate

Kadir Has University sees no Ethical conflict about this study and approves
this manuscript to be submitted to scientific journals. Besides the need for
University consent was waived.

## References

  * Liebovitch [1998] Larry S. Liebovitch.  _Fractals and chaos simplified for the life sciences_.  New York: Oxford University Press, 1998.
  * van der Pol Jun Docts. Sc. and van der Mark [1928] Balth van der Pol Jun Docts. Sc. and J. van der Mark.  Lxxii. the heartbeat considered as a relaxation oscillation, and an electrical model of the heart.  _Philosophical Magazine Series 1_ , 6:763–775, 1928.
  * Katholi et al. [1977] Charles R. Katholi, F. Urthaler, J. Macy, and T.N. James.  A mathematical model of automaticity in the sinus node and av junction based on weakly coupled relaxation oscillators.  _Computers and Biomedical Research_ , 10(6):529–543, 1977.  ISSN 0010-4809.  doi:https://doi.org/10.1016/0010-4809(77)90011-8.  URL https://www.sciencedirect.com/science/article/pii/0010480977900118.
  * West et al. [1985] Bruce J. West, Ary L. Goldberger, Galina Rovner, and Valmik Bhargava.  Nonlinear dynamics of the heartbeat: I. the av junction: Passive conduit or active oscillator?  _Physica D: Nonlinear Phenomena_ , 17(2):198–206, 1985.  ISSN 0167-2789.  doi:https://doi.org/10.1016/0167-2789(85)90004-1.  URL https://www.sciencedirect.com/science/article/pii/0167278985900041.
  * Babloyantz et al. [1985] Agnessa Babloyantz, J.M. Salazar, and C. Nicolis.  Evidence of chaotic dynamics of brain activity during the sleep cycle.  _Physics Letters A_ , 111(3):152–156, 1985.  ISSN 0375-9601.  doi:https://doi.org/10.1016/0375-9601(85)90444-X.  URL https://www.sciencedirect.com/science/article/pii/037596018590444X.
  * Babloyantz and Destexhe [1986] Agnessa Babloyantz and Alain Destexhe.  Low-dimensional chaos in an instance of epilepsy.  _Proceedings of the National Academy of Sciences of the United States of America_ , 83(10):3513–3517, 1986.
  * L.Goldberger et al. [1985] Ary L.Goldberger, Valmik Bhargava, Bruce J. West, and Arnold J. Mandell.  On a mechanism of cardiac electrical stability. the fractal hypothesis.  _Biophysical journal_ , 48(3):525–528, 1985.  doi:https://doi.org/10.1016/S0006-3495(85)83808-X.
  * Babloyantz and Destexhe [1988] Agnessa Babloyantz and Alain Destexhe.  Is the normal heart a periodic oscillator?  _Biological Cybernetics_ , 58(3):203–211, 1988\.  doi:https://doi.org/10.1007/BF00364139.
  * Bergé et al. [1984] Pierre Bergé, Yves Pomeau, and Christian Vidal.  _L’ordre dans le chaos: Vers une approche déterministe de la turbulence_.  Hermann, 1984.
  * Grassberger and Procaccia [1983a] Peter Grassberger and Itamar Procaccia.  Measuring the strangeness of strange attractors.  _Physica D: Nonlinear Phenomena_ , 9(1):189–208, 1983a.
  * Grassberger and Procaccia [1983b] Peter Grassberger and Itamar Procaccia.  Estimation of the kolmogorov entropy from a chaotic signal.  _Physical Review A_ , 28:2591–2593, 1983b.
  * Rovere et al. [2003] Maria Teresa La Rovere, Gian Domenico Pinna, Roberto Maestri, Andrea Mortara, Soccorso Capomolla, Oreste Febo, Roberto Ferrari, Mariella Franchini, Marco Gnemmi, Cristina Opasich, Pier Giorgio Riccardi, Egidio Traversi, and Franco Cobelli.  Short-term heart rate variability strongly predicts sudden cardiac death in chronic heart failure patients.  _Circulation_ , 107(4):565–570, 2003.
  * Mackey and Glass [1977] Michael C. Mackey and Leon Glass.  Oscillation and chaos in physiological control systems.  _Science_ , 197(4300):287–289, 1977.
  * Mackey and Glass [1979] Michael C. Mackey and Leon Glass.  Pathological conditions resulting from instabilities in physiological control systems.  _Annals of the New York Academy of Sciences_ , 316:214–235, 1979.
  * Wessel et al. [2009] Niels Wessel, Maik Riedl, and Jurgen Kurths.  Is the normal heart rate ”chaotic” due to respiration?  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 19(2), 2009.
  * Wilson and Haueisen [2017] James D. Wilson and Jens Haueisen.  Separation of physiological signals using minimum norm projection operators.  _IEEE Transactions on Biomedical Engineering_ , 64(4):904–916, 2017.
  * Kiranyaz et al. [2016] Serkan Kiranyaz, Turker Ince, and Moncef Gabbouj.  Real-time patient-specific ecg classification by 1-d convolutional neural networks.  _IEEE Transactions on Biomedical Engineering_ , 63(3):664–675, 2016.
  * Goldberger et al. [2000] Ary L. Goldberger, Luis A. N. Amaral, Leon Glass, Jeffrey M. Hausdorff, Plamen Ch. Ivanov, Roger G. Mark, Joseph E. Mietus, George B. Moody, Chung-Kang Peng, and H. Eugene Stanley.  Physiobank, physiotoolkit, and physionet: Components of a new research resource for complex physiologic signals.  _Circulation_ , 101(23), 2000.
  * Liebovitch and Toth [1989] Larry S. Liebovitch and Tibor Toth.  A fast algorithm to determine fractal dimensions by box counting.  _Physics Letters A_ , 141(8):386–390, 1989.
  * Burden and Faires [1993] Richard L. Burden and J. Douglas Faires.  _Numerical Analysis_.  PWS-Kent Pub. Co., Boston, 1993.
  * Khan and Ohba [1999] Ishtiaq Rasool Khan and Ryoji Ohba.  Closed form expressions for the finite difference approximations of first and higher derivatives based on taylor series.  _J. Comp. Appl. Math._ , 107:179–193, 1999.
  * Ronco et al. [1999] Eric Ronco, Taner Arsan, and Peter J. Gawthrop.  Open-loop intermittent feedback control: Practical continuous-time gpc.  _IEE Proceedings - Control Theory and Applications_ , 146(5):426–434, 1999.