# A Novel Approach in Solving Stochastic Generalized Linear Regression via
Nonconvex Programming

Vu Duc Anh1 , Tran Anh Tuan2∗, Tran Ngoc Thang2, Nguyen Thi Ngoc Anh2
1Nanyang Technological University, Singapore
2Faculty of Mathematics and Informatics, Hanoi University of Science and
Technology, Hanoi, Vietnam
<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>
Co-first authors, contributed equally Corresponding author

###### Abstract

Generalized linear regressions, such as logistic regressions or Poisson
regressions, are long-studied regression analysis approaches, and their
applications are widely employed in various classification problems. Our study
considers a stochastic generalized linear regression model as a stochastic
problem with chance constraints and tackles it using nonconvex programming
techniques. Clustering techniques and quantile estimation are also used to
estimate random data’s mean and variance-covariance matrix. Metrics for
measuring the performance of logistic regression are used to assess the
model’s efficacy, including the F1 score, precision score, and recall score.
The results of the proposed algorithm were over 1 to 2 percent better than the
ordinary logistic regression model on the same dataset with the above
assessment criteria.

_K_ eywords nonconvex programming, stochastic regression, generalized linear
regression, clustering, quantile

##  1 Introduction

In the realm of regression analysis, generalized linear models (GLMs) have
carved out a significant niche, particularly in solving classification
challenges. This paper delves into an advanced application of GLMs, focusing
on the stochastic generalized linear regression model. Unlike traditional
methods, this study approaches the model as a stochastic problem, integrated
with chance constraints. Our methodology is rooted in the use of nonconvex
programming techniques, a departure from conventional convex approaches.

We integrate clustering techniques and quantile estimation to further enhance
the model’s robustness. These methods are instrumental in accurately
estimating the mean and variance-covariance matrix of random data, a crucial
step in analyzing stochastic problems. Our study applies these techniques and
innovatively combines them with the stochastic generalized linear regression
model.

A critical aspect of our research is evaluating the model’s performance. In
this context, we employ various metrics specifically tailored for logistic
regression. These include the F1 score, precision score, and recall score.
These metrics comprehensively evaluate the model’s effectiveness in real-world
scenarios.

Logistic regression stands out as a crucial analytical method, as demonstrated
by [1], who illustrated the construction and updating of a logistic model for
tackling classification issues. [2] applied logistic regression in the context
of classifying multiple image classes, serving as a classification algorithm
for predicting labels within a discrete set of classes. Its widespread
applications include identifying email spam (as investigated by [3]),
detecting internet fraud (explored by [4]), and forecasting short-term
electricity load (explored by [5]).

To enhance the effectiveness of logistic regression in handling large datasets
with speed and precision, robust tools have been employed. Various convex
programming algorithms play a role in finding optimal solutions, as discussed
by [6], who also delved into the theoretical foundations of optimization
problem-solving methods. Challenges arise in maintaining the consistency of
problem constraints in the face of stochastic events affecting attribution
data. To address this, [7], [8], and [9] introduced chance-constrained
programming specifically designed to overcome such difficulties. In the
pursuit of solving optimization problems related to logistic regression, the
CVXOPT package by L. Vandenberghe for convex optimization problems has proven
instrumental.

The structure of our research paper is methodically organized as follows: In
the second section of our investigation, we provide some general preliminary
findings. Section 3 contains the specifics of our suggested technique. Within
Section 4, we provide the findings, and within Section 5, we talk about the
advantages of our approach. In the sixth section, we provide our concluding
observations as well as instructions for further research.

##  2 Preliminaries

The general loss function of the logistic regression model which is proposed
by [2]:

$\displaystyle-\log\bigg{(}P(\mathbf{y}|\mathbf{X},\mathbf{w})\bigg{)}=-$
$\displaystyle\sum_{i=1}^{\mathbf{N}}\bigg{(}y_{i}\log(\widehat{y_{i}})+(1-y_{i})\log(1-\widehat{y_{i}})\bigg{)}$
(1)

where
$\mathbf{X}=(x^{(1)},x^{(2)},\dots,x^{(\mathbf{N})}),x^{(i)}\in\mathbb{R}^{\mathbf{d}}$
$\mathbf{y}=(y_{1},y_{2},\dots,y_{N}),y_{i}\in\\{0,1\\}$,
$\mathbf{w}=(w_{1},w_{2},\dots,w_{d}),\mathbf{w}\in\mathbb{R}^{\mathbf{d}}$
$\widehat{y}_{i}=\theta(\mathbf{w}^{T}x^{(i)}),\theta$ is activation function
$\mathbf{d}$ is the number of variables, and $\mathbf{N}$ is the number of
observations.

Loss function (1) represent the cross-entropy loss between $y_{i}$ and
$\widehat{y_{i}}$. Furthermore, [2] showed also the loss function of the
logistic sigmoid regression model which is modelized by sigmoid function

$\displaystyle J(\mathbf{w})=$
$\displaystyle-\sum_{i=1}^{\mathbf{N}}\bigg{(}y_{i}\log\left(\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\right)+(1-y_{i})\log\left(1-\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\right)\bigg{)},$
(2)

where
$\mathbf{X}=(x^{(1)},x^{(2)},\dots,x^{(\mathbf{N})}),x^{(i)}\in\mathbb{R}^{\mathbf{d}}$
$\mathbf{y}=(y_{1},y_{2},\dots,y_{N}),y_{i}\in\\{0,1\\}$
$\mathbf{w}=(w_{1},w_{2},\dots,w_{d}),\mathbf{w}\in\mathbb{R}^{\mathbf{d}}$
$\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)=\dfrac{1}{1+e^{(-\mathbf{w}^{T}x^{(i)})}},$
is sigmoid function $\mathbf{d}$ is the number of variables $\mathbf{N}$ is
the number of observations.

Figure 1: The general architecture of proposed model.

######  Definition 2.1.

([10]). Consider a convex set $\Omega\subseteq\mathbb{R^{n}}$, a continuously
function $J(w)$ is said to be convex on $\Omega$ if, with $\forall
w^{(1)},w^{(2)}\in\Omega\text{ and }\alpha\in\left[0,1\right]$, the following
inequality holds true:

$\displaystyle J\left(\alpha w^{(1)}+(1-\alpha)w^{(2)}\right)\leq\alpha
J(w^{(1)})+(1-\alpha)J(w^{(2)}),$

then, $-J(w)$ is a concave function.

######  Definition 2.2.

([11]). A continuously function $J:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is
said to be semi-strictly quasi-convex on a convex set $\Omega$ if for every
pair of distinct $w^{(1)},w^{(2)}\in\Omega,\lambda\in(0,1)$, we have

$\displaystyle J\left(\lambda w^{(1)}+(1-\lambda)w^{(2)}\right)\geq
J(w^{(1)})\Rightarrow J(w^{(2)})\geq J(w^{(1)}),$

then, $-J(w)$ is a semi-strictly quasi-concave function.

######  Definition 2.3.

([11]). With $\alpha\in\mathbb{R}$, a convex set
$\Omega\subseteq\mathbb{R^{n}}$, let

$\displaystyle U^{(\alpha)}(J)=\\{w\in\Omega|J(w)\geq\alpha\\}\text{ is upper-
level set of function J},$ $\displaystyle
L^{(\alpha)}(J)=\\{w\in\Omega|J(w)\leq\alpha\\}\text{ is lower-level set of
function J}.$

######  Corollary 2.1.

([10]). Let $J$ be a convex function on $\Omega\subset\mathbb{R}^{n}$. The
lower-level set of function $J$ is a convex set with
$\forall\alpha\in\mathbb{R}$.

######  Theorem 2.1.

([12]). A continuous function $J$ is semi-strictly quasi-convex on a convex
set $\Omega$ if and only if lower-level set $L^{(\alpha)}(J)$ is convex set
with $\forall\alpha\in\mathbb{R}$.

######  Theorem 2.2.

([11]). Let $\Phi_{0}$ be a positive convex function on
$\Omega\subset\mathbb{R}^{n}$ and $\Phi_{k}$ be non-negative and concave
functions on $\Omega$. Assume that $b\geq\sum_{i=1}^{k}a_{i}$ where $a_{i}>0.$
Then,

$\displaystyle F(x)=\Pi_{i=1}^{k}[\Phi_{i}(x)]^{a_{i}}/[\Phi_{0}(x)]^{b},$

is semistrictly quasiconcave on $\Omega$.

######  Definition 2.4.

([13]). A quadratic form $\xi(w)=w^{T}Qw$ is convex function on convex set
$\Omega\subseteq\mathbb{R}^{n}$, if only if $Q$ is positive semi-definite,
i.e.

$\displaystyle\xi(w)\text{ is convex function on }\Omega\text{ if only if
}\xi(w)\geq 0,\forall w\in\Omega.$

######  Theorem 2.3.

([14]). A function $J(w)$ which is twice-differentiable is convex if and only
if its hessian matrix (matrix of second-order partial derivatives) is positive
semi-definite, i.e.

$\displaystyle\text{with }\forall v:v^{T}\left(\nabla^{T}_{w}J(w)\right)v\geq
0,$

where, $\nabla^{T}_{w}J(w)$ is the hessian matrix.

######  Proposition 2.1.

([14]). Let $J(w)$ and $G(w)$ be two convex functions,
$\lambda_{1},\lambda_{2}$ are non-negative number. Then any linear combination
of these two functions

$\displaystyle\left(\lambda_{1}J+\lambda_{2}G\right)(w)=\lambda_{1}J(w)+\lambda_{2}G(w),$

is also a convex function.

######  Proposition 2.2.

([15]). Suppose $v:\mathbb{R}^{k}\to\mathbb{R}$,
$l:\mathbb{R}^{n}\to\mathbb{R}^{k}$, consider $f=v\circ
l:\mathbb{R}^{n}\to\mathbb{R}$. If both $v$ and $l$ are convex functions, $l$
is a non-decreasing function, then $f=v\circ l$ is convex function.

######  Lemma 2.1.

$J(\mathbf{w})$ in (2) is a convex function relatives to $\mathbf{w}$

###### Proof.

To prove $J(\mathbf{w})$ is convex function, we’ll prove two component
functions

$\displaystyle-\log\left(\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\right)\text{
and }-\log\left(1-\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\right),$

are convex functions relatives to $\mathbf{w}$.
Firstly, we have

$\displaystyle-\log\bigg{(}\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\bigg{)}$
$\displaystyle=-log\left(\dfrac{1}{1+e^{(-\mathbf{w}^{T}x^{(i)})}}\right)$
$\displaystyle=log\left(1+e^{(-\mathbf{w}^{T}x^{(i)})}\right),$

has the hessian matrix

$\displaystyle\nabla^{2}_{\mathbf{w}}\bigg{[}-\log\bigg{(}\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\bigg{)}\bigg{]}$
$\displaystyle=\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\left(1-\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\right)x^{(i)}\left(x^{(i)}\right)^{T}.$

Now, we use Theorem 2.3

$\displaystyle\forall
v:v^{T}\nabla^{2}_{\mathbf{w}}\bigg{[}-\log\bigg{(}\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\bigg{)}\bigg{]}v$
$\displaystyle=\bigg{[}\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\left(1-\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\right)x^{(i)}\left(x^{(i)}\right)^{T}\bigg{]}$
$\displaystyle=\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\left(1-\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\right)\bigg{[}\left(x^{(i)}\right)^{T}v\bigg{]}^{2}\geq
0,$

hence, $-\log\bigg{(}\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\bigg{)}$ is
convex function.
Secondly, we have also

$\displaystyle-\log\bigg{(}1-\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\bigg{)}=-log\left(1-\dfrac{1}{1+e^{(-\mathbf{w}^{T}x^{(i)})}}\right)$
$\displaystyle=\mathbf{w}^{T}x^{(i)}+log\left(1+e^{(-\mathbf{w}^{T}x^{(i)})}\right),$

has the hessian matrix

$\displaystyle\nabla^{2}_{\mathbf{w}}\bigg{[}-\log\bigg{(}1-\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\bigg{)}\bigg{]}$
$\displaystyle=\nabla^{2}_{\mathbf{w}}\bigg{[}-\log\bigg{(}\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\bigg{)}\bigg{]}$
$\displaystyle=\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\left(1-\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\right)x^{(i)}\left(x^{(i)}\right)^{T},$

use above proof, we have
$-\log\bigg{(}1-\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\bigg{)}$ is also
convex function.
Last, use Proposition 2.1 with $y_{i}\geq 0$ and $1-y_{i}\geq 0$, we obtain

$\displaystyle-\sum_{i=1}^{\mathbf{N}}\bigg{(}y_{i}\log\left(\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\right)+(1-y_{i})\log\left(1-\sigma\left(-\mathbf{w}^{T}x^{(i)}\right)\right)\bigg{)},$

is convex function relatives to $\mathbf{w}$. ∎

##  3 Proposed Methodology

###  3.1 Main Methodology

The fundamental idea is to combine additional information from the mean and
variance of data, then utilize probability formulas to transform a
deterministic problem into a problem with chance constraints.
Let $\varepsilon_{i}=\mathbf{w}^{T}x^{(i)}$, modelize minimize loss function
(2) into problem ($P_{1}$):

$\displaystyle-\sum_{i=1}^{\mathbf{N}}\bigg{(}y_{i}\log\left(\sigma\left(-\varepsilon_{i}\right)\right)+(1-y_{i})\log\left(1-\sigma\left(-\varepsilon_{i}\right)\right)\bigg{)}$
($P_{1}$) $\displaystyle
s.t:\begin{cases}P\bigg{(}|\mathbf{w}^{T}\overline{x}^{(i)}-\varepsilon_{i}|\leq\alpha_{i}\bigg{)}\geq\beta_{i}\\\
E\bigg{[}\mathbf{w}^{T}\overline{x}^{(i)}-\varepsilon_{i}\bigg{]}=0\end{cases},$

where $\alpha_{i},\dots,\alpha_{N}$ are the non-negative real numbers and
$\beta_{1},\beta_{2},\dots,\beta_{N}$ are the probability level of each
constraint occurring .
Use probability formulas, we have:

$\displaystyle
P\bigg{(}|\mathbf{w}^{T}\overline{x}^{(i)}-\varepsilon_{i}|\leq\alpha_{i}\bigg{)}=P\bigg{(}-\alpha_{i}\leq\mathbf{w}^{T}\overline{x}^{(i)}-\varepsilon_{i}\leq\alpha_{i}\bigg{)}$
$\displaystyle=P\bigg{(}\mathbf{w}^{T}\overline{x}^{(i)}-\varepsilon_{i}\leq\alpha_{i}\bigg{)}-P\bigg{(}\mathbf{w}^{T}\overline{x}^{(i)}-\varepsilon_{i}\leq-\alpha_{i}\bigg{)}$
$\displaystyle=2P\bigg{(}\mathbf{w}^{T}\overline{x}^{(i)}-\varepsilon_{i}\leq\alpha_{i}\bigg{)}-1\geq\beta_{i},$

hence,

$\displaystyle
P\bigg{(}\mathbf{w}^{T}\overline{x}^{(i)}-\varepsilon_{i}\leq\alpha_{i}\bigg{)}\geq\dfrac{\beta_{i}+1}{2}$
$\displaystyle
P\bigg{(}\mathbf{w}^{T}\overline{x}^{(i)}-\varepsilon_{i}\leq-\alpha_{i}\bigg{)}\leq\dfrac{1-\beta_{i}}{2},$

with,

$\displaystyle-\sum_{i=1}^{\mathbf{N}}\bigg{(}y_{i}\log\left(\sigma\left(-\varepsilon_{i}\right)\right)+(1-y_{i})\log\left(1-\sigma\left(-\varepsilon_{i}\right)\right)\bigg{)}$
$\displaystyle=-\sum_{i=1}^{\mathbf{N}}\bigg{(}y_{i}\varepsilon_{i}-\log{\left(e^{\varepsilon_{i}}+1\right)}\bigg{)}.$

We rewrote problem ($P_{1}$) into problem ($P_{2}$):

$\displaystyle\underset{[\mathbf{w},\varepsilon]}{\min}J(\mathbf{w},\varepsilon)=-\sum_{i=1}^{\mathbf{N}}\bigg{(}y_{i}\varepsilon_{i}-\log{\left(e^{\varepsilon_{i}}+1\right)}\bigg{)}$
($P_{2}$) $\displaystyle
s.t:\begin{cases}P\bigg{(}\mathbf{w}^{T}\overline{x}^{(i)}-\varepsilon_{i}\leq\alpha_{i}\bigg{)}\geq\dfrac{\beta_{i}+1}{2}\\\
\\\
P\bigg{(}\mathbf{w}^{T}\overline{x}^{(i)}-\varepsilon_{i}\leq-\alpha_{i}\bigg{)}\leq\dfrac{1-\beta_{i}}{2}\\\
E\big{[}\mathbf{w}^{T}\overline{x}^{(i)}-\varepsilon_{i}\big{]}=0.\end{cases}$

Next, we assumed $\overline{x}_{j}^{(i)}$ is a normal random variable with
mean is $m_{\overline{x}_{j}^{(i)}}$, we denoted:

  * •

$\mathbf{m}_{\overline{x}^{(i)}}=\bigg{(}m_{\overline{x}_{1}^{(i)}},\dots,m_{\overline{x}_{d}^{(i)}}\bigg{)}$
is sample mean of $\overline{x}^{(i)}$.

  * •

$V_{\overline{x}^{(i)}}$ is sample variance-covariance symmetric matrix of
$\overline{x}^{(i)}$ with shape $d\times d$.

  * •

$\phi$ is the cumulative distribution function of the standard
Gaussian$(0,1)$.

Then, we transformed chance constraints of problem ($P_{2}$) and utilize the
monotonic of the function $\phi$, it follows that:

$\displaystyle
P\bigg{(}\mathbf{w}^{T}\overline{x}^{(i)}-\varepsilon_{i}\leq\alpha_{i}\bigg{)}$
$\displaystyle=P\Bigg{(}\dfrac{\mathbf{w}^{T}\overline{x}^{(i)}-\mathbf{w}^{T}\mathbf{m}_{\overline{x}^{(i)}}}{\sqrt{\mathbf{w}^{T}V_{\overline{x}^{(i)}}\mathbf{w}}}\leq\dfrac{\alpha_{i}+\varepsilon_{i}-\mathbf{w}^{T}\mathbf{m}_{\overline{x}^{(i)}}}{\sqrt{\mathbf{w}^{T}V_{\overline{x}^{(i)}}\mathbf{w}}}\Bigg{)}$
$\displaystyle=\phi\Bigg{(}\dfrac{\alpha_{i}+\varepsilon_{i}-\mathbf{w}^{T}\mathbf{m}_{\overline{x}^{(i)}}}{\sqrt{\mathbf{w}^{T}V_{\overline{x}^{(i)}}\mathbf{w}}}\Bigg{)}\geq\dfrac{\beta_{i}+1}{2}$
$\displaystyle\Leftrightarrow\dfrac{\alpha_{i}+\varepsilon_{i}-\mathbf{w}^{T}\mathbf{m}_{\overline{x}^{(i)}}}{\sqrt{\mathbf{w}^{T}V_{\overline{x}^{(i)}}\mathbf{w}}}\geq\phi^{-1}\Bigg{(}\dfrac{\beta_{i}+1}{2}\Bigg{)}$
$\displaystyle\Leftrightarrow\phi^{-1}\Bigg{(}\dfrac{\beta_{i}+1}{2}\Bigg{)}-\dfrac{\alpha_{i}+\varepsilon_{i}-\mathbf{w}^{T}\mathbf{m}_{\overline{x}^{(i)}}}{\sqrt{\mathbf{w}^{T}V_{\overline{x}^{(i)}}\mathbf{w}}}\leq
0,$

similarly,

$\displaystyle
P\bigg{(}\mathbf{w}^{T}\overline{x}^{(i)}-\varepsilon_{i}\leq-\alpha_{i}\bigg{)}$
$\displaystyle=P\Bigg{(}\dfrac{\mathbf{w}^{T}\overline{x}^{(i)}-\mathbf{w}^{T}\mathbf{m}_{\overline{x}^{(i)}}}{\sqrt{\mathbf{w}^{T}V_{\overline{x}^{(i)}}\mathbf{w}}}\leq\dfrac{-\alpha_{i}+\varepsilon_{i}-\mathbf{w}^{T}\mathbf{m}_{\overline{x}^{(i)}}}{\sqrt{\mathbf{w}^{T}V_{\overline{x}^{(i)}}\mathbf{w}}}\Bigg{)}$
$\displaystyle=\phi\Bigg{(}\dfrac{-\alpha_{i}+\varepsilon_{i}-\mathbf{w}^{T}\mathbf{m}_{\overline{x}^{(i)}}}{\sqrt{\mathbf{w}^{T}V_{\overline{x}^{(i)}}\mathbf{w}}}\Bigg{)}\leq\dfrac{1-\beta_{i}}{2}$
$\displaystyle\Leftrightarrow\dfrac{-\alpha_{i}+\varepsilon_{i}-\mathbf{w}^{T}\mathbf{m}_{\overline{x}^{(i)}}}{\sqrt{\mathbf{w}^{T}V_{\overline{x}^{(i)}}\mathbf{w}}}-\phi^{-1}\Bigg{(}\dfrac{1-\beta_{i}}{2}\Bigg{)}\leq
0,$

and,

$\displaystyle
E\bigg{[}\mathbf{w}^{T}\overline{x}^{(i)}-\varepsilon_{i}\bigg{]}=\mathbf{w}^{T}\mathbf{m}_{\overline{x}^{(i)}}-\varepsilon_{i}=0.$

Thus, problem ($P_{2}$) equivalent with problem ($P_{3}$):

$\displaystyle\underset{[\mathbf{w},\varepsilon]}{\min}J(\mathbf{w},\varepsilon)=-\sum_{i=1}^{\mathbf{N}}\bigg{(}y_{i}\varepsilon_{i}-\log{\left(e^{\varepsilon_{i}}+1\right)}\bigg{)}$
($P_{3}$) $\displaystyle
s.t:\begin{cases}\phi^{-1}\Bigg{(}\dfrac{\beta_{i}+1}{2}\Bigg{)}-\dfrac{\alpha_{i}+\varepsilon_{i}-\mathbf{w}^{T}\mathbf{m}_{\overline{x}^{(i)}}}{\sqrt{\mathbf{w}^{T}V_{\overline{x}^{(i)}}\mathbf{w}}}\leq
0\\\ \\\
\dfrac{-\alpha_{i}+\varepsilon_{i}-\mathbf{w}^{T}\mathbf{m}_{\overline{x}^{(i)}}}{\sqrt{\mathbf{w}^{T}V_{\overline{x}^{(i)}}\mathbf{w}}}-\phi^{-1}\Bigg{(}\dfrac{1-\beta_{i}}{2}\Bigg{)}\leq
0\\\
\mathbf{w}^{T}\mathbf{m}_{\overline{x}^{(i)}}-\varepsilon_{i}=0.\end{cases}$

Or,

$\displaystyle\underset{[\mathbf{w},\varepsilon]}{\min}J(\mathbf{w},\varepsilon)=-\sum_{i=1}^{\mathbf{N}}\bigg{(}y_{i}\varepsilon_{i}-\log{\left(e^{\varepsilon_{i}}+1\right)}\bigg{)}$
($P_{3}$) $\displaystyle
s.t:\begin{cases}\phi^{-1}\Bigg{(}\dfrac{\beta_{i}+1}{2}\Bigg{)}-\dfrac{\alpha_{i}}{\sqrt{\mathbf{w}^{T}V_{\overline{x}^{(i)}}\mathbf{w}}}\leq
0\\\
\dfrac{-\alpha_{i}}{\sqrt{\mathbf{w}^{T}V_{\overline{x}^{(i)}}\mathbf{w}}}-\phi^{-1}\Bigg{(}\dfrac{1-\beta_{i}}{2}\Bigg{)}\leq
0.\\\ -\mathbf{w}^{T}\mathbf{m}_{\overline{x}^{(i)}}+\varepsilon_{i}\leq 0\\\
\mathbf{w}^{T}\mathbf{m}_{\overline{x}^{(i)}}-\varepsilon_{i}\leq
0.\end{cases}$

To make it convenient to follow, we let:

$\displaystyle\mathbf{Z}$
$\displaystyle=[Z_{1},Z_{2},\dots,Z_{N+d}]^{T}=[w_{0},\dots,w_{d},\varepsilon_{d+1},\dots,\varepsilon_{N+d}]^{T}(Z\neq
0),$

so, problem ($P_{3}$) was rewritten:

$\displaystyle\underset{\mathbf{Z}}{\min}J(\mathbf{Z})=-\sum_{i=1}^{\mathbf{N}}\Bigg{(}y_{i}Z_{d+i}-\log{\left(e^{Z_{d+i}}+1\right)}\Bigg{)}$
($P_{4}$) $\displaystyle s.t:$
$\displaystyle\begin{cases}\phi^{-1}\left(\dfrac{1+\beta_{i}}{2}\right)-\dfrac{\alpha_{i}}{\sqrt{\mathbf{Z}^{T}V_{[\overline{x}^{(i)},\mathbf{K}^{(i)}]}\mathbf{Z}}}\leq
0\\\ \\\
\dfrac{-\alpha_{i}}{\sqrt{\mathbf{Z}^{T}V_{[\overline{x}^{(i)},\mathbf{K}^{(i)}]}\mathbf{Z}}}-\phi^{-1}\left(\dfrac{1-\beta_{i}}{2}\right)\leq
0\\\ \mathbf{P_{i}Z}\leq 0\\\ -\mathbf{P_{i}Z}\leq 0,\end{cases}$

where,

  * •

$i=\overline{1,N}$.

  * •

$\mathbf{P_{i}}=\bigg{[}\mathbf{m}_{\overline{x}^{(i)}},0_{d+1},\dots,-1_{d+i},\dots,0_{d+N}\bigg{]}$
is the $i^{th}$ coefficient vectors.

  * •

$\phi$ is the cumulative distribution function of the standard Gaussian(0,1).

  * •

$\beta_{1},\dots,\beta_{N}$ are the probability level of each constraint
occurring.

  * •

$\mathbf{m}_{\overline{x}^{(i)}}=\left(m_{\overline{x}_{1}^{(i)}},\dots,m_{\overline{x}_{d}^{(i)}}\right)$
is sample mean of $\overline{x}^{(i)}$.

  * •

$V_{\bigg{[}\overline{x}^{(i)},\mathbf{K}^{(i)}_{(1\times
N)}\bigg{]}}=\begin{bmatrix}V_{\overline{x}^{(i)}}&\textbf{O}_{(N\times
N)},\\\ \textbf{O}_{(N\times N)}&\textbf{O}_{(N\times N)}\end{bmatrix}$ is
sample variance-covariance symmetric matrix of $\overline{x}^{(i)}$ with shape
$(N+d)\times(N+d)$.

  * •

$\mathbf{K}^{(i)}$ is vector $(0_{d+1},\dots,1_{d+i},\dots,0_{N+d})$.

  * •

O is zero matrix with shape $N\times N$.

######  Proposition 3.1.

([16]). A matrix is a variance-covariance matrix if only if it is positive
semi-definite.

######  Lemma 3.1.

The constraints set of problem ($P_{4}$) are convex set.

###### Proof.

Consider $\Omega=\\{Z\in\mathbb{R}^{N+d}:G_{i}^{(k)}(Z)\leq
0,k=\overline{1,4}\\}$ is feasible region (the constraints set) of problem
($P_{4}$).
Let

$\displaystyle
G_{i}^{(1)}(Z)=\phi^{-1}\left(\dfrac{1+\beta_{i}}{2}\right)-\dfrac{\alpha_{i}}{\sqrt{\mathbf{Z}^{T}V_{[\overline{x}^{(i)},\mathbf{K}^{(i)}]}\mathbf{Z}}}$
$\displaystyle
G_{i}^{(2)}(Z)=\dfrac{-\alpha_{i}}{\sqrt{\mathbf{Z}^{T}V_{[\overline{x}^{(i)},\mathbf{K}^{(i)}]}\mathbf{Z}}}-\phi^{-1}\left(\dfrac{1-\beta_{i}}{2}\right)$
$\displaystyle G_{i}^{(3)}(Z)=\mathbf{P_{i}Z}$ $\displaystyle
G_{i}^{(4)}(Z)=-\mathbf{P_{i}Z}.$

Now, we prove $G_{i}^{(1)}(Z)$ and $G_{i}^{(2)}(Z)$ are the semi-strictly
quasi-convex functions on $\Omega$.

In fact,
$\dfrac{\alpha_{i}}{\sqrt{\mathbf{Z}^{T}V_{[\overline{x}^{(i)},\mathbf{K}^{(i)}]}\mathbf{Z}}}=\dfrac{\sqrt[3]{\alpha_{i}^{3}}}{\sqrt{\mathbf{Z}^{T}V_{[\overline{x}^{(i)},\mathbf{K}^{(i)}]}\mathbf{Z}}}$
is a semi-strictly quasi-concave function on $\Omega$, cause of
$\alpha_{i}^{3}$ is a non-negative constant function. In addition, from
Proposition 3.1 and $Z\neq 0$ then $V_{[\overline{x}^{(i)},\mathbf{K}^{(i)}]}$
is positive definite, then use Definition 2.4 we implied
$\mathbf{Z}^{T}V_{[\overline{x}^{(i)},\mathbf{K}^{(i)}]}\mathbf{Z}$ is a
positive convex function. Thus, use Theorem 2.2, choose
$a_{1}=\dfrac{1}{3},b=\dfrac{1}{2}$, and collaborate with
$\phi^{-1}\left(\dfrac{1+\beta_{i}}{2}\right)$ is a constant function, we
imply $G_{i}^{(1)}(Z)$ is a semi-strictly quasi-convex functions on $\Omega$.
Similarity, we obtain $G_{i}^{(2)}(Z)$ is a semi-strictly quasi-convex
functions on $\Omega$.
From Theorem 2.3, we have

$\displaystyle
L^{(0)}\bigg{(}G_{i}^{(1)}(Z)\bigg{)}=\Bigg{\\{}\mathbf{Z}\in\Omega\Bigg{|}\phi^{-1}\left(\dfrac{1+\beta_{i}}{2}\right)-\dfrac{\alpha_{i}}{\sqrt{\mathbf{Z}^{T}V_{[\overline{x}^{(i)},\mathbf{K}^{(i)}]}\mathbf{Z}}}\leq
0\Bigg{\\}},$

and

$\displaystyle
L^{(0)}\bigg{(}G_{i}^{(2)}(Z)\bigg{)}=\Bigg{\\{}\mathbf{Z}\in\Omega\Bigg{|}\dfrac{-\alpha_{i}}{\sqrt{\mathbf{Z}^{T}V_{[\overline{x}^{(i)},\mathbf{K}^{(i)}]}\mathbf{Z}}}-\phi^{-1}\left(\dfrac{1-\beta_{i}}{2}\right)\leq
0\Bigg{\\}},$

are two convex sets.
Furthermore, $G_{i}^{(3)}(Z)$ and $G_{i}^{(4)}(Z)$ are two linear functions,
thus $G_{i}^{(3)}(Z)$ and $G_{i}^{(4)}(Z)$ are two convex functions.
Consequently, based on Corollary 2.1, then

$\displaystyle
L^{(0)}\bigg{(}G_{i}^{(3)}(Z)\bigg{)}=\\{Z\in\Omega|\mathbf{P_{i}Z}\leq 0\\},$

and

$\displaystyle
L^{(0)}\bigg{(}G_{i}^{(4)}(Z)\bigg{)}=\\{Z\in\Omega|-\mathbf{P_{i}Z}\leq
0\\},$

are two convex sets.
In short, the feasible region (the constraints set) $\Omega$ of problem
($P_{4}$) is a convex set. ∎

From Lemma 3.1 and Lemma 2.1, we implied problem ($P_{4}$) is a convex
programming problem. Hence, we used convex programming tools such as the
Lagrange multiplier method to tackle the solving of optimal solutions.

###  3.2 The main algorithm

In this section, we describe the main algorithms for our proposed model. But
before applying the main algorithms, we scored $y$ data by the original
Logistic Regression model and utilized clustering algorithms to estimate the
sample mean and variance-covariance matrix.

####  3.2.1 The sample mean and variance-covariance matrix estimation

K-means clustering and Q-quantiles clustering were used in the experimentation
with real datasets. We divided $X$ into groups using K-means clustering and
then calculated the mean and variance-covariance matrices for each group. We
also estimated the sample mean of $y$, which was in the appropriate group, as
part of this addition.

A. K-Means clustering algorithm

In the first approach, we used the K-means clustering algorithm. Clustering is
a popular tool for finding groups or clusters which have the same feature in
multivariate data and has found lots of applications in biology (see [17]),
medicine (see [18]), psychology, and economics ([19]).

We had difficulty finding the number of data in each cluster because of the
randomness of the cluster centers at initialization. An obligatory way for the
clustering algorithm is to ask for input on the number of clusters in advance,
which demands knowledge of the underlying datasets. K-Means is a simple
unsupervised learning algorithm that solves difficult clustering problems.
Hence, we must provide deterministic k clusters for solving a problem (see
[20]).

The variance-covariance matrix will not be computed if a cluster contains just
one data point. As a result, we chose a predetermined number of clusters to be
employed in the clustering process, ensuring that the number of data points in
each cluster is more than one. It is possible to estimate the median value and
the variance-covariance matrices by grouping the data from $X$ into a cluster.
The original logistic regression model will be used to translate $y$ into a
value in the interval $(0,1)$. Then, the value of $y$ is commensurate with $X$
in a cluster together is calculated by mean.

B. Q-Quantiles algorithm

The definition of quantile was proposed in the paper of [21]. On the contrary,
the K-Means clustering algorithm follows as Quantile estimation is divided
into quantile levels for data of $y$. Then, at the same level, the mean and
variance-covariance matrix of $X$ will be estimated. To estimate Large-scale
Data, we need to define Q-Quantiles, which are values that divide a finite set
into Q subsets of the same size. In our algorithm, we use a uniform
probability distribution, so Q-Quantiles have the values
$\\{1/q,2/q,\dots,q-1/q\\}$.

Moreover, instead of using a uniform probability distribution, we approached a
new method based on the ideas of [22] to divide quantile levels compatibly.
Similar to the K-Means clustering algorithm, the Q-Quantiles algorithm will be
implemented ineffectively if any quantile level contains the unit observation.

####  3.2.2 The pseudocode and procedure of the main algorithm

The pseudocode of the main algorithm was demonstrated by Algorithm 3.
Algorithm 1 is our proposed algorithm for solving the solution using convex
programming, its inputs include the sample mean, the sample variance-
covariance matrix of $X$ data, the sample mean of $y$ data, and
$\alpha,\beta$. Its output is an optimal solution to Problem $P_{4}$.

Algorithm 1 The pseudocode of the Stochastic Logistic Regression model (SLR)

1:Data input: $meanX,meanY,covX,\alpha,\beta$

2:The weight: $w$

3:Define the objective function ($F$) of Problem $P_{4}$;

4:Calculate constraint functions ($G$), and the coefficient of freedom ($h$)
from $meanX,meanY,covX$;

5:Use $CVXOPT$ package  $\displaystyle w=CVXOPT(F,G,h,dim);$

6:return $w$.

Algorithm 2 is the pseudocode of the Simulated Annealing algorithm for solving
the optimal parameters $\alpha,\beta$, its inputs consist of the sample mean,
the sample variance-covariance matrix of $X$ data, the sample mean of $y$
data, and $X\\_train,y\\_train$. Its outputs return the optimal $\alpha,\beta$
parameters, and the corresponding optimal solution.

Algorithm 2 The pseudocode of the Simulated Annealing algorithm (SA)

1:Data input: $meanX,meanY,covX,X\\_train,y\\_train$

2:$\alpha_{0},\beta_{0},w_{0}$

3:Set initial  $\displaystyle\alpha_{0}\geq 0,\beta_{0}\in(0,1),$
$\displaystyle w_{0}=SLR(meanX,covX,meanY,\alpha_{0},\beta_{0}),$
$\displaystyle f(w_{0})=-log\\_loss(w_{0},X\\_train,y\\_train);$

4:while The iteration isn’t satisfied do

5: Generate new initial  $\displaystyle\alpha_{new},\beta_{new},$
$\displaystyle w_{new}=SLR(meanX,covX,meanY,\alpha_{new},\beta_{new}),$
$\displaystyle f(w_{new})=-log\\_loss(w_{new},X\\_train,y\\_train);$

6: Calculate $\Delta f=f(w_{new})-f(w_{0})$;

7: if $\Delta f\leq 0,\alpha_{new}\geq 0,\beta_{new}\geq 0$ then

8: $w_{0}=w_{new},\alpha_{0}=\alpha_{new},\beta_{0}=\beta_{new}$,

9: if The stopping criteria is met then

10: Terminate the computation process;

11: else

12: The iteration isn’t satisfied;

13: end if

14: else

15: Accept the new solution based on Metropolis,

16: The iteration isn’t satisfied;

17: end if

18:end while

19:return $\alpha_{0},\beta_{0},w_{0}$.

Algorithm 3 The pseudocode of the main algorithm

1:Data input: $X,y$

2:$w_{opt},\alpha_{opt},\beta_{opt}$

3:Split, normalize data into train, and validation set.

4:Scoring $y\\_{train}$ data  $\displaystyle
coef=LogisticRegression(X\\_train,y\\_train),$ $\displaystyle
y\\_score=\dfrac{1}{1+e^{-coef^{T}X\\_train}};$

5:if Kmean then

6: Enter $n_{c}$ clusters;

7: Cluster $X$ data
$\displaystyle[X,\text{scope}]=Kmeans\left(X\\_train,n_{c}\right);$

8:else Quantile

9: Enter $n_{q}$ levels;

10: Divide $y$ data
$\displaystyle[y,\text{scope}]=Quantile\left(y\\_score,n_{q}\right);$

11:end if

12:Estimate $covX,meanX,meanY$ $\displaystyle
covX={\big{(}X.group(\text{by=['scope']})\big{)}.\text{cov}()},$
$\displaystyle
meanX={\big{(}X.group(\text{by=['scope']})\big{)}.\text{mean}()},$
$\displaystyle
meanY={\big{(}y.group(\text{by=['scope']})\big{)}.\text{mean}()};$

13:Use Simulated Annealing algorithm  $\displaystyle
w_{0}=SA(meanX,meanY,covX,X\\_train,y\\_train);$

14:return $w_{opt}=w_{0},\alpha_{opt}=\alpha_{0},\beta_{opt}=\beta_{0}$.

In the training process for the best optimal cluster, we set the maximum
cluster to less than half the entire observations, then conduct training from
cluster 1 to the max cluster. The best solution will be preserved with the
corresponding $\alpha,\beta,$ and the cluster if the accuracy score on the
validation set is greater than the original model is. The training process was
illustrated in Figure 2.

##  4 Results

###  4.1 Data

We used the Heart Failure Clinical Records dataset, and the Rice Osmancik
Cammeo dataset to perform our experiments. UCI Machine Learning Repository
site was used to get all datasets.

Figure 2: The training process.

Heart Failure Clinical Records

In April and December of 2015, the Faisalabad Institute of Cardiology and the
Allied Hospital gathered medical records of 299 heart failure patients who
were treated here. The patients, who ranged in age from 40 to 95 years old,
were divided into two groups: 105 females and 194 males, respectively. They
were all classified as classes III or IV of the New York Heart Association
categorization of the stages of heart failure because they had left
ventricular systolic dysfunction and prior heart failures. The dataset
contains 13 attributes including Age, Anaemia, Creatinine phosphokinase,
Diabetes, Ejection fraction, High blood pressure, Platelets, Serum creatinine,
Serum sodium, Sex, Smoking, Time, and the prediction variable Death event.

Rice Osmancik Cammeo

For this research, the Osmancik and Cammeo rice varieties, which have been
farmed in Turkey since 1997 and 2014, respectively, have been chosen. Osmancik
species have broad, lengthy, glassy, and drab appearances when seen as a
whole. The typical features of the Cammeo species are broad and long, glassy,
and dull in appearance. Photos of rice grains were collected from the two
species, which were then processed to get feature inferences from 3810 images.
Each rice grain had its morphology analyzed for seven distinct morphological
traits.

###  4.2 Main results

In the Heart Failure Clinical Records dataset, we used 70 percent of patients
to train a stochastic logistic regression, and 30 percent of the remainder to
verify it on the complete dataset. Furthermore, we compare our model’s
performance with results, which were mentioned by [23] including metrics such
as F1 score, Accuracy score, True positive rate, True negative rate,
Precision-Recall curve, and Roc-AUC score. The formulas are shown in Table 2.
Figure 3, 4, 5, LABEL:fig:re perform the change of Accuracy score, F1 score,
Precision score, and Recall score over the number of clusters (level
quantiles) with the K-Means clustering algorithm, and Q-Quantiles algorithm
correspondingly. In the 10th and 25th clusters, the accuracy score of our
proposed model by both K-Means clustering and Q-Quantiles is outstanding
compared to the original model. Besides, in the same above clusters, the F1
score, Precision score, and recall score are all greater than the metrics of
the original model.

Figure 3: Accuracy score on Heart Failure Clinical Records dataset.  Figure 4:
F1 score on Heart Failure Clinical Records dataset. Figure 5: Precision score
on Heart Failure Clinical Records dataset. Figure 6: Recall score on Heart
Failure Clinical Records dataset.

With the Rice Osmancik Cammeo dataset, we used 75 percent to train a
stochastic logistic regression, and 25 percent to validate the dataset. In
addition, we compare our model’s performance with results, which were
mentioned by [24] including metrics in Table 1. The results are shown in Table
3, we trained the follow-up of 4 executions and calculated the mean of these
evaluation metrics.

Table 4 and Table 5 show the best of the evaluation metrics in $\alpha,\beta$,
and the clusters (level quantiles) correspondingly, these metrics were greater
and improved the ability prediction than the normal method was.

Table 1: Evaluation Metrics were based on TP: True Positive, TN: True negative, FP: False positive, and FN: False negative. Metric | Formula | Metric | Formula
---|---|---|---
Accuracy | $\dfrac{TP+TN}{TP+FP+TN+FN}$ | False Positive | $\dfrac{FP}{TN+FP}$
Sensivitiy (TPR) | $\dfrac{TP}{TP+FN}$ | False Discovery | $\dfrac{FP}{TP+FP}$
Specificity (TNR) | $\dfrac{TN}{TN+FP}$ | False Negative | $\dfrac{FN}{TP+FN}$
Precision | $\dfrac{TP}{TP+FP}$ | MCC | $\dfrac{TP\times TN-FP\times FN}{\sqrt{(TP+FP)(TP+FN)(TN+FP)(TN+FN)}}$
F1 score | $\dfrac{2TP}{2TP+FP+FN}$ | PR AUC | Precision-recall under the curve
Negative Predictive |  |  |
Value | $\dfrac{TN}{TN+FN}$ | ROC AUC | Receiver operating characteristic under the curve
Table 2: 2 top rows: Mean follow-up time of 100 executions for survival prediction results by [23]. 2 bottom rows: Top 2 models which have the best performance on all of the metrics follow-up time of 100 executions. The prediction results of our proposed model using the K-means clustering algorithm ($n_{c}=[1,30]$). Model | MCC | F1 score | Accuracy | TPR | TNR | PR AUC | ROC AUC
---|---|---|---|---|---|---|---
Logistic regression |  $+$0.616 | 0.719 | 0.838 | 0.785 | 0.860 | 0.617 | 0.822
(3 selected features) |  |  |  |  |  |  |
Logistic regression |  $+$0.607 | 0.714 | 0.833 | 0.780 | 0.856 | 0.612 | 0.818
(all features) |  |  |  |  |  |  |
SLR (ours) | $+$0.685 | 0.790 | 0.850 | 0.785 | 0.891 | 0.709 | 0.834
(K-Means algorithm, all features) |  |  |  |  |  |  |
SLR (ours) |  $+$0.677 | 0.793 | 0.850 | 0.809 | 0.872 | 0.699 | 0.841
(Q-Quantiles algorithm, all features) |  |  |  |  |  |  |
Table 3: 7 top rows: Mean follow-up time of 4 executions for survival prediction results by [24]. 2 bottom rows: Top 2 models which have the best performance on all of the metrics follow-up time of 4 executions. The prediction results of our proposed model using the K-means clustering (Q-Quantiles) algorithm ($n_{c}(n_{q})=[1,30]$). Model | Accuracy | Sensivitiy | Specificity | Precision | F1 score | NPV | FPR | FDR | FNR
---|---|---|---|---|---|---|---|---|---
LR | 0.930 | 0.923 | 0.936 | 0.914 | 0.918 | 0.942 | 0.064 | 0.087 | 0.077
MLP | 0.929 | 0.922 | 0.934 | 0.910 | 0.916 | 0.942 | 0.066 | 0.090 | 0.078
SVM | 0.928 | 0.917 | 0.937 | 0.915 | 0.916 | 0.938 | 0.063 | 0.085 | 0.083
DT | 0.925 | 0.912 | 0.935 | 0.913 | 0.912 | 0.934 | 0.065 | 0.087 | 0.088
RF | 0.924 | 0.914 | 0.932 | 0.908 | 0.911 | 0.936 | 0.069 | 0.092 | 0.086
NB | 0.917 | 0.909 | 0.923 | 0.896 | 0.902 | 0.933 | 0.077 | 0.104 | 0.091
k-NN | 0.886 | 0.864 | 0.903 | 0.871 | 0.867 | 0.897 | 0.097 | 0.129 | 0.136
SLR-K (ours) | 0.934 | 0.954 | 0.908 | 0.918 | 0.944 | 0.934 | 0.093 | 0.065 | 0.046
SLR-Q (ours) | 0.932 | 0.963 | 0.889 | 0.910 | 0.942 | 0.946 | 0.111 | 0.077 | 0.037
Table 4: Results of our proposed model on Heart Failure Clinical Records validation data. Model | n-clusters | $\alpha$ | $\beta$ | Accuracy | F1 score | Precision score | Recall score
---|---|---|---|---|---|---|---
| (level-quantiles) |  |  |  |  |  |
SLR (ours) | 10 | 4.5 | 0.855 | 0.856 | 0.723 | 0.810 | 0.692
(K-Means algorithm) |  |  |  |  |  |  |
SLR (ours) | 25 | 4.5 | 0.855 | 0.844 | 0.720 | 0.750 | 0.654
(Q-Quantiles algorithm) |  |  |  |  |  |  |
LR | - | - | - | 0.833 | 0.681 | 0.762 | 0.615
Table 5: Results of our proposed model on Rice Osmancik Cammeo validation data. Model | n-clusters | $\alpha$ | $\beta$ | Accuracy | F1 score | Precision score | Recall score
---|---|---|---|---|---|---|---
| (level-quantiles) |  |  |  |  |  |
SLR (ours) | 12 | 0.3 | 0.62 | 0.938 | 0.947 | 0.944 | 0.969
(K-Means algorithm) |  |  |  |  |  |  |
SLR (ours) | 31 | 4.5 | 0.855 | 0.937 | 0.946 | 0.935 | 0.960
(Q-Quantiles algorithm) |  |  |  |  |  |  |
LR | - | - | - | 0.933 | 0.942 | 0.937 | 0.948

##  5 Discussion

Although numerous numerical approaches have been employed to solve the
Original Logistic Regression issue, they were unable to eliminate the
stochastic occurrences of data in practice. We may entirely explore a
technique to change the chance constraints into the determined constraints by
using the previously described chance constraints programming. The
experiment’s research results show that the effectiveness of the proposed
model yields based on the evaluation metrics with the corporeal dataset.

We overcome the stochastic distribution of data that impacted the inability to
occur with the acquired findings. Despite the extensive computations in the
estimate procedure for the mean vectors or the covariance-variance matrix, the
suggested approach outperforms contemporary machine learning algorithms.

To date, no techniques have addressed the solution of the logistic model by
addressing the programming issue; instead, those studies have shown only a
physical iterative formula that yields a proximal solution. With our novel
concept, we were able to solve many of the different regression situations
without knowing their solution formula.

##  6 Conclusion and Future Work

In this paper, we investigated the incorporation of stochastic elements into
logistic regression through our proposed stochastic logistic regression model.
To address the stochastic nature of the data, we employed the chance
constraint introduced by Charnes and Cooper, treating each data point as a
normal random variable. The key parameters governing the handling of the
chance constraint are denoted as $\alpha$ and $\beta$, representing the
acceptable range and probability. These parameters significantly influence the
outcomes, as evident in the presented graphs. Despite some known asymptotic
properties, selecting optimal values for these parameters remains a
challenging task, especially considering the uniqueness of each dataset,
leaving it as an open problem.

The importance of data scaling cannot be overstated when dealing with large-
scale datasets, offering advantages such as computational speedup. Stochastic
elements in the data are prevalent in real-world scenarios, either inherently
or emerging post-scaling, as we treat groups of data points as new entities
using methods like K-Means clustering or the Q-Quantiles algorithm discussed
earlier.

The capability to navigate uncertainties and randomness provides a crucial
advantage over deterministic methods. Through meticulous analysis and
adjustments, the outcomes can be enhanced, making the approach applicable to a
broader array of problems. We intend to enhance the proposed algorithm for
solving the multiobjective variant of this problem in future endeavors.
Numerous previous investigations have studied this multiobjective problem,
such as optimizing a bicriteria convex programming problem’s efficient set
[25] and implementing an outcome-based branch and bound algorithm [26].
Furthermore, substantial contributions have been made to developing an outcome
space algorithm for generalized multiplicative problems and optimization over
the efficient set [27]. Other related works include the research on optimizing
over the Pareto set of semistrictly quasiconcave vector maximization and its
application in stochastic portfolio selection [28], as well as the framework
for controllable Pareto front learning with completed scalarization functions
and its applications [29].

## References

  * [1] Christopher Bishop.  Pattern Recognition and Machine Learning (Information Science and Statistics).  10 2007.
  * [2] Daniel Jurafsky and James Martin.  Speech and Language Processing: An Introduction to Natural Language Processing, Computational Linguistics, and Speech Recognition, volume 2.  02 2008.
  * [3] Emmanuel Gbenga Dada, Joseph Stephen Bassi, Haruna Chiroma, Shafi’i Muhammad Abdulhamid, Adebayo Olusola Adetunmbi, and Opeyemi Emmanuel Ajibuwa.  Machine learning for email spam filtering: review, approaches and open research problems.  Heliyon, 5(6):e01802, 2019.
  * [4] Mehrbod Sharifi, Eugene Fink, and Jaime Carbonell.  Detection of internet scam using logistic regression.  pages 2168–2172, 10 2011.
  * [5] Nguyen Xuan Tung, Nguyen Quang Dat, Tran Ngoc Thang, Vijender Kumar Solanki, and Nguyen Thi Ngoc Anh.  Analysis of temperature-sensitive on short-term electricity load forecasting.  In 2020 IEEE-HYDCON, pages 1–7, 2020.
  * [6] Miguel Sousa Lobo, Lieven Vandenberghe, Stephen Boyd, and Hervé Lebret.  Applications of second-order cone programming.  Linear Algebra and its Applications, 284(1):193–228, 1998.  International Linear Algebra Society (ILAS) Symposium on Fast Algorithms for Control, Signals and Image Processing.
  * [7] A. Charnes and W. W. Cooper.  Chance-constrained programming.  Manage. Sci., 6(1):73–79, oct 1959.
  * [8] A. Charnes and W. W. Cooper.  Deterministic equivalents for optimizing and satisficing under chance constraints.  Operations Research, 11(1):18–39, 1963.
  * [9] Tran Anh Tuan, Tran Ngoc Thang, Vu Viet Hoang, Do Manh Dung, and Nguyen Thi Ngoc Anh.  A stochastic logistic sigmoid regression using convex programming and clustering.  In 2021 International Conference on Technologies and Applications of Artificial Intelligence (TAAI), pages 205–210, 2021.
  * [10] Werner Fenchel.  Convex cones, sets, and functions.  1953\.
  * [11] M. Avriel, W.E. Diewert, S. Schaible, and I. Zang.  Generalized Concavity.  Classics in Applied Mathematics. Society for Industrial and Applied Mathematics, 1988.
  * [12] Luc Dinh The.  Generalized Convexity in Vector Optimization, pages 195–236.  Springer New York, New York, NY, 2005.
  * [13] Richard W. Cottle.  Letter to the editor - on the convexity of quadratic forms over convex sets.  Oper. Res., 15:170–172, 1967.
  * [14] Olvi L. Mangasarian.  6\. Differentiable Convex and Concave Functions, pages 83–91.  1969\.
  * [15] Stephen P. Boyd and Lieven Vandenberghe.  Convex optimization.  IEEE Transactions on Automatic Control, 51:1859–1859, 2006.
  * [16] N.H. Timm.  Applied multivariate analysis.  Springer texts in statistics. Springer, 2002.
  * [17] M.B. Eisen, P.T Spellman, P.O Brown, and D Botstein.  Cluster analysis and display of genome-wide expression patterns.  Proceedings of the National Academy of Sciences - PNAS, 95(25):14863–14868, 1998.
  * [18] Fazli Can and Esen Ozkarahan.  A dynamic cluster maintenance system for information retrieval.  pages 123–131, 01 1987.
  * [19] Nataliya Boyko, Hanna Komarnytska, Yurii Kryvenchuk, and Yuriy Malynovskyy.  Clustering algorithms for economic and psychological analysis of human behavior.  In Solomiia Fedushko, Sergiy Gnatyuk, Andriy Peleshchyshyn, Zhengbing Hu, Roman Odarchenko, and Igor Korobiichuk, editors, Proceedings of the International Workshop on Conflict Management in Global Information Networks (CMiGIN 2019) co-located with 1st International Conference on Cyber Hygiene and Conflict Management in Global Information Networks (CyberConf 2019), Lviv, Ukraine, November 29, 2019, volume 2588 of CEUR Workshop Proceedings, pages 614–626. CEUR-WS.org, 2019.
  * [20] Trupti Kodinariya and Prashant Makwana.  Review on determining of cluster in k-means clustering.  International Journal of Advance Research in Computer Science and Management Studies, 1:90–95, 01 2013.
  * [21] Ichiro Takeuchi, Quoc V. Le, Timothy D. Sears, and Alexander J. Smola.  Nonparametric quantile estimation.  Journal of Machine Learning Research, 7(45):1231–1264, 2006.
  * [22] Zhiwei Chen and Aoqian Zhang.  A survey of approximate quantile computation on large-scale data.  IEEE Access, 8:34585–34597, 2020.
  * [23] Davide Chicco and Giuseppe Jurman.  Machine learning can predict survival of patients with heart failure from serum creatinine and ejection fraction alone.  BMC Medical Informatics and Decision Making, 20, 02 2020.
  * [24] Ilkay Cinar and Murat Koklu.  Classification of rice varieties using artificial intelligence methods.  International Journal of Intelligent Systems and Applications in Engineering, 7:188–194, 09 2019.
  * [25] Nguyen Thi Bach Kim and Tran Ngoc Thang.  Optimization over the efficient set of a bicriteria convex programming problem.  Pac. J. Optim., 9(1):103–115, 2013.
  * [26] Tran Ngoc Thang.  Outcome-based branch and bound algorithm for optimization over the efficient set and its application.  In Some Current Advanced Researches on Information and Computer Science in Vietnam: Post-proceedings of The First NAFOSTED Conference on Information and Computer Science, pages 31–47. Springer, 2015.
  * [27] Tran Ngoc Thang and Nguyen Thi Bach Kim.  Outcome space algorithm for generalized multiplicative problems and optimization over the efficient set.  2016\.
  * [28] Nguyen Duc Vuong and Tran Ngoc Thang.  Optimizing over pareto set of semistrictly quasiconcave vector maximization and application to stochastic portfolio selection.  Journal of Industrial and Management Optimization, 19(3):1999–2019, 2023.
  * [29] Tran Anh Tuan, Long P Hoang, Dung D Le, and Tran Ngoc Thang.  A framework for controllable pareto front learning with completed scalarization functions and its applications.  Neural Networks, 169:257–273, 2024.


††thanks: A.M. and M.J.K contributed equally to this work.††thanks: A.M. and
M.J.K contributed equally to this work.

# Experimental Demonstration of Broadband Reconfigurable Mechanical
Nonreciprocity

Amin Mehrvarz  Department of Mechanical and Industrial Engineering,
Northeastern University, Boston, MA, 02115, USA.  Mohammad Javad Khodaei
Department of Mechanical and Industrial Engineering, Northeastern University,
Boston, MA, 02115, USA.  Amir Darabi<EMAIL_ADDRESS>Woodruff School of
Mechanical Engineering, Georgia Institute of Technology, Atlanta, GA, 30332,
USA.  Ahmad Zareei<EMAIL_ADDRESS>Harvard John A. Paulson School of
Engineering and Applied Sciences, Harvard University, Cambridge, MA, 02138,
USA.  Nader Jalili<EMAIL_ADDRESS>Department of Mechanical Engineering,
University of Alabama, Tuscaloosa, AL, 35487, USA.  Department of Mechanical
and Industrial Engineering, Northeastern University, Boston, MA, 02115, USA.

###### Abstract

Breaking reciprocity has recently gained significant attention due to its
broad range of applications in engineering systems. Here, we introduce the
first experimental demonstration of a broadband mechanical beam waveguide,
which can be reconfigured to represent wave nonreciprocity. This is achieved
by using spatiotemporal stiffness modulation with piezoelectric patches in a
closed-loop controller. Using a combination of analytical methods, numerical
simulations, and experimental measurements, we show that contrary to the
conventional shunted piezoelectrics or nonlinearity based methods, our setup
is stable, less complicated, reconfigurable, and precise over a broad range of
frequencies. Our reconfigurable nonreciprocal system has potential
applications in phononic logic, wave diodes, energy trapping, and
localization.

Reciprocity is a fundamental property of various physical systems, where the
transmission of a physical quantity, such as waves, between two points in
space is symmetrical. Breaking this reciprocity offers an enhanced control
over wave signal transmission and has recently become of interest in many
branches of physics such as optics Miri  _et al._ (2017); Sounas and Alù
(2017), electronics Dobson (1995), thermodynamics Torrent  _et al._ (2018);
Nakai and Nagaosa (2019), electromagnetism Mahmoud  _et al._ (2015); Caloz
_et al._ (2018), acoustics Nassar  _et al._ (2020); Fleury  _et al._ (2014),
and classical mechanics Sugino  _et al._ (2020); Marconi  _et al._ (2020). In
mechanical systems, breaking wave reciprocity has numerous applications in
trapping waves for efficient energy harvesting devices and designing
transistors in mechanical logic circuits Cummer  _et al._ (2016); Cullen
(1958); Hadad  _et al._ (2016); Felsen and Whitman (1970); Auld  _et al._
(1968).

Breaking reciprocity in mechanical systems can be achieved either by (i)
passively employing nonlinear elements in an asymmetric structure or (ii)
actively varying material properties in space and time periodically Casimir
(1945). Passive methods of breaking reciprocity require high wave amplitudes,
and as a result are impractical in compact devices Sugino  _et al._ (2020).
Additionally, since passively designed structures cannot be reconfigured or
reprogrammed due to their static design, creating robust devices for broad
frequency ranges arises additional complexities Darabi  _et al._ (2019); Fronk
_et al._ (2019); Wu  _et al._ (2018); Boechler  _et al._ (2011). Active
metamaterials, on the other hand, are reprogrammable and tunable by leveraging
the active spatiotemporal modulations. Such active metamaterials offer an
effective platform for breaking reciprocity Riva  _et al._ (2019); Trainiti
and Ruzzene (2016); Nassar  _et al._ (2017), and as a result, have found many
applications in (i) increasing the width of bandgaps Airoldi and Ruzzene
(2011), (ii) focusing or redirecting wave propagation Celli  _et al._ (2017);
Darabi  _et al._ (2018a); Zareei  _et al._ (2018), (iii) changing the
amplitude and phase of transmitted and reflected waves Chen  _et al._ (2018),
and (iv) one-way wave blocking and cloaking Ning  _et al._ (2019); Darabi  _et
al._ (2018b).

In active metamaterials, active elements such as magnetoelastic Chen  _et al._
(2016); Korivand  _et al._ (2021), photosensitive Yannopoulos and Trunov
(2009), or piezoelectric Airoldi and Ruzzene (2011); Mehrvarz  _et al._
(2019), are used to modulate physical properties. In the photosensitive and
magnetoelastic materials, the stiffness can be varied by changing the magnetic
field Chen  _et al._ (2016); Korivand  _et al._ (2021) and temperature
Yannopoulos and Trunov (2009), respectively. In piezoelectric transducers,
however, the stiffness is modulated by connecting these patches to shunted
circuits Airoldi and Ruzzene (2011); Marconi  _et al._ (2020) with negative
capacitance Tateo  _et al._ (2015). Although shunted piezoelectric patches are
precise Darabi  _et al._ (2020); Sugino  _et al._ (2020); Trainiti  _et al._
(2019); Chen  _et al._ (2019) and able to function in a wide frequency range,
a dramatic change in the equivalent Young’s modulus only happens when the
system operates very close to the unstable zones of the circuit Airoldi and
Ruzzene (2011). As such, the shunted piezoelectric patches are prone to an
error where a small variation of the applied negative capacitance can make the
system unstable and result in a large deviation from desired values Tateo  _et
al._ (2015). Additionally, implementing a system with shunted piezoelectric
patches is extremely challenging since each piezoelectric element is
separately connected to its own independent circuit and this requires working
with too many elements and connections.

Here, we present a closed-loop feedback control system connected to two sets
of parallel piezoelectric patches 111more specifically, the PZT-5J type, for
convenience, PZT refers to this material in the rest of this manuscript bonded
on a host beam to modulate the beam’s stiffness. The controller here
continually measures the voltage from one set of PZTs and applies the required
analog signals to the other set, and as such, changes the beam’s effective
stiffness Li  _et al._ (2017); Sugino  _et al._ (2018). Contrary to shunted
PZTs with negative capacitance circuits, a closed-loop stiffness modulation
system is stable and precise with reduced complexity and has functionality
over a broad range of frequencies. As a result, stiffness modulation with a
closed-loop control system offers a robust platform for breaking reciprocity.
Here, we first analytically/numerically show that spatiotemporal (or spatial)
stiffness modulation in beam results in directional (or total) bandgaps. We
then experimentally exhibit the bandgaps using a closed-loop controlled
system. We show that different behaviors (i.e., wave transmission, blockage,
and nonreciprocal transmission) are all attainable on the same reconfigurable
beam setup by only tuning the controller’s parameters.

Figure 1: Nonreciprocal wave transmission: (a) Schematic of a beam structure
with a spatiotemporally modulated section in the middle. The spatiotemporal
modulation results in a nonreciprocal wave propagation allowing wave to
propagate from right to left and blocking waves from left to right. (b) The
stiffness modulation diagram for the structure is shown in space and time. The
spatial/temporal modulation wavelength/period are $\lambda_{m}$ and $T_{m}$.

## Results

Band-diagrams of a spatio-temporally modulated beam. We consider elastic wave
propagation in a thin beam with linear mass density $\rho_{l}$, rigidity
$E_{0}$, second moment of inertial $I_{0}$, and stiffness $D_{0}=E_{0}I_{0}$
where the displacement $w$ is governed by the Euler-Bernoulli equation

$\rho_{l}\frac{{{\partial}^{2}}w}{\partial{{t}^{2}}}+\frac{{{\partial}^{2}}}{\partial{{x}^{2}}}\left(D(x,t)\frac{{{\partial}^{2}}w(x,t)}{\partial{{x}^{2}}}\right)=0.$
(1)

Assuming a spatiotemporal modulation in stiffness with spatial modulation
wavelength $\lambda_{m}$ and temporal modulation period $T_{m}$, the stiffness
is obtained as
$D(x,t)=D_{0}\left(1+\alpha_{m}\cos\left(\omega_{m}t-k_{m}x\right)\right)$,
where $k_{m}=2\pi/\lambda_{m}$ is the spatial modulation wavenumber,
$\omega_{m}=2\pi/T_{m}$ is the temporal modulation angular frequency, and
$\alpha_{m}$ is the modulation amplitude (Fig. 1b). Inserting the modulation
stiffness, $D(x,t)$ into Eq. (1) and taking the Fourier transform, we find the
characteristic equation

$\displaystyle\begin{split}&\frac{{{\alpha}_{m}}\gamma}{2}{{\left(k+\left(n-1\right){{k}_{m}}\right)}^{2}}{{\left(k+n{{k}_{m}}\right)}^{2}}{{{\hat{w}}}_{n-1}}\\\
&+\frac{{{\alpha}_{m}}\gamma}{2}{{\left(k+\left(n+1\right){{k}_{m}}\right)}^{2}}{{\left(k+n{{k}_{m}}\right)}^{2}}{{{\hat{w}}}_{n+1}}\\\
&+\left[\gamma{{\left(k+n{{k}_{m}}\right)}^{4}}-{{\left(\omega+n{{\omega}_{m}}\right)}^{2}}\right]{{{\hat{w}}}_{n}}=0,\end{split}$
(2)

where $\gamma=D_{0}/\rho_{l}$, and $\hat{w}_{n}$ is the $n^{\text{th}}$ mode
of the wave amplitude Fourier transform given by
$w(x,t)={{e}^{\mathrm{i}\left(\omega
t-kx\right)}}\sum\limits_{p=-\infty}^{+\infty}{{{{\hat{w}}}_{n}}{{e}^{\mathrm{i}p\left({{\omega}_{m}}t-{{k}_{m}}x\right)}}}$.
The band diagram of a spatiotemporally modulated beam is the solution of Eq.
(2) and can be found for different modulation parameters $k_{m}$,
$\omega_{m}$, and $\alpha_{m}$. Achieving different behaviors is possible by
selecting proper modulation functions and simply changing the controller’s
parameters accordingly (see Supplementary Information for the effects of
density and stiffness modulation with or without a phase difference between
them.) To avoid nonlinear effects, we assume a small fixed value of
$\alpha_{m}=0.15$ in our modulation. Initially assuming only a spatial
modulation, i.e., $k_{m}\neq 0,\omega_{m}=0$, the band diagram exhibits a full
bandgap as shown in Fig. 2a. Next, by introducing and varying temporal
modulation, i.e., $\omega_{m}\neq 0$, we find that the bandgap becomes
asymmetrical (Fig. 2b). Furthermore, increasing the temporal modulation
results in a fully asymmetrical bandgap, which has a nonreciprocal wave
behavior with unidirectional wave propagation (Fig. 2c). In the case of only
temporal modulation, i.e., $k_{m}=0,\omega_{m}\neq 0$, we find that the
bandgaps are shifted to the wavenumbers where certain wavelengths are not
allowed in the system (see Fig. 2d) Trainiti  _et al._ (2019). (A continuous
changing of the modulation parameters (i.e., $k_{m}$ and $\omega_{m}$) can
provide a better insight into spatiotemporal modulation, shown in
supplementary movies SV1 and SV2.)

Figure 2: Band diagrams of an aluminium beam ($E=69.9$ GPa and $\rho=2700$
kg.m-3) with $1$ mm thickness and $2.5$ cm width. (a) Complete bandgap for
only space modulated beam $\alpha_{m}=0.15$, $k_{m}=80$ m-1, and
$\omega_{m}=0$. (b) Directional and complete bandgaps in the beam with
$\alpha_{m}=0.15$, $k_{m}=80$ m-1, and $\omega_{m}=400$ Hz. (c) Directional
bandgaps in the beam with $\alpha_{m}=0.15$, $k_{m}=80$ m-1, and
$\omega_{m}=800$ Hz. (d) Wavenumber bandgap with $\alpha_{m}=0.15$, $k_{m}=0$,
and $\omega_{m}=800$ Hz (for more information on the effect of $k_{m}$ and
$\omega_{m}$ on band diagram and bandgaps, readers are referred to
supplementary videos A and B).

Stiffness modulation with PZT patches and closed-loop control. We use a
closed-loop circuit with two parallel PZT patches (one actuator and one
sensor) on both sides of the beam to implement spatiotemporal modulation in
stiffness. Based on the Euler-Bernoulli theory, the stiffness in a beam is
defined $M=-Dw_{xx}$, where $M$ is the moment, and $w_{xx}$ is the curvature
of the beam. When a PZT actuator and a PZT sensor are attached on each side of
a beam extending between $x_{l}$ and $x_{r}$, the stiffness relation is
modified to Jalili (2010)

$D=\frac{\int_{{{x}_{l}}}^{{{x}_{r}}}{M}dx}{-\int_{x_{l}}^{x_{r}}w_{xx}dx}=D_{0}+D_{p}-\frac{K_{p}V_{a}}{\int_{x_{l}}^{x_{r}}w_{xx}dx},$
(3)

where $D_{p}$ is the stiffness of the PZT patches, $V_{a}$ is the PZT actuator
voltage, and $K_{p}$ is a factor that depends on PZT parameters (see
Supplementary Information for more details). If the actuator’s voltage,
$V_{a}$, becomes

$V_{a}=\frac{\left(D_{p}-\alpha
D_{0}\right)}{K_{p}}\int_{x_{l}}^{x_{r}}w_{xx}dx,$ (4)

then the beam’s rigidity changes to $D=D_{0}(1+\alpha)$. Note that $\alpha$ is
the modulation amplitude and can depend on the location of PZT, $x$, and also
time $t$. To apply the actuator’s voltage in Eq. (4), the only unknown
parameter on the right-hand side is the change in the beam’s slope, i.e.,
$\int_{x_{l}}^{x_{r}}w_{xx}dx$, where it can be obtained based on the PZT
sensor and its dielectric permittivity Yi  _et al._ (2019). In a PZT sensor
$V_{s}=C_{p}\int_{x_{l}}^{x_{r}}w_{xx}dx$, where $V_{s}$ is the sensor’s
voltage, and $C_{p}$ is a factor that depends on PZT parameters (see
Supplementary Information for the detail). As a result, given sensor voltage
$V_{s}$, if the actuator voltage is applied as

$V_{a}=\frac{\left(D_{p}-\alpha D_{0}\right)}{K_{p}C_{p}}V_{s},$ (5)

the flexural rigidity of the beam becomes $D=D_{0}(1+\alpha)$. Modulating the
coefficient $\alpha$ using $\alpha=\alpha_{m}\cos(\omega_{m}t+\phi)$ for each
PZT set (sensor and actuator) and applying the phase $\phi$ based on the
location of the PZT set, a spatiotemporal modulation can be achieved.

Figure 3: (a) The schematic of the spatiotemporally modulated beam with PZT
actuators (blue) and sensors (red) bonded on both sides of the beam. (b) The
stiffness modulation coefficient $\alpha_{m}$ of three consecutive PZT patches
over time. Note that the phase difference between consequent patches results
in the spatial modulation. (c) The schematic of the closed-loop feedback
controller to implement the active modulation. The PZT sensor’s voltage is
used in Eq. (5) to obtain actuating voltage, which is applied to the actuator
PZT. (d) The numerical transfer function (ratio between transmitted wave to
incident wave) without any modulation (black-line) and with spatial modulation
with $k_{m}=80.5$ m-1 for waves propagating left-to-right (red line) and waves
propagating right-to-left (blue line). Since there is no time modulation, the
results are identical for both LR and RL propagations. The yellow boxes
indicate the bandgaps. (e) Numerically calculated transmission ratio for
spatiotemporally modulated beam with $f_{m}=400$ Hz and $k_{m}=80.5$ m-1 for
waves propagating RL (blue line) and LR (red line). The red boxes indicate the
nonreciprocal bandgaps.

Numerical simulations of spatially and spatio-temporally modulated beam. We
first numerically test the effect of spatiotemporal modulation of the PZTs on
an aluminum beam using Finite Element simulations in COMSOL Multiphysics. We
consider an aluminum beam with the thickness of $1$ mm, width of $2.54$ mm,
and length of $1.5$ m. We further assume twelve pairs of PZT patches perfectly
bonded on both sides of the beam. Each PZT covers an area of $21$ mm$\times
21$ mm, thickness $0.55$ mm, and we have a distance of $5$ mm in-between them
(see Fig. 3a). We choose these numbers based on an experimentally feasible
setup. A low reflection boundary is used on both sides of the beam to avoid
reflecting the wave. The closed-loop parameters are selected so that the
stiffness of the piezoelectric actuators vary by $15$ percent, i.e.,
$\alpha_{m}=0.15$. The closed-loop feedback controller based on Eq. (5) is
used to modulate the PZT (see Fig. 3c). The phase difference between
subsequent PZTs is set to $\phi=2\pi/3$, which results in three PZT pairs per
spatial wavelength $\lambda_{m}$ (Fig. 3b). It is to be noted that the
equivalent stiffness of each PZT pair, contrary to the shunted circuits, can
vary continuously based on the continuous modulation signal (See Fig. 3b). To
test the setup, the wave is initiated on one side of the structure, and
measurements are done on both sides of the spatiotemporally modulated section.
Finally, a sweep over a range wave frequencies is done to obtain the transfer
function (i) in the absence of any modulation, or in the presence of (ii)
spatial modulation (i.e., $\omega_{m}=0,k_{m}\neq 0$), or (iii) spatiotemporal
modulation (i.e., $\omega_{m}\neq 0,k_{m}\neq 0$).

In Fig. 3 d-e, we plot the transfer function, i.e., the ratio of transmitted
wave amplitude to incident wave amplitude, in the frequency range
$2-10\leavevmode\nobreak\ $kHz for the spatial/spatiotemporal modulation of
the beam. Here, the transmission ratio below $r_{t}\leq 0.1$ is considered
blockage because the energy level, which is proportional to the deflection
square, drops more than 0.01 Khelif and Adibi (2015). Note that this
transmission level corresponds to $-20dB$ which is commonly considered as
bandgap Yi  _et al._ (2017); Zouari  _et al._ (2018); Tang and Cheng (2017).
First, in the absence of modulation, no bandgap is observed over the frequency
range of interest (black line in Fig. 3d). In the presence of spatial
modulation, however, the transfer function changes and results in identical
transfer functions for left-to-right (LR) and right-to-left (RL) wave of
propagation (see red and blue lines in Fig. 3d respectively). Additionally,
the results show two identical bandgaps for opposite propagation direction (LR
and RL) bounded between $2.5-4.5$ kHz, marked with yellow boxes in Fig. 3d,
and $7-8.5$ kHz. Interestingly, the first bandgap contains a localized mode at
$3.16$ kHz, which dominates the transfer function response, and results in a
spike in the transmission ratio in the bandgap. Next, we run the
spatiotemporal modulated simulations (i.e., $\omega_{m}\neq 0,k_{m}\neq 0$)
and report the results in Fig. 3e for the wave propagation in LR and RL
directions. This figure shows the presence of three directional bandgaps
(marked with red rectangles), in which the magnitude of the wave for LR is
different from RL. As demonstrated, for two of these bandgaps ($2.65-2.86$
kHz, and $6.46-7.48$ kHz), waves only travel for LR, while for the
($7.97-8.69$ kHz) waves propagate in the opposite direction. Furthermore,
bandgaps move to lower frequencies for the RL and higher frequencies for the
LR, respectively. Nonetheless, this change is not visible for the first
bandgap due to the presence of localized mode in this frequency range.

Figure 4: (a) Schematic of the experimental setup showing the beam, absorbing
patches, PZTs sensors and actuators, piezodrives, MFC actuators, controller,
and the laser. MFCs on both sides are used to introduce the wave to the
system. The wave amplitude is measured before/after the first/last PZTs using
a laser vibrometer. (b) Experimental setup. (c) The experimental measurements
of transmission ratio without any modulation (solid black line) and with
spatial modulation with $k_{m}=80.5$ m-1 for wave propagating LR (red) and RL
(blue). Note that the transmission ratio for RL and LR wave are identical. The
bandgaps are indicated with the yellow box. (d) Transmission ratio in presence
of spatiotemporal modulation with $k_{m}=80.5$ m-1 and $f_{m}=400$ Hz for LR
(red) and RL (blue) waves. The nonreciprocal bandgaps are identified with the
red box.

Experimental frequency response measurement for spatially and spatio-
temporally modulated beam. We implement the time-periodic stiffness modulation
of the elastic waveguide beam using an array of PZT actuators and sensors
attached on a flexible aluminum beam controlled using closed-loop circuits
(Fig. 4). The aluminum beam has a rectangular cross-section similar to the
numerical simulation, and layers of Butyl rubber at its boundaries damp the
wave reflection on the structure. The waves are produced using Macro Fiber
Composite (MFC) actuators that are bonded on two ends of the beam, and a laser
vibrometer (Polytec CLV-2544) measures the system’s response at a distance of
$5$ mm before and after the modulated section (Fig. 4a). The voltages of PZT
sensors are measured using a DS1006 R&D controller board. The measured
voltages are used in the closed-loop circuits in Simulink-MATLAB with
reconfigurable parameters of Eq. (5) (see Supplementary Information for more
details). As can be seen in Fig. 4a, ControlDesk software with DS1006 R&D
controller board applies the control voltages to the piezodrive amplifier, and
the outputs are applied to the PZT actuators. See Fig. 4b for the experimental
setup. The experimental results for the transformation ratio without any
modulation are shown in Fig. 4c with a solid black line.

## Discussion

Introducing spatial modulation with $k_{m}=80.5$ m-1, the transmission ratio
changes and results in two identical bandgaps for RL and LR wave propagation
direction at frequencies $2.86-3.33$ kHz and $3.18-3.6$ kHz (shown with yellow
boxes in Fig. 4c). Note that in the analytical results, this bandgap is
continuous from 2.86-3.6 kHz without any break in the middle. However, in the
experiments similar to the numerical simulations, this bandgap breaks into two
separate bandgaps due to an internal localized mode of the system. Introducing
a time modulation in addition to the spatial modulation with $f_{m}=400$ Hz
results in asymmetrical bandgaps that depend on the direction of wave
propagation. Figure 4d shows the presence of four directional bandgaps (red
boxes), in which, same as the numerical simulations, the magnitude of the wave
for LR is different from RL. Interestingly, two of the bandgaps (specifically
$2.84-2.9$ kHz, and $3.16-3.27$ kHz) allow only LR wave propagation (RL
nonreciprocal bandgaps), while on the contrary, the other two band gaps (i.e.,
$3.01-3.13$ kHz and $3.57-3.6$ kHz) allow only RL waves propagation (LR
nonreciprocal bandgaps). Moreover, similar to the numerical results, bandgaps
move to lower/higher frequencies for the RL/LR bandgaps. It should be noted
that two of the directional bandgaps are observed only because the bandgaps
are in the vicinity of the interior localized mode. Additionally, we expect
the experimental result to defer from the numerics due to the following
sources of error, (i) In the numerical simulations, the contacts between the
PZTs and beam are assumed to be perfect contacts; however, there is an
adhesive layer between the PZTs and aluminum beam that changes the effective
thickness of the system but also changes the corresponding density of the
system Rabinovitch and Vinson (2002). Additionally, (ii) our experiments
include structural damping and nonlinear effects, which are ignored in the
simulations and theory. Lastly, (iii) including the equipment’s error (e.g.,
signal generator, piezodrive amplifier, and controller) would further defer
our experimental result from the numerical simulation.

In summary, we experimentally implemented a beam structure to achieve
total/nonreciprocal bandgaps. For the first time, an active closed-loop
feedback control system was employed using an array of PZTs bonded on both
sides of an aluminum host layer to implement the stiffness modulation. At
first, we investigated the changing effects of spatiotemporal parameters on
the bandgaps. Then, we numerically simulated spatial/spatiotemporal stiffness
modulation using the closed-loop system and validated our analytical results.
Finally, we experimentally showed that actively controlled stiffness
modulation using closed-loop circuits is reconfigurable and can be used for
both spatial and spatiotemporal modulation, and has potential applications in
waveguides, diodes, phononic logic circuits, or energy localization.

## Data availability

The data that support the plots within this paper and other findings of this
study are available from the corresponding author upon request.

## References

  * Miri  _et al._ (2017) M.-A. Miri, F. Ruesink, E. Verhagen, and A. Alù, Physical Review Applied 7, 064014 (2017).
  * Sounas and Alù (2017) D. L. Sounas and A. Alù, Nature Photonics 11, 774 (2017).
  * Dobson (1995) I. Dobson, IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications 42, 517 (1995).
  * Torrent  _et al._ (2018) D. Torrent, O. Poncelet, and J.-C. Batsale, Physical review letters 120, 125501 (2018).
  * Nakai and Nagaosa (2019) R. Nakai and N. Nagaosa, Physical Review B 99, 115201 (2019).
  * Mahmoud  _et al._ (2015) A. M. Mahmoud, A. R. Davoyan, and N. Engheta, Nature communications 6, 1 (2015).
  * Caloz  _et al._ (2018) C. Caloz, A. Alù, S. Tretyakov, D. Sounas, K. Achouri, and Z.-L. Deck-Léger, Physical Review Applied 10, 047001 (2018).
  * Nassar  _et al._ (2020) H. Nassar, B. Yousefzadeh, R. Fleury, M. Ruzzene, A. Alù, C. Daraio, A. N. Norris, G. Huang, and M. R. Haberman, Nature Reviews Materials 5, 667 (2020).
  * Fleury  _et al._ (2014) R. Fleury, D. L. Sounas, C. F. Sieck, M. R. Haberman, and A. Alù, Science 343, 516 (2014).
  * Sugino  _et al._ (2020) C. Sugino, M. Ruzzene, and A. Erturk, Physical Review B 102, 014304 (2020).
  * Marconi  _et al._ (2020) J. Marconi, E. Riva, M. Di Ronco, G. Cazzulani, F. Braghin, and M. Ruzzene, Physical Review Applied 13, 031001 (2020).
  * Cummer  _et al._ (2016) S. A. Cummer, J. Christensen, and A. Alù, Nature Reviews Materials 1, 16001 (2016).
  * Cullen (1958) A. Cullen, Nature 181, 332 (1958).
  * Hadad  _et al._ (2016) Y. Hadad, J. C. Soric, and A. Alu, Proceedings of the National Academy of Sciences 113, 3471 (2016).
  * Felsen and Whitman (1970) L. Felsen and G. Whitman, IEEE Transactions on Antennas and Propagation 18, 242 (1970).
  * Auld  _et al._ (1968) B. Auld, J. H. Collins, and H. Zapp, Proceedings of the IEEE 56, 258 (1968).
  * Casimir (1945) H. B. G. Casimir, Reviews of Modern Physics 17, 343 (1945).
  * Darabi  _et al._ (2019) A. Darabi, L. Fang, A. Mojahed, M. D. Fronk, A. F. Vakakis, and M. J. Leamy, Physical Review B 99, 214305 (2019).
  * Fronk  _et al._ (2019) M. D. Fronk, S. Tawfick, C. Daraio, S. Li, A. Vakakis, and M. J. Leamy, Journal of Vibration and Acoustics 141 (2019).
  * Wu  _et al._ (2018) Z. Wu, Y. Zheng, and K. Wang, Physical Review E 97, 022209 (2018).
  * Boechler  _et al._ (2011) N. Boechler, G. Theocharis, and C. Daraio, Nature materials 10, 665 (2011).
  * Riva  _et al._ (2019) E. Riva, J. Marconi, G. Cazzulani, and F. Braghin, Journal of Sound and Vibration 449, 172 (2019).
  * Trainiti and Ruzzene (2016) G. Trainiti and M. Ruzzene, New Journal of Physics 18, 083047 (2016).
  * Nassar  _et al._ (2017) H. Nassar, X. Xu, A. Norris, and G. Huang, Journal of the Mechanics and Physics of Solids 101, 10 (2017).
  * Airoldi and Ruzzene (2011) L. Airoldi and M. Ruzzene, New Journal of Physics 13, 113010 (2011).
  * Celli  _et al._ (2017) P. Celli, S. Gonella, V. Tajeddini, A. Muliana, S. Ahmed, and Z. Ounaies, Smart Materials and Structures 26, 035001 (2017).
  * Darabi  _et al._ (2018a) A. Darabi, A. Zareei, M.-R. Alam, and M. J. Leamy, Scientific reports 8, 1 (2018a).
  * Zareei  _et al._ (2018) A. Zareei, A. Darabi, M. J. Leamy, and M.-R. Alam, Applied Physics Letters 112, 023901 (2018).
  * Chen  _et al._ (2018) Y. Chen, X. Li, H. Nassar, G. Hu, and G. Huang, Smart Materials and Structures 27, 115011 (2018).
  * Ning  _et al._ (2019) L. Ning, Y.-Z. Wang, and Y.-S. Wang, International Journal of Mechanical Sciences 153, 287 (2019).
  * Darabi  _et al._ (2018b) A. Darabi, A. Zareei, M.-R. Alam, and M. J. Leamy, Physical review letters 121, 174301 (2018b).
  * Chen  _et al._ (2016) Y. Chen, J. Hu, and G. Huang, Journal of Intelligent Material Systems and Structures 27, 1337 (2016).
  * Korivand  _et al._ (2021) S. Korivand, A. Mehrvarz, N. Candelino, M. J. Khodaei, and N. Jalili, ASME Letters in Dynamic Systems and Control 1 (2021), 041007.
  * Yannopoulos and Trunov (2009) S. Yannopoulos and M. Trunov, physica status solidi (b) 246, 1773 (2009).
  * Mehrvarz  _et al._ (2019) A. Mehrvarz, A. Najafi Ardekani, M. J. Khodaei, and N. Jalili, Journal of Vibration and Control 25, 1393 (2019).
  * Tateo  _et al._ (2015) F. Tateo, M. Collet, M. Ouisse, M. Ichchou, K. Cunefare, and P. Abbe, Journal of Intelligent Material Systems and Structures 26, 952 (2015).
  * Darabi  _et al._ (2020) A. Darabi, X. Ni, M. Leamy, and A. Alù, Science advances 6, eaba8656 (2020).
  * Trainiti  _et al._ (2019) G. Trainiti, Y. Xia, J. Marconi, G. Cazzulani, A. Erturk, and M. Ruzzene, Physical review letters 122, 124301 (2019).
  * Chen  _et al._ (2019) Y. Chen, X. Li, H. Nassar, A. N. Norris, C. Daraio, and G. Huang, Physical Review Applied 11, 064052 (2019).
  * Note (1) More specifically, the PZT-5J type, for convenience, PZT refers to this material in the rest of this manuscript.
  * Li  _et al._ (2017) F. Li, C. Zhang, and C. Liu, Journal of Sound and Vibration 393, 14 (2017).
  * Sugino  _et al._ (2018) C. Sugino, M. Ruzzene, and A. Erturk, IEEE/ASME Transactions on Mechatronics 23, 2144 (2018).
  * Jalili (2010) N. Jalili, “Piezoelectric-based systems modeling,” in  _Piezoelectric-Based Vibration Control: From Macro to Micro/Nano Scale Systems_ (Springer US, Boston, MA, 2010) pp. 183–232.
  * Yi  _et al._ (2019) K. Yi, M. Ouisse, E. Sadoulet-Reboul, and G. Matten, Smart Materials and Structures 28, 065025 (2019).
  * Khelif and Adibi (2015) A. Khelif and A. Adibi,  _Phononic crystals_ (Springer, 2015).
  * Yi  _et al._ (2017) K. Yi, M. Collet, and S. Karkar, Physical review B 96, 104110 (2017).
  * Zouari  _et al._ (2018) S. Zouari, J. Brocail, and J.-M. Génevaux, Journal of Sound and Vibration 435, 246 (2018).
  * Tang and Cheng (2017) L. Tang and L. Cheng, Journal of Applied Physics 121, 194901 (2017).
  * Rabinovitch and Vinson (2002) O. Rabinovitch and J. R. Vinson, Journal of intelligent material systems and structures 13, 689 (2002).

## Acknowledgments

The authors declare that they have no known competing for financial interests
or personal relationships that could have appeared to inﬂuence the work
reported in this paper.