\section{Introduction}
\label{sec:Introduction}

The celebrated filtered-x least mean square (FxLMS) algorithm and
its variants have been successfully applied in active noise control
(ANC) systems. However, in practical situations, the primary path
$P(z)$ and secondary path $S(z)$ may be nonlinear
\cite{das2004active}. Moreover, the reference noise $d(n)$ arises
from dynamic systems and as such, the noise may be a nonlinear and
deterministic or stochastic, colored, and non-Gaussian signal
\cite{george2013review}. In these cases, linear ANC techniques
cannot fully utilize the coherence in the noise and achieve
suboptimal performance.


\begin{table}[h]
    \scriptsize
    \centering \doublerulesep=0.05pt
    \caption{Time line of NLANC.}
    \begin{tabular}{l|l|l|l}
        \hline\hline
        \textbf{Years} & \textbf{Authors} & \textbf{Contributions} & \textbf{References} \\ \hline
        1994      &Wangler and Hansen         &\begin{tabular}[c]{@{}l@{}}Genetic algorithm for nonlinear\\active noise and vibration control\end{tabular}               & \cite{wangler1994genetic}            \\ \hline
        1995      &Snyder and Tanaka         &ANN for NLANC               &\begin{tabular}[c]{@{}l@{}} \cite{snyder1995active}\end{tabular}            \\ \hline
        1996      &Tang et al.          &\begin{tabular}[c]{@{}l@{}}Genetic algorithm for ANC\end{tabular}                &\cite{tang1996application}            \\ \hline
        1997      &Tokhi and Wood          &\begin{tabular}[c]{@{}l@{}}Radial basis function networks for NLANC\end{tabular}                &\cite{tokhi1997active}            \\ \hline
        1998      &Strauch and Mulgrew          &\begin{tabular}[c]{@{}l@{}}A combined linear and
            nonlinear controller to\\compensate for the nonlinear effects\end{tabular}                & \cite{strauch1998active}            \\ \hline
        1999      &\begin{tabular}[c]{@{}l@{}}Bouchard, Paillard,\\and Dinh\end{tabular}          &\begin{tabular}[c]{@{}l@{}}Five new algorithms for the learning of a\\multilayer perceptron neural network for NLANC\end{tabular}                &\cite{bouchard1999improved}            \\ \hline
        1997, 2001    &Tan and Jiang         &VFxLMS algorithm        & \cite{tan1997filtered,tan2001adaptive}            \\ \hline
        2001    &Bouchard         &\begin{tabular}[c]{@{}l@{}}Filtered-neuron-level extended Kalman\\algorithm (FxNEKA) and filtered-x extended\\Kalman filter (FxEKF) for NLANC\end{tabular}        &\cite{bouchard2001new}            \\ \hline
        2004      &Das and Panda         &FsLMS algorithm               &\cite{das2004active}            \\ \hline
        2004      &Sicuranza and Carini         &\begin{tabular}[c]{@{}l@{}}Volterra filtered-x affine projection\\ (VFxAP) algorithm\end{tabular}                & \cite{sicuranza2004filtered}            \\ \hline
        2005      &Sicuranza and Carini         &\begin{tabular}[c]{@{}l@{}}Implementation of  the VFxAP algorithm\\based on V-vector algebra\end{tabular}     & \cite{sicuranza2005multichannel}            \\ \hline
        2005      &Kuo and Wu         &\begin{tabular}[c]{@{}l@{}}Bilinear FxLMS algorithm\end{tabular}                &\cite{kuo2005nonlinear}           \\ \hline
        2006      &\begin{tabular}[c]{@{}l@{}}Modares, Ahmadyfard,\\ and Hadadzarif \end{tabular}          &\begin{tabular}[c]{@{}l@{}}PSO for NLANC\end{tabular}                & \cite{modares2006pso}            \\ \hline
        2004, 2006      &Zhang et al.          &\begin{tabular}[c]{@{}l@{}}Fuzzy neural networks for NLANC\end{tabular}     & \cite{zhang2004active,zhang2006adaptive}            \\ \hline
        2006      &Sicuranza and Carini         &\begin{tabular}[c]{@{}l@{}}Generalized Hammerstein model for NLANC\end{tabular}     & \cite{sicuranza2006accuracy}            \\ \hline
        2007      &Chang and Luoh          &\begin{tabular}[c]{@{}l@{}}Neural-based FxLMS algorithm\end{tabular}                &\cite{chang2007enhancement}            \\ \hline
        2008      &Bambang         &\begin{tabular}[c]{@{}l@{}}Recurrent neural networks for NLANC\end{tabular}     & \cite{bambang2008adjoint}            \\ \hline
        2008      &\begin{tabular}[c]{@{}l@{}}Reddy, Das, \\and
            Prabhu\end{tabular}          &Fast FsLMS algorithm             &\cite{reddy2008fast} \\ \hline\hline
    \end{tabular}
    \label{Table000}
\end{table}

To address these limitations, several nonlinear ANC (NLANC)
algorithms were proposed, yielding a more stable performance and
lower noise level
\cite{wangler1994genetic,sicuranza2004filtered,reddy2008fast}.  It
should be emphasized that some NLANC methods mentioned in this paper
first appeared before 2009. The research on NLANC dates back to
1995, 59 years after the concept of ANC was proposed by Lueg
\cite{paul1936process}. In this year, Snyder and Tanaka developed
ANN for NLANC problem, and so far a large number of NLANC algorithms
have been proposed. The milestones in the progress of NLANC before
2009 are presented in Table \ref{Table000}. Survey papers on NLANC
techniques have been reported \cite{george2013review}. However,
these surveys do not cover the literature since 2013. To complete
the review of NLANC techniques and include the advances, in this
paper we review the research progress of NLANC models in recent
years.

The Volterra filter is gaining importance in NLANC. According to the
\textit{Stone-Weierstrass approximation theorem}, the Volterra
filter is a `universal approximator', i.e., the filter can
arbitrarily approximate discrete-time with fading memory,
time-invariant, continuous, finite-memory, nonlinear system
\cite{burt2018efficient}. The Volterra filter belongs to
linear-in-the-parameters (LIP) filters, where the output depends
linearly on the parameters of the filter itself. Note that this type
of filter becomes computationally expensive when a large number of
coefficients are required \cite{burt2018efficient}. To solve this
problem, the second-order Volterra (SOV) series and the third-order
Volterra (TOV) series are often used. For example, in
\cite{lu2016adaptive,he2019an}, the Volterra filtered-x-based
algorithms that use of SOV were developed for active impulsive noise
control (AINC) problem.

The NLANC based on an artificial neural network (ANN) method was
first studied, which provides the capacity to maintain causality
within the control scheme. From 1997 to 2008, many papers discussed
the usage of radial basis function (RBF) networks, recurrent neural
networks (RNNs), and fuzzy neural networks for NLANC (see Table
\ref{Table000}). In the last decade, some attempts were made by
using an ANN as the controller
\cite{le2017adaptive,zhang2010adaptive,azadi2012filtered}. A neural
controller equipped with the filtered-u least mean square (FuLMS)
algorithm with correction terms momentum was proposed in
\cite{chang2010neural}. The correction terms momentum can
efficiently prevent the system from unstable poles and can enhance
the performance as compared with the conventional leaky FuLMS
algorithm.  Regrettably, the above mentioned ANNs have heavy
computational complexity which may hinder their practical
applications \cite{george2013review}.

The functional link ANN (FLANN) exploits the single layer of ANN to
obtain reliable estimation performance with acceptable error level
\cite{reddy2009fast}. The FLANN filter is also a type of LIP filter,
which shares a unified implementation with the Volterra filter.
Moreover, the linear adaptive filtering algorithms can easily extend
to this framework with favorable implications on the algorithm
stability and computational complexity \cite{sicuranza2011a}. Over
the past decade, a number of filtered-s least mean square (FsLMS)
algorithms were proposed by employing FLANN, such as fast FsLMS
\cite{reddy2009fast}, robust FsLMS (RFsLMS) \cite{george2012a}. Most
of these algorithms use trigonometric functional expansion as the
basis function, and a small number of algorithms that consider using
Chebyshev \cite{chen2019nonlinear}, Fourier \cite{carini2012new},
and Legendre expansions \cite{george2013active} were also developed.

To compensate for the strong nonlinearity of  the system with
saturated signals, several bilinear algorithms with reduced
computational complexity were proposed, which employ the cross-terms
based on both feedback and feedforward polynomials and they can
model NLANC systems with shorter filter length
\cite{le2019m-max,guo2018adaptive}. In \cite{le2018bilinear}, a new
bilinear filter was proposed, which uses FLANN to expand the input
vector and then adapts coefficients via bilinear filter to obtain
excellent nonlinear modeling capability.

In this context, the particle swarm optimization (PSO) algorithm is
an effective biologically inspired meta-heuristic algorithm, which
performs by a competition and cooperation mechanism between the
particle swarm and thus it retains the global search strategy of
population \cite{rout2016particle}. Such algorithm has been
extensively studied in ANC technique, especially for nonlinear
system \cite{rout2012particle}. The most obvious advantage of the
PSO-ANC algorithm is that it can reduce noise without online
estimation of the secondary path
\cite{rout2012particle,rout2019pso}. In addition, other heuristic
ANC algorithms were also developed, see \cite{khan2019a} and
references therein.

Since 2009, several attractive methods were proposed for NLANC. For
performance improvement, the kernel adaptive filter (KAF) algorithm
and spline adaptive filter were introduced
\cite{bao2009active,scarpiniti2013nonlinear}. The KAF recasts the
input data into a high-dimensional feature space via a reproducing
kernel Hilbert space (RKHS) and as such, the linear adaptive filter
is applied in the feature space \cite{bao2009active}. In contrast,
the spline adaptive filter establishes a nonlinear mathematics model
by an adaptive look-up table (LUT) in which the control points are
interpolated by a local low-order polynomial spline curve
\cite{scarpiniti2013nonlinear}. An important merit of the spline
adaptive filter is that its computational complexity is much lower
as compared to other traditional NLANC algorithms
\cite{patel2016compensating}. Motivated by distributed algorithms,
some nonlinear distributed algorithms were proposed by integrating
FLANN into distributed ANC algorithms, which is important for
ensuring that the wireless acoustic sensor networks (WASNs) can
combat nonlinear distortions \cite{kukde2017on}.

Initial studies have shown that the KAF is limited by high
computational complexity. Specifically, the input data of KAF must
be a dictionary, where every new data input that arrives is used to
calculate the filter output \cite{Weifeng2010Kernel}. To overcome
this limitation, the quantized scheme \cite{chen2011quantized} and
the set-membership scheme \cite{flores2019set-membership} can be
taken into account to start new research on the field of NLANC. For
example, the adaptive step size scheme and data-selective update
strategy proposed in \cite{flores2019set-membership} can be directly
combined with the KAF in NLANC. On the other front, Internet of
Things (IoT) has been shown to be compatible with communication and
signal processing techniques for creating a smart world. In
\cite{shen2018mute}, the concept of ANC was introduced to IoT,
generating \textit{MUTE leverages}, for environmental noise
cancellation. In \cite{galambos2015active}, the concept system of
ANC-IoT is tested. All these works have demonstrated the feasibility
and great potential of the combination of the IoT and ANC.

The goal of this survey is to provide a comprehensive overview of
the latest NLANC research results and innovations. We also discuss
the applications of ANC techniques that have been introduced in many
fields, such as functional magnetic resonance imaging (fMRI)
acoustic noise control \cite{reddy2011hybrid} and transformer noise
control \cite{zhang2012active}. Moreover, we present some remarks on
future research directions for ANC. The papers that have been
published since 2009 are considered in this survey. In particular,
the main problems that the NLANC algorithms should address are
presented in Section \ref{sec:b2}. The modeling method of NLANC is
discussed in Section \ref{sec:2}, where Volterra ANC algorithms,
FLANN-based algorithms, and bilinear ANC algorithms are used as
examples. Section \ref{sec:3} focuses on
evolutionary-computing-based ANC algorithms. In Section \ref{sec:4},
some novel NLANC methods emerging in the past decade are
investigated. Some important implementation issues and future
research challenges of ANC are discussed in Sections \ref{sec:5} and
\ref{sec:6}, respectively. Finally, conclusions are provided in
Section \ref{sec:7}. All the brackets follow the order $\{ [
(\cdot)]\}$.

\section{Problem formulation}
\label{sec:b2}

In practical situations, three types of nonlinearities may exist in
ANC systems. First of all, the primary noise at the cancellation
point $d(n)$ may exhibit nonlinear effects. For instance, when the
primary noise propagates in a duct with very high sound pressure,
the transfer function of $P(z)$ can be well modeled by nonlinear
function \cite{chen2019nonlinear,luo2016efficient}.

On the other hand, the secondary plant $S(z)$ models the signal
converters (A/D and D/A), power amplifiers, and transducers
(actuators or speakers). The nonlinear effect may be caused by
overdriving the electronics or the speakers/transducers in the
secondary path \cite{napoli2009nonlinear}, and the amplifier may be
in saturation \cite{kuo2010active}. In these cases, the nonlinear
characteristics should be considered in the mathematical model for
enhanced performance. In research experiments, the transfer function
between the speaker and the error microphone is often modeled by
\textit{minimum} or \textit{non-minimum phase} response.

Apart from the above two cases, some nonlinear distortions may come
from the components of the control system. The most common source of
nonlinearity in ANC systems is the actuator. The control actuator
has typically nonlinear response characteristics, when it excites
the frequency of interest and its related harmonics at the same
time, and it is hoped that the ANC system will compensate and
overcome such adverse effect. Furthermore, the noise arises from
dynamic systems and as such, the noise may be a nonlinear and
deterministic signal \cite{lu2016adaptive}. For example, it is known
that the noise from blowers, grinders, airfoils, and fans is shown
to be chaotic \cite{das2012nonlinear,behera2014nonlinear}. The type
and examples of nonlinearities present in ANC systems are summarized
in Table \ref{Table00}.

\begin{table}[]
    \scriptsize
    \centering \doublerulesep=0.05pt
    \caption{The type and examples of nonlinearities present in ANC systems.}
    \begin{tabular}{l|l}
        \hline
        \textbf{Type of nonlinearities}                                                                & \textbf{Examples} \\ \hline
        Nonlinearities in the primary path                                                            &  \begin{tabular}[c]{@{}l@{}}The primary noise propagates in a duct with very high \\sound pressure. \end{tabular}        \\ \hline
        Nonlinearities in the secondary path                                                            &  \begin{tabular}[c]{@{}l@{}} (1)  Overdriving the electronics or the speakers/transducers \\in the secondary path.\\ (2) The amplifier may be in saturation.   \end{tabular}  \\ \hline
        \begin{tabular}[c]{@{}l@{}}Nonlinearities in the components\\ of the ANC systems\end{tabular}
        & \begin{tabular}[c]{@{}l@{}} (1) Actuator excites the frequency of interest and its \\related harmonics at the same time.\\ (2) Actuator operates with the input signal whose \\amplitude is close to or higher than the nominal input\\ signal value.\\ (3) The working frequency of the actuator is beyond \\the normal frequency range (or close to the limit).\\ (4) Noise from the blowers, grinders, airfoils, and fans\\ is shown to be chaotic. \end{tabular}       \\ \hline
    \end{tabular}
\label{Table00}
\end{table}

In these cases, the linear ANC techniques cannot fully utilize the
coherence in the noise and achieve suboptimal performance even
though most adaptive algorithms
\cite{jidf,spa,intadap,mbdf,jio,jiols,jiomimo,sjidf,ccmmwf,tds,mfdf,l1stap,mberdf,jio_lcmv,locsme,smtvb,ccmrls,dce,itic,jiostap,aifir,ccmmimo,vsscmv,bfidd,mbsic,wlmwf,bbprec,okspme,rdrcb,smce,armo,wljio,saap,vfap,saalt,mcg,sintprec,stmfdf,1bitidd,jpais,did,rrmber,memd,jiodf,baplnc,als,vssccm,doaalrd,jidfecho,dcg,rccm,ccmavf,mberrr,damdc,smjio,saabf,arh,lsomp,jrpaalt,smccm,vssccm2,vffccm,sor,aaidd,lrcc,kaesprit,lcdcd,smbeam,ccmjio,wlccm,dlmme,listmtc,smcg,mfsic,cqabd,rmmse,rsthp,dmsmtc,dynovs,dqalms,detmtc,1bitce,mwc,dlmm,rsbd,rdcoprime,rdlms,lbal,wlbd,rrser},
can be adapted to NLANC scenarios. To tackle these limitations, both
the models and the algorithms utilized in ANC must be properly
re-designed, resulting in NLANC. There are two problems which should
be addressed in NLANC: what are the practical performance
characteristics of nonlinear model in this implementation, and what
are the details of an algorithm which can be used to adapt a
nonlinear architecture for NLANC implementation.

\section{NLANC algorithms}
\label{sec:2}

The modeling accuracy of NLANC systems is crucial for meeting the
performance requirements of noise reduction level. This section
presents different approaches for the design of NLANC.

\begin{figure}[!htb]
    \centering
    \includegraphics[scale=0.8]{Fig01}
    \caption{\label{01} Diagram of feedforward NLANC system.}
    \label{Fig01}
\end{figure}

\subsection{Volterra ANC algorithms}

Nonlinear autoregressive with exogenous variables (NARX) models have been preliminarily considered for NLANC techniques. Such methods are flexible and satisfy the representation capability of nonlinear systems. In addition, the truncated Volterra filter and bilinear filter can be subsumed into the NARX model. The potential superiority of the narrowband ANC (NANC) model for NLANC applications has been emphasized in \cite{napoli2009nonlinear} and the improved version that considers the commutation of the control filter block and the secondary path dynamics was proposed in \cite{delvecchio2014dual}.

The Volterra ANC algorithms have been studied by several researchers. The preeminent merit of this modeling approaches is that they are `universal approximators' that models diverse systems with satisfied performance. The block diagram of the NLANC system is shown in Fig. \ref{Fig01}, where $e(n)$ denotes the residual noise, $\hat S(z)$ is the estimate of $S(z)$, $x(n)$ stands for the reference signal, $d(n)$ denotes the undesired signal, and $\bm x_e(n)$ denotes the expanded input vector. The Volterra ANC system utilizes truncated Volterra series as a nonlinear expansion. The filter output $y(n)$ of these algorithms can be represented by
\begin{equation}
y(n) = \sum\limits_{q = 1}^Q {\sum\limits_{m_1 = 0}^{M-1} {\cdots\sum\limits_{m_q = 0}^{M-1} {{h_q}(m_1,...,m_q)\prod\limits_{i = 1}^q {x(n - m_i)}}}}
\label{002}
\end{equation}
where $Q$ denotes the order of nonlinearity, $M$ is the memory length, and $h_q(\cdot)$ is the $q$th-order Volterra kernel. The filter output of (\ref{002}) can be compactly rewritten as
\begin{equation}
y(n) = \bm h^{\mathrm T}\left[\bm x^{\mathrm T}_1(n),\bm x^{\mathrm T}_2(n),\ldots,\bm x^{\mathrm T}_Q(n)\right]
\label{002b}
\end{equation}
where $\bm h$ denotes the filter coefficient vector, $(\cdot)^{\mathrm T}$ is the transposition, $\bm x_1(n)=[x(n),x(n-1),\ldots,x(n-M+1)]^{\mathrm T}$, and the th$j$-order input vector $\bm x_j(n)$ can be gained recursively by
\begin{equation}
\bm x_j(n) = \bm x_1(n) \otimes \bm x_{j-1}(n),\;\;j=2,3,\ldots,Q
\label{002c}
\end{equation}
where $\otimes$ represents the Kronecker product.

In the general case with memory size $M$ and nonlinearity order $Q$, the number of coefficients  can be computed by $N_c = \frac{(M+Q)!}{M!Q!}-1$. We can observe that the number of model coefficients grows in $M^Q$, causing an exponential dimension increase. Therefore, using a truncation strategy to limit the value of $Q$ is indispensable in practice. For $Q=2$, the SOV-ANC algorithm is obtained, whose filter output can be expressed as follows \cite{lu2016adaptive}:
\begin{equation}
\begin{aligned}
y(n) = {\bm h}^\mathrm{T}{\bm x}_e(n)  = \sum\limits_{m_1 = 0}^{M-1} h_{1}(m_1)x(n-m_1) \\
+ \sum\limits_{m_1 = 0}^{M-1}\sum\limits_{m_2 = m_1}^{M-1} h_{2}(m_1,m_2)x(n-m_2)x(n-m_1)
\end{aligned}
\label{003}
\end{equation}
where the expanded input vector and filter coefficient vector can be defined by
\begin{equation}
\begin{aligned}
{\bm x}_e(n) \triangleq&\; \left[\bm x^\mathrm{T}_1(n),\bm x^\mathrm{T}_2(n)\right]^\mathrm{T}\\
=&\; [x(n),...,x(n-M+1),\\
&x^2(n),x(n)x(n-1),...,x^2(n-M+1)]^\mathrm{T},
\label{004}
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
{\bm h} \triangleq&\; [h_{1}(0),h_{1}(1),...,h_{1}(M-1),\\
&h_{2}(0,0),h_{2}(0,1),...,h_{2}(M-1,...,M-1)]^\mathrm{T}.
\label{005}
\end{aligned}
\end{equation}

To obtain a stable performance in a high impulsive process, for example, $\alpha=$1.1 or 1.2, two novel algorithms, termed as VFxlogLMP and VFxlogCLMP, were proposed \cite{lu2016adaptive}. These algorithms are derived from the $l_p$-norm of logarithmic cost and outperform than the conventional Volterra FxLMS (VFxLMS) and Volterra filtered-x least mean $l_p$-norm (VFxLMP) algorithms. Furthermore, it is shown in \cite{he2019an} that, by using the maximum correntropy criterion (MCC), the Volterra filtered-x recursive maximum correntropy (VFxRMC) algorithm can obtain good nonlinear capability with normalized and adaptive kernel size schemes.

\begin{figure}[!htb]
    \centering
    \subfigure[]{
    \includegraphics[scale=0.85] {Fig02a}}
    \label{fig:subfig:aa}
    \subfigure[]{
    \includegraphics[scale=0.85] {Fig02b}}
    \label{fig:subfig:bb}
    \subfigure[]{
    \includegraphics[scale=0.85] {Fig02c}}
    \label{fig:subfig:cc}
    \caption{\label{02} Structure of the nonlinear model, where N denotes the nonlinear part and L denotes the linear system. (a) Hammerstein model; (b) Hammerstein-Wiener model; (c) Wiener-Hammerstein model.}
    \label{Fig02}
\end{figure}

\subsection{Hammerstein ANC algorithms}

As the type of point-wise nonlinear expansions, the Hammerstein system model is regarded as the class of truncated Volterra filters. The basic Hammerstein system (N-L) model is comprised of a cascade connection of a memoryless nonlinearity in series with a linear system (see Fig. \ref{Fig02} (a)). The input-output mapping of the memoryless nonlinearity can be expressed in the form of polynomial expansion
\begin{equation}
\begin{aligned}
g(n)=p_1 x(n) + p_2 x^2(n) + \ldots + p_M x^M(n)
\label{006}
\end{aligned}
\end{equation}
where $p_j$, $j=1,2,\ldots,M$ denotes the $j$th order coefficient. Combining the Wiener model (L-N), Hammerstein-Wiener model (N-L-N) is closer to the actual nonlinear system than Hammerstein model. The structure of Hammerstein-Wiener model is shown in Fig. \ref{Fig02} (b). Although Hammerstein-Wiener type nonlinear system has one more nonlinear link than Hammerstein and Wiener models, its essential feature is still a series of static nonlinear link and dynamic linear link. Changing the cascading order of the static nonlinear part and the dynamic linear part, the Wiener-Hammerstein model (L-N-L) can be obtained, as shown in Fig. \ref{Fig02} (c). The above mentioned models have been applied to NLANC system, see \cite{sra2017nl}.

\subsection{FLANN-based algorithms}

\begin{figure}[!htb]
    \centering
    \includegraphics[scale=0.8] {Fig03}
    \caption{\label{03} Structure of the FLANN.}
    \label{Fig03}
\end{figure}

\textit{1) FsLMS-based algorithms:} Since the classical FsLMS algorithm was proposed in 2004, many FsLMS-type algorithms integrated with the FLANN controller were presented in the past decade, of which the trigonometric expansion is the most commonly used in FLANN. Fig. \ref{Fig03} shows the structure of the FLANN, where $h_i(n)$ is the coefficient of the weight vector at time instant $n$. The reference signal and weight vector can be defined as follows \cite{reddy2009fast}:
\begin{equation}
\begin{aligned}
{\bm x}_e(n) \triangleq&\; \{x(n),\mathrm{sin}[\pi x(n)],\mathrm{cos}[\pi x(n)],...,\\
&\mathrm{sin}[b\pi x(n)],\mathrm{cos}[b\pi x(n)],...,x(n-M+1),\\
&\mathrm{sin}[\pi x(n-M+1)],\mathrm{cos}[\pi x(n-M+1)],...\\
&\mathrm{sin}[b\pi x(n-M+1)],\mathrm{cos}[b\pi x(n-M+1)]\}^\mathrm{T},
\label{007}
\end{aligned}
\end{equation}
\begin{equation}
{\bm h}(n) \triangleq [h_1(n),h_2(n),...,h_{M(2b+1)}(n)]^\mathrm{T}
\label{008}
\end{equation}
where $b$ is the order of the FLANN filter. The number of coefficients of the FLANN filter can be calculated by $N_c = M(2b+1)$. From this result, we can see that the main advantage of FLANN-based algorithms is that it owns less computational burden than the Volterra ANC algorithm.

In \cite{reddy2009fast}, the FsLMS algorithm was devised for multi-channel NLANC. To achieve moderate computational complexity, the data block processing is usually preferred. In \cite{kumar2010block}, a fast Fourier transform (FFT)-based block FsLMS algorithm was developed. A RFsLMS algorithm was proposed in \cite{george2012a}, which minimizes a new logarithmic cost to obtain enhanced performance in both Gaussian and impulsive noise environments. Moreover, it extended the estimation in \cite{bergamasco2012active} to online update parameter. To further enhance performance, two novel NLANC algorithms were proposed based on the $q$-gradient method \cite{yin2018functional}. One of the algorithms was designed for Gaussian noise with time-varying $q$ scheme; the other algorithm was intentionally designed for AINC.

\textit{2) Generalized FLANN:} Most of the works focus on developing new FLANN algorithms, and little attention is paid for improving FLANN itself. In \cite{sicuranza2011a}, an interesting trial was conducted by adding suitable cross-terms between the input signal $x(n)$ and ${\rm {sin}}(b\pi x(n))$ ${\rm {cos}}(b\pi x(n))$ with different time shifts, resulting in generalized FLANN (GFLANN). Following this main idea, other improved GFLANNs were also developed in \cite{luo2018improved,le2019hierarchical}.

\textit{3) Exponential FLANN:} Inspired by Taylor series expansion that contains all odd powers of sine signal and includes all even powers for cosine signal, an exponential FLANN (EFLANN) filter was proposed for NLANC, which adds the exponential factor in trigonometric functional expansion to accelerate convergence speed \cite{patel2016design}. After that, some algorithms were developed for further enhancement of the EFLANN, such as generalized EFLANN \cite{le2018generalized} and approximated EFLANN \cite{deb2019design}.

\begin{figure}[h]
    \centering
    \includegraphics[scale=0.8] {Fig04}
    \caption{\label{04} Structure of the RFLANN, where $\bm \beta$, $\bm \delta$, and $\bm \eta$ are the parameter vector and $N_y$ is the recent samples of the output signal.}
    \label{Fig04}
\end{figure}

\textit{4) Recursive FLANN:} Recalling linear ANC systems, the filtered-u LMS (FuLMS) algorithm is capable of updating coefficients of the infinite impulse response (IIR) filter. In \cite{george2012development}, the IIR filter with FuLMS algorithm was convexly combined with an FsLMS filter, resulting in the FsuLMS algorithm and allowing performance enhancement with respect to its component filters. Based on a bounded-input bounded-output (BIBO) stability condition, the related work \cite{sicuranza2011bibo} has confirmed that the FuLMS algorithm can also apply for NLANC. In this case, the recursive FLANN (RFLANN) filter was proposed that uses past output as part of the input signal. The block diagram of the RFLANN is shown in Fig. \ref{Fig04}, where $\bm h^\mathrm T=[\bm \beta^\mathrm T, \bm \delta^\mathrm T,\ldots,\bm \eta^\mathrm T]$.  As a flexible extension of LIP filters, the RFLANN has become another important method. An application of the even mirror Fourier nonlinear (EMFN) filter to RFLANN has been considered for NLANC \cite{sicuranza2013new}.

\subsection{Orthogonal basis function-based algorithms}

The above Volterra ANC algorithms are derived by Volterra series  that do not satisfy
orthogonal properties. In a comparison, the orthogonal property of these LIP filters can ensure fast convergence of gradient descent algorithms, thereby improving performance over non-orthogonal LIP filters. In what follows, we review the advances of some orthogonal basis function-based algorithms since 2009.

\textit{1) Chebyshev ANC algorithms:} The Chebyshev nonlinear (CN) filter based on Chebyshev polynomials, which hold two-fold characteristics as follows. I) Chebyshev polynomials are a family of orthogonal polynomials. II) Chebyshev polynomials converge faster than expansions in other sets of polynomials. The Chebyshev polynomials are given by the generating function \cite{chen2019nonlinear}
\begin{equation}
T_{q+1}(x) = 2xT_q(x) - T_{q-1}(x),
\label{009}
\end{equation}
where $x \in [-1,1]$, $T_0(x)=1$, and $q$ is the order of Chebyshev polynomials expansion. The closed-form expression for Chebyshev polynomials of any order is expressed as
\begin{equation}
T_q(x) = \sum\limits_{i = 0}^{[\frac{q}{2}]} {(-1)^i} \binom{q}{2i}x^{q-2i}(1-x^2)^i
\label{010}
\end{equation}
where $[q/2]$ denotes the integer part of $(q/2)$. The corresponding function expansion using Chebyshev polynomials (increased from 1 to $Q$) is expressed as (\ref{011}).

The novel CN filter was proposed in \cite{chen2019nonlinear}, whose input of CN filter is derived from  correlated empirical mode decomposition (CEMD) method. As shown in \cite{chen2019nonlinear}, such Chebyshev ANC algorithm outperforms the known FsLMS and VFxLMS algorithm in noise reduction and convergence rate.

\textit{2) Fourier ANC algorithms:} Despite the FLANN and GFLANN filter significantly reduce the computational complexity, they do not satisfy the condition of the Stone-Weierstrass theorem. Therefore, they cannot be universal approximators of Volterra filters, i.e., these types of filters are unable to model arbitrary systems. In \cite{carini2012new}, a Fourier nonlinear (FN) filter was investigated for NLANC, where trigonometric basis functions are designed to avoid repetition or cancellation terms. However, as discussed in \cite{carini2014fourier}, the discontinuities on the borders of the unit hypercube cause relevant oscillations in the Fourier expansion of continuous function.

To tackle this challenge, the EMFN filter was proposed for approximation, by formulating the truncation of a multidimensional even symmetric generalized Fourier series \cite{carini2014fourier}. For the second-order EMFN filter, the input signal $\bm x_e(n)$ is expanded to $N_c = \frac{M(M+3)}{2}$ coefficients as (\ref{012})\footnote{The length of SOV filter has the same calculation as the second-order EMFN filter.}.
\begin{figure*}[!t]
\begin{equation}
\begin{aligned}
\bm x_e(n) =&\;\left\{1,T_1[x(n)],T_1[x(n-1)],\ldots,T_1[x(n-M+1)], \right.T_1[x(n)]T_1[x(n-1)],\\
&\ldots ,T_1[x(n-M+2)]T_1[x(n-M+1)], T_2[x(n)],T_2[x(n-1)],\ldots,T_2[x(n-M+1)],\\
&\ldots \left. ,T_Q[x(n)],T_Q[x(n-1)],\ldots,T_Q[x(n-M+1)]\right\}^{\mathrm T}
\end{aligned}
\label{011}
\end{equation}
\begin{equation}
\begin{aligned}
\bm x_e(n) = &\left\{ \sin\left[\frac{\pi}{2}x(n)\right],\sin\left[\frac{\pi}{2}x(n-1)\right],\ldots, \right.
\sin\left[\frac{\pi}{2}x(n-M+1)\right],\cos\left[\pi x(n)\right],\\
&\ldots,\cos\left[\pi x(n-M+1)\right],\sin\left[\frac{\pi}{2}x(n)\right]\sin\left[\frac{\pi}{2}x(n-1)\right],\\
&\ldots,\sin\left[\frac{\pi}{2}x(n-M+1)\right]\sin\left[\frac{\pi}{2}x(n-M+1)\right],\ldots,\\
&\sin\left[\frac{\pi}{2}x(n-M+3)\right]\sin\left[\frac{\pi}{2}x(n-M+1)\right],\ldots,
\left. \sin\left[\frac{\pi}{2}x(n)\right]\sin\left[\frac{\pi}{2}x(n-M+1)\right]\right\}^{\mathrm T}
\end{aligned}
\label{012}
\end{equation}
\begin{equation}
\begin{aligned}
\bm x_e(n) =&\;\{L_0[x(n)],L_1[x(n)],\ldots,L_Q[x(n)],\ldots,L_1[x(n-M+1)],L_2[x(n-M+1)],\\
&\ldots,L_Q[x(n-M+1)]\}^{\mathrm T}
\end{aligned}
\label{013}
\end{equation}
\end{figure*}
The EMFN filter was successfully employed to solve NLANC problems. The work \cite{guo2019improved} embedded the EMFN filter inside of the RFLANN structure under the BIBO stability condition. Numerical simulation results indicated that above modified EMFN filters are feasible and effective for solving NLANC problems. The work in \cite{guo2018sparse} developed two novel EMFN and CN filters with sparse modeling of nonlinear secondary path, which can significantly reduce computational load without sacrificing the control performance. Moreover, it showed that the EMFN filter outperforms the CN filter with sparse modeling.

\textit{3) Legendre ANC algorithms:} The Legendre nonlinear (LN) filter is also an orthogonal universal approximators, which shares closed properties with FN and EMFN filters. The Legendre polynomial can be obtained by the following formula \cite{george2013active}:
\begin{equation}
L_{q+1}(x) = \frac{1}{q+1}\left[(2q+1)xL_{q}(x) - qL_{q-1}(x)\right]\;\;q>2,
\nonumber
\end{equation}
where $L_q(x)$ denotes the $q$th order Legendre polynomial derived from the solution of the differential equation $\frac{\rm d}{{\rm d}x}\left[(1-x^2)\frac{{\rm d}y}{{\rm d}x}\right]+q(q+1)y=0$ and $L_0(x)=1$. The expanded input vector of the LN filter is shown in (\ref{013}). In 2013, a new algorithm was suggested in \cite{george2013active}, which uses a cascade structure of FLANN and LN filters. Since the LN filter requires $x \in [-1,1]$ for stability, a function ${\mathrm {tanh}(\cdot)}$ is used to link the output of FLANN and the input of the LN filter. Simulation results demonstrated improved performance as compared to the FsLMS algorithm \cite{george2013active}.

\begin{figure}[!htb]
    \centering
    \includegraphics[scale=0.4] {Simu}
    \caption{Averaged noise reductions (ANRs) of the representative NLANC algorithms, where the filter length of all algorithms is set to 77 and simulation settings are the same as Example 1 in \cite{yin2018functional}.}
    \label{Simu}
\end{figure}

We compare representative NLANC algorithms in Fig. \ref{Simu} for a scenario described in \cite{yin2018functional}, where the LMS algorithm is utilized to update weight vector, and $\alpha=2$ is used as the reference signal. The results show that the SOV filter has the best performance, followed by the LN filter and the CN filter, which have comparable noise residual performance.

\subsection{Bilinear ANC algorithms}

\begin{figure}[!htb]
    \centering
    \includegraphics[scale=0.8] {Fig05}
    \caption{\label{05} Structure of the bilinear filter.}
    \label{Fig05}
\end{figure}

The deteriorated performance of the SOV filter is achieved when strong nonlinear distortions occurred in the ANC system \cite{le2019m-max}. The bilinear filter is developed to address this problem with shorter filter length. The input-output relationship of the bilinear filter is given by
\begin{equation}
\begin{aligned}
y(n) =&\; \sum\limits_{i=0}^{N_1} a_i(n)x(n-i) + \sum\limits_{t=1}^{N_1} b_t(n)y(n-t) \\
&\sum\limits_{i=0}^{N_1}\sum\limits_{t=1}^{N_1}c_{i,t}(n)x(n-i)y(n-t)
\end{aligned}
\label{014}
\end{equation}
where $a_i(n)$, $b_t(n)$, and $c_{i,t}(n)$ are the coefficients with $N_1$, $N_1+1$ and $N_1(N_1+1)$ elements, respectively. The block diagram of the convention bilinear filter for ANC problem is shown in Fig. \ref{Fig05}.

As can be seen, the bilinear filter has $x(n)$, $y(n)$, and $x(n-i)y(n-t)$ with different time shifts. An update method, namely diagonal-channel structure, was proposed in \cite{tan2014nonlinear}. This structure updates the weights of these terms independently to acquire improved performance. By using this structure, the FxLMS and filtered-x recursive least square (FxRLS) algorithms were proposed and a counterpart with sequential channel update was developed \cite{tan2014nonlinear}. Using the properties of the diagonal-channel can also reduce the computational complexity of the bilinear filter \cite{guo2020bibo,dong2020diagonal}. In \cite{tan2016implementation}, two diagonal strategies were devised for the bilinear filter in the context of NLANC system, which omits some diagonal elements to significantly reduce the computational complexity by sacrificing performance of the algorithm.

In recent years, the FLANN-inspired bilinear filters \cite{guo2018adaptive,zhu2019reweighted,le2018bilinear,le2019m-max,luo2018novel} received much attention owing to their good modeling capability. For the ANC system with saturation nonlinearity, a general function expansion bilinear filter with a 3-D diagonal-channel structure was proposed, where a trigonometric-based FLANN expansion diagonal-channel structure bilinear filter is presented \cite{guo2018adaptive}. Aimed at improving the performance of a control system with a nonlinear secondary path, a FxLMS algorithm using reweighted bilinear filters was proposed and a generalized filter to model different nonlinear secondary paths was explored in \cite{zhu2019reweighted}. In \cite{le2018bilinear}, the FLANN was introduced to the bilinear filter for NLANC, which adapts the coefficients by diagonal-channel structure and achieves refined performance than the FLANN and GFLANN. Then, the leaky bilinear filter with filtered-error LMS (FeLMS) and PU scheme was subsequently developed to enhance stability and reduce computational complexity \cite{le2019m-max}. The number of coefficients for the NLANC algorithms based on nonlinear expansion is outlined in Table \ref{Tablebb2}, where $L_F$ denotes the FIR memory size and $\mathcal N_c$ denotes the number of coefficients of the order $Q$.

\begin{table}[]
    \scriptsize
    \centering \doublerulesep=0.05pt
    \caption{Number of coefficients for the algorithms.}
    \begin{tabular}{l|l|l}
        \hline
        \textbf{Nonlinear model}    & \textbf{Number of coefficients} & \textbf{Special instruction} \\ \hline
        Volterra filter    &$N_c = \frac{(M+Q)!}{M!Q!}-1$                       & SOV filter: $N_c = \frac{M(M+3)}{2}$           \\ \hline
        Hammerstein filter &$N_c = M+L_F$                          &\begin{tabular}[c]{@{}l@{}}The nonlinear part of Hammerstein filter does\\ not require historical input and output information \end{tabular}                   \\ \hline
        FLANN filter             &$N_c = M(2b+1)$                      &Based on trigonometric expansions            \\ \hline
        GFLANN filter          &$N_c = 2M+1$                          &Generally, $b=1$ in the GFLANN filter                  \\ \hline
        EFLANN filter            &$N_c = M(2b+1)+1$                           &Require exponential operation                     \\ \hline
        CN filter   &$N_c = 1+Q(M+1)+M$                          &Generally, $Q \leq 3$ in the CN filter                      \\ \hline
        FN filter     &\begin{tabular}[c]{@{}l@{}}$\mathcal N_c = 2\sum\limits_{j=0}^{Q-1}\tbinom{Q-1}{j}\tbinom{M+j}{Q}$\\(for $Q$th order)\end{tabular}
                                  &\begin{tabular}[c]{@{}l@{}}Second-order ENFN filter: $N_c = \frac{M(M+3)}{2}$ \\  ENFN filter is suitable for modeling strong nonlinearities  \end{tabular}             \\ \hline
        LN filter    &$N_c = QM+1$                       &\begin{tabular}[c]{@{}l@{}}Suitable for modeling mild and medium\\ nonlinearities\end{tabular}            \\ \hline
        Biliner filter     &$N_c = N_1^2+3N_1+1$          &\begin{tabular}[c]{@{}l@{}}Suitable for modeling strong
        nonlinearities\\ with saturated signals \end{tabular}                     \\ \hline
    \end{tabular}
\label{Tablebb2}
\end{table}

\subsection{Discussion}

In this subsection, we discuss some existing strategies in NLANC systems. Most of these schemes have been reviewed in Part I of this work, and can be extended to solve the ANC problem under nonlinear distortions.

The remote microphone technique (RMT) has been used for NLANC in the past decade. In \cite{spiriti2014gradient}, the RMT is incorporated into the empirical weight update (EWU) algorithm to circumvent the slow convergence problem of stochastic approximation methods and PSO algorithms. Based on FLANN controller, a novel FeLMS algorithm with RMT was proposed in \cite{das2012adjoint}, which exhibits its computational advantages and reliable performance for the sinusoidal signal. Its FsLMS counterpart also has been proposed in \cite{das2012nonlinear}.

As stated in Part I, the convex combination scheme can substantially improve the performance of linear ANC. Such a strategy is also applicable to NLANC systems. The first attempt of convex combination for NLANC can be found in \cite{george2014convex}. In this work, two different nonlinear modeling methods, FLANN and Volterra, were performed by convex combination to obtain enhanced performance.

The MCC not only has been used to the Volterra FxRLS algorithm \cite{he2019an}, but also has been applied to normalized versions of the FsLMS algorithm \cite{kurian2017robust}. As a generalization of the minimum error entropy (MEE) criterion, R{\' e}nyi's entropy has been applied to NLANC. A generalized filtered-x integrated-evaluation-criterion-based gradient descent (GFx-IECGD) algorithm was presented, whose update equation is mixed by the mean square error (MSE) criterion and R{\' e}nyi's entropy \cite{zhang2010adaptive}.

The chaotic noise commonly exists in NLANC system, and some researchers developed various algorithms to control such type of noise. The work in \cite{behera2014functional} considers ANC of a mixture of tonal and chaotic noise via the hybrid FLANN method. The effort in \cite{behera2014nonlinear} first explored the characteristic of logistic chaotic noise and applied FLANN for active chaotic noise control. It was observed that the logistic chaotic noise 1 is highly broadband, while the logistic chaotic noises 2$\sim$6 are mostly broadband with some frequency notches \cite{behera2014nonlinear}. To enhance readability, we summarize these schemes in Table \ref{Table01}.

\begin{table}[htp]
    \scriptsize
    \centering \doublerulesep=0.05pt
    \caption{Contributions of the algorithms in Section \ref{sec:2}.}
    \begin{tabular}{c|c|c}
        \hline\hline
        \textbf{RMT}                    & \begin{tabular}[c]{@{}l@{}}\textbf{Convex Combination}\\ \textbf{algorithms}\end{tabular} & \begin{tabular}[c]{@{}l@{}}\textbf{Entropy-based}\\ \textbf{algorithms}\end{tabular} \\ \hline \cite{spiriti2014gradient,das2012adjoint,das2012nonlinear}
        &\cite{george2012development,george2014convex,guo2019convex}                 &\cite{he2019an,kurian2017robust,zhang2010adaptive}                                                                    \\ \hline
        \begin{tabular}[c]{@{}l@{}}\textbf{Chaotic signal}\\ \textbf{control}\end{tabular}
        & \textbf{Leaky algorithms}
        & \begin{tabular}[c]{@{}l@{}}\textbf{PU-based} \\ \textbf{algorithms}\end{tabular}     \\ \hline
        \begin{tabular}[c]{@{}l@{}}\cite{guo2019convex,le2019m-max,le2018bilinear,luo2018novel}\\ \cite{das2012nonlinear,behera2014functional,behera2014nonlinear,luo2017novel}\end{tabular}
        &\cite{le2019m-max,luo2018novel}                    &\cite{le2019m-max}                                                                    \\ \hline\hline
    \end{tabular}
\label{Table01}
\end{table}

\section{Heuristic-based ANC algorithms}
\label{sec:3}
Both linear and nonlinear ANC tasks can be solved by transforming into a global optimization problem. Especially for NLANC systems, such problems are well-suitable for considering global optimization problems, because the nonlinear secondary path may cause the conventional adaptive filter to suffer from the local minima trapping phenomenon in non-convex or non-deterministic polynomial (NP)-hard problems, and it turns out the degraded reduction performance. The first heuristic algorithm used to optimize ANC was the genetic algorithm (GA) in 1994 \cite{wangler1994genetic}. Afterward, the PSO algorithm was applied to the NLANC problem in 2006. The relevant time line can be referred to Table \ref{Table000}.

\textit{1) Genetic algorithms:} The GAs are inspired by Darwin's evolution theory, which seeks the optimal solution by studying the evolutionary process of simulation. In \cite{chang2010active}, the adaptive GA (AGA)-based ANC system was designed for various environments. Besides, this work confirmed that the AGA can suppress the nonlinear noise interference without estimating the secondary path. The use of the interior-point method (IPM) for GA-ANC has also been proved to be an effective technique \cite{raja2018bio}. The IPM searches for the optimal solution by traversing the internal feasible region, and it shows excellent performance for solving linear or quadratic constraint optimization problems. Benefiting from this, GA-ANC systems can further enhance reduction performance for various noises.

\textit{2) Backtracking search algorithms:} Closer to the GA is the backtracking search algorithms (BSA), which are a population-based iterative evolutionary algorithm (EA) designed to prevent the algorithm from falling into a local optimum. In \cite{khan2018backtracking}, the BSA integrated with sequential quadratic programming (SQP) was proposed for ANC systems with sinusoidal signal and complex random signal, which can obtain improved performance without modeling $S(z)$.

\textit{3) PSO algorithms:} Based on the competition and cooperation mechanism between the particle swarm, the PSO can produce a swarm intelligence to guide the optimization search. By exploiting the soft bound ${\rm {tanh}}\{\cdot\}$ as the saturation nonlinearity module of the secondary path, the PSO has been applied to ANC in \cite{rout2016particle}. Specifically, it can tackle the cases with varying saturation grades of both reference and secondary paths. Similar to the GA, the PSO-ANC can also reduce noise without exact estimation of the secondary path, see \cite{rout2012particle,rout2019pso,raja2019design}. As a nonlinear FLANN extension of \cite{rout2012particle}, the algorithm in \cite{george2012particle} was trained by PSO for multi-channel NLANC systems, which is subsequently developed into a decentralized version. In the PSO-literature, a less common exploration strategy is the Wilcoxon norm \cite{george2012robust2}, which has been proved to have effectiveness for data contaminated by outliers. It has been shown through simulation experiments in \cite{george2012robust2} that considerable reduction can be obtained as compared with the MSE criterion.

\textit{4) Bacterial foraging optimization algorithms:} The bacterial foraging optimization (BFO) is based on the foraging behavior of \textit{escherichia coli} bacteria in human intestine. Bacteria look for nutrients to maximize the energy gained in a unit of time. They go through different stages of chemotaxis, swarming, reproduction and elimination-and-dispersal. Such BFO algorithms have been successfully employed to adapt the coefficients of the NLANC controller, which outperforms the standard GA-ANC by nearly 5dB in steady-state \cite{gholami2012non,gholami2014active1}.

In addition to the above mentioned heuristic algorithms, other heuristic algorithms have been studied in the ANC system. The firefly (FF) algorithm was a classical optimization and has been used to NLANC. The new algorithm cascades the FLANN and FIR filter, and adopts the FF algorithm to obtain coefficients \cite{walia2018design}. The fireworks algorithm (FWA) is also a swarm intelligence algorithm inspired by the fireworks explosion. The research in \cite{khan2019a} reveals the profound optimization capacity of the three variants of FWA for NLANC.

\section{Novel NLANC methods emerging in past decade}
\label{sec:4}
In this section, we present some novel NLANC approaches, such as spline ANC algorithms, KAF ANC algorithms, and nonlinear distributed ANC algorithms, proposed in the last ten years.

\subsection{Spline ANC algorithms}
\begin{figure}[!htb]
    \centering
    \includegraphics[scale=0.8] {Fig06}
    \caption{\label{06} Structure of the NLANC system with spline adaptive filter.}
    \label{Fig06}
\end{figure}

In 2013, the spline adaptive filter was originally proposed for nonlinear system identification, whose adaptation is based on a spline function \cite{scarpiniti2013nonlinear}. As shown in \cite{scarpiniti2013nonlinear} and references therein, the spline adaptive filter enjoys good nonlinear modeling performance in diverse applications and owns low computational cost. Such spline adaptive filter was extended to address NLANC problem in \cite{patel2015nonlinear}. The schematic of the NLANC system with a spline adaptive filter is shown in Fig. \ref{Fig06}, where $\bm U(n)$ is the filtered signal of $\bm u(n)=[u^3(n),u^2(n),u(n),1]^{\mathrm T}$, and $\varepsilon$ denotes the span index, and $u(n)$ is the local parameter, which can be computed by
\begin{subequations}
    \begin{equation}
    \varepsilon = \left\lfloor \frac{y_s(n)}{\Delta x} \right\rfloor + \frac{C-1}{2}
    \end{equation}
    \begin{equation}
    u(n) = \frac{y_s(n)}{\Delta x} - \left\lfloor \frac{y_s(n)}{\Delta x} \right\rfloor
    \label{015ab}
    \end{equation}
\end{subequations}
where $y_s(n)$ denotes the output of the finite impulse response (FIR) filter, $\Delta x$ denotes the gap between the control points, $C$ denotes the total number of control points in the activation function, and $\lfloor \cdot \rfloor$ stands for the floor operation. Then, the generated parameters pass through the adaptive look-up table (LUT) and are interpolated by a local low order polynomial spline curve. Hence, the output of the spline adaptive filter is
\begin{equation}
y(n) = \bm u^{\mathrm T}(n) \bm \Omega \bm q
\label{016}
\end{equation}
where $\bm \Omega$ is the pre-computed spline basis matrix and $\bm q$ is the control point vector.

As an initial step towards multi-channel spline adaptive filters, the algorithm proposed in \cite{patel2020multi}  has obvious advantages in terms of computational complexity and final MSE level when compared to the multi-channel VFxLMS and FsLMS algorithms. Furthermore, the advantages of the FIR spline adaptive filter can also refer to \cite{patel2016design2}. In 2016, a novel infinite impulse response (IIR) spline adaptive filter was proposed \cite{patel2016compensating}, which corresponds to the use of a feedback scheme and the FuLMS algorithm. Compared with the FIR filters, the main merit of using IIR filters is their implementation efficiency.

\subsection{KAF-ANC algorithms}

\begin{figure}[!htb]
    \centering
    \includegraphics[scale=0.8] {Fig07}
    \caption{\label{07} Structure of the NLANC system with KAF.}
    \label{Fig07}
\end{figure}
The kernel trick has become rather popular and has been successfully applied to adaptive filtering. The key idea of the KAF is to recast the input data from input space $\mathbb U$ into a high-dimensional feature space $\mathbb F$ via a reproducing kernel Hilbert space (RKHS). As a result, the inner product operations in linear adaptive filters are translated to the calculation of a kernel function $\kappa\left(\cdot\right)$ in the feature space without knowing the exact nonlinear mapping.

The KAF was first applied to NLANC system in 2009 \cite{bao2009active}. The basic block diagram of the ANC with KAF is shown in Fig. \ref{Fig07}, where $\bm \phi: \mathbb U \to \mathbb F$ denotes feature vector. By using \textit{Mercer theorem}, the inner products in RKHS can be calculated as
\begin{equation}
\kappa\left(\bm {X},\bm {X'}\right) = \bm \phi\left(\bm {X}\right)\bm \phi^{\mathrm T}\left(\bm {X'}\right)
\label{017}
\end{equation}
where ${\bm {X}},{\bm {X'}} \in \mathbb U$ and the kernel function $\kappa\left(\cdot\right)$ often uses Gaussian kernel, which is expressed as $\kappa\left(\bm {X},\bm {X'}\right) = {\exp}\left(-\varsigma\left\|{\bm {X}}-{\bm {X'}}\right\|^2\right)$, where $\varsigma$ is the kernel size and $\left\|\cdot\right\|$ denotes the $l_2$-norm. Finally, the filter output can be expressed as
\begin{equation}
y(n) = \sum\limits_{j = 1}^n {{a_j} \langle \bm \phi (\bm X),} \bm \phi (\bm X_j)\rangle
\label{019}
\end{equation}
where $a_j$ denotes the coefficient and $\langle\cdot\rangle$ denotes the inner product operation. The work in \cite{bao2009active} considered the scenario that chaotic noise and only the primary path $P(z)$ is nonlinear, which has certain limitations. The reference \cite{liu2018kernel} made up for this shortcoming, which took the nonlinear secondary path into account. Moreover, some multi-tonal and sinusoidal noise sources are simulated. The choice of the kernel function is crucial to the performance of the kernel algorithm, and most works employed Gaussian kernel as a default choice due to its universal approximating capability, numerical stability, and desirable smoothness. The work in \cite{raja2019novel} specifically discussed the selection of different kernel functions. Three novel kernel functions, namely logistic kernel, tan-sigmoid kernel, and inverse-tan kernel, were proposed and subsequently applied to KAF for NLANC. Its simulation results demonstrated the improved performance by using new kernel functions.

It should be noted that the above mentioned KAF and spline adaptive filters have been exhibited improved performance as compared to the conventional FsLMS and VFxLMS algorithms \cite{patel2015nonlinear,liu2018kernel}. Hence, we can conclude that the KAF and spline adaptive filters are good candidates for NLANC.

\subsection{Nonlinear distributed ANC algorithms}
\begin{figure}[!htb]
    \centering
    \includegraphics[scale=0.7] {Fig08}
    \caption{\label{08} Structure of the nonlinear distributed multi-channel ANC system.}
    \label{Fig08}
\end{figure}

The idea of a nonlinear distributed ANC system is basically derived from a linear distributed ANC network. The difference lies in the use of function expansion strategies. In \cite{kukde2017on}, a nonlinear distributed multi-channel ANC algorithm was proposed followed by the diffusion strategy, where the LN filter is employed and both primary and secondary paths are considered to be nonlinear. Fig. \ref{Fig08} shows the structure of the nonlinear distributed multi-channel ANC algorithm, where $i=1,2,\dots,I$, $I$ denotes the number of loudspeaker, $k=1,2,\dots,K$, and $K$ denotes the number of error sensors. Very recently, a new nonlinear distributed ANC system was presented by using incremental implementation and adaptive EFLANN expansion, which exhibits lower computational cost than the distributed ANC system with the FsLMS, VFxLMS, and conventional EFLANN controller \cite{kukde2020incremental}.

\section{Recent implementations and applications of ANC}
\label{sec:5}


    \begin{table}[]
        \scriptsize
        \centering \doublerulesep=0.05pt
        \caption{Time line of ANC applications.}
            \begin{tabular}{l|l|l|l}
                \hline\hline
                \textbf{Years} & \textbf{Authors} & \textbf{Contributions} & \textbf{References} \\ \hline
                1993      &Craig and Angevine          &\begin{tabular}[c]{@{}l@{}}Active transformer noise \\control\end{tabular}                &\begin{tabular}[c]{@{}l@{}} \cite{craig1993active}\end{tabular}            \\ \hline
                1994      &Sutton et al.          &\begin{tabular}[c]{@{}l@{}}ANC for road noise\end{tabular}                &\begin{tabular}[c]{@{}l@{}} \cite{sutton1994active}\end{tabular}            \\ \hline
                1997      &McJury et al.          &\begin{tabular}[c]{@{}l@{}}ANC for fMRI\end{tabular}                &\begin{tabular}[c]{@{}l@{}} \cite{mcjury1997use}\end{tabular}            \\ \hline
                2002      &Gan and Kuo         &\begin{tabular}[c]{@{}l@{}}Integrated feedback ANC\\for headset\end{tabular}                &\begin{tabular}[c]{@{}l@{}} \cite{gan2002integrated}\end{tabular}            \\ \hline
                2002      &Qiu et al.          &\begin{tabular}[c]{@{}l@{}}Waveform synthesis method\\for active transformer noise\\ control\end{tabular}                &\begin{tabular}[c]{@{}l@{}} \cite{qiu2002waveform}\end{tabular}            \\ \hline
                2003      &Jakob and M{\"o}ser          &\begin{tabular}[c]{@{}l@{}}ANC with double-glazed \\windows\end{tabular}                &\begin{tabular}[c]{@{}l@{}} \cite{jakob2003active,jakob2003active2}\end{tabular}            \\ \hline
                2003      &\begin{tabular}[c]{@{}l@{}}Pinte, Desmet,\\ and Sas\end{tabular}           &\begin{tabular}[c]{@{}l@{}}Iterative learning control \\algorithm for ANC systems\end{tabular}                &\begin{tabular}[c]{@{}l@{}} \cite{pinte2003active}\end{tabular}            \\ \hline
                2005      &Zhou et al.          &\begin{tabular}[c]{@{}l@{}}Back-propagation (BP)\\-FeLMS algorithm\end{tabular}                &\begin{tabular}[c]{@{}l@{}} \cite{zhou2005analysis}\end{tabular}            \\ \hline
                2006      &\begin{tabular}[c]{@{}l@{}}Kuo, Mitra,\\ and Gan\end{tabular}          &\begin{tabular}[c]{@{}l@{}}Feedback ANC for headphone\end{tabular}                &\begin{tabular}[c]{@{}l@{}} \cite{kuo2006active}\end{tabular}            \\ \hline
        \hline\hline
        \end{tabular}
        \label{Table001}
    \end{table}

In this section, the application development of ANC in the past decade is discussed. Before we start, it should be highlighted that many applications have been developed before last decade (see time line of ANC applications in Table \ref{Table001}). These works have laid the foundation for the implementations and applications of ANC in the past ten years.

\subsection{ANC implementation}

\textit{1) DSP implementation:} ANC techniques have been implemented in real-time systems such as digital signal processor (DSP), field programmable gate array (FPGA), and very large scale integration (VLSI). Several hardware implementations dedicated to studying lower-power consumption and small area realization have been reported. The compensated-FxLMS algorithm implemented by fixed-point DSP was presented in \cite{chang2009efficient}. However, such method requires a specialized measurement of the devices in the control system. In addition, the round-off and quantization errors will change with the programmer's code implementation, and the correct values of the phase and amplitude compensation factors need to be determined \cite{kim2017enhancement}. An improvement of fixed-point DSP method was proposed in \cite{kim2017enhancement}, which added feedback and \textit{communication error}\footnote{The definition of communication error can be found in Eq. (12), Part I of this work.} Considering the implementation of different structures in fixed-point DSP, a detailed performance comparison of FxLMS, FxNLMS, IIR, and subband ANC algorithms was presented in \cite{siravara2003comparative}. The research results showed that the FxLMS and its variant are more suitable for ANC implementation and the fixed-point DSP has about 2dB degradation as compared with its floating-point counterpart.

\textit{2) FPGA implementation:} The DSP owns limited processing capability and it is intended for low sampling rate applications. In contrast, the performance of the FPGA implementation largely depends on the selection of FPGA device. Several algorithms and architectures were proposed for performance enhancement with high-speed computational architecture \cite{ardekani2011filtered,shi2017multiple,leva2011fpga}. The novel algorithm, called filtered-weight FxLMS (FwFxLMS), was proposed \cite{ardekani2011filtered}. The FwFxLMS algorithm can efficiently increase the upper-bound of the step size and as such, accelerating the convergence rate of algorithm without affecting the stability. Devoted to reducing computational load of multi-channel ANC system, a novel architecture of FPGA was proposed to solve the conflicted requirement between throughput and hardware resource consumption by using parallel execution and folding \cite{shi2017multiple}.

\textit{3) VLSI implementation:} The VLSI provides a large amount of calculation in a small area, and the resource assumptions is significantly reduced compared to FPGA and DSP. The effort in \cite{mohanty2017hardware} focused on VLSI implementation on the FxLMS and FsLMS algorithms. Prior to this, few literature considered the implementation of NLANC algorithms. By devising folded structures of the algorithm, the intentionally proposed delay FsLMS algorithm becomes easy to use in practice.

\subsection{ANC through open window}

ANC through open windows technique can attenuate noise while maintaining natural ventilation in dwellings \cite{kwon2013interior,huang2011active,wang2019boundary,wang2019boundary2,murao2019hybrid}. The finite-element-based ANC methods have been used to establish the control performance and the passive attenuation offered by the open window \cite{lam2018active,lam2018physical}. In \cite{lam2018active}, an ANC system with 16-channel and secondary sources distributed evenly across the opening was implemented with a full-sized window, yielding more than 5dB below 2000Hz overall attenuation. However, the size of an opening influences the control performance and the opening windows may also introduce noise at the same time \cite{elliott2018wavenumber}. To avoid this issue, a scheme of small opening in the wall was implemented in \cite{zhang2020secondary}. To design an ANC system with reduced computational complexity of DSP implementation and acceptable level of noise reduction, a novel multi-channel ANC method was proposed for opening window. By preprocessing the error signals, the number of inputs to the controller can be decreased and a good performance can be achieved when it employs 4-channel ANC system with 0.2$\times$0.2m opening window \cite{murao2017mixed}.

\subsection{Headphone application}

Many headphone applications of ANC were developed based on DSP, FPGA and VLSI implementations \cite{vu2015low}. Some algorithms focus on implementing ANC headphones with less resources, head-tracking problems, multi-channel, and multi-reference \cite{jung2017combining,behera2017head,cheer2019application}. The most popular algorithm used in ANC headphone is the feedback FxLMS algorithm. A number of (in-ear) headphones based on the feedback ANC system were implemented \cite{chang2010active2,zhang2013intuitive}. Their performance comparison of feedback ANC system with different circuit implementations can be seen in \cite{vu2017high}. In \cite{schumacher2011active}, a hybrid feedback ANC scheme was presented by making use of the combination of classical and adaptive feedback ANC techniques. Considering that system delay will limit implementation performance in practice, this work also proposed a mixed analog-digital realization to maintain very low system delay in the classic part and enable adaptive ANC components to adapt to the digital domain.

In comparison, a feedforward structure was adopted for binaural ANC system in \cite{tanaka2017binaural}. It is worth pointing out that the performance improvement of \cite{tanaka2017binaural} also benefits from the use of the parametric array loudspeakers (PAL) as the secondary sources. Similar techniques have been applied to ANC problems such as steerable PAL \cite{tanaka2010active} and curved-type PAL \cite{tanaka2011mathematically}. Also considering the use of a binaural ANC system, a combined bilateral-binaural ANC (CBBANC) system was proposed for closed-back headphones \cite{belyi2020combined}. By switching between the binaural and bilateral implementation, the improved performance can be obtained as compared to the conventional bilateral ANC headphones.

The ANC headphone application has also embedded some techniques discussed above. For example, in \cite{das2013nonlinear}, the virtual sensing method was adopted. Particularly, in \cite{das2013nonlinear}, the FsLMS algorithm with nonlinear virtual microphone was developed for active headrest application, which considers the nonlinear distortion in the presence of primary noise. The prediction filter was another method applied in feedback ANC headphones, which does not need any real-time coefficient updates and is therefore very economical \cite{guldenschuh2013prediction}.

\subsection{Zone-of-quiet based approach}

The traditional ANC technique can only minimize the noise level at the location where the error microphone is placed. Correspondingly, a zone-of-quiet (ZoQ) is created around this error microphone. However, unless multiple sensors are used, it is not technically possible to control or monitor noise levels in ZoQ. The applications of spherical loudspeakers arrays for ANC has been presented to extend ZoQ \cite{rafaely2009spherical}. Following the idea of \cite{rafaely2009spherical,peleg2011investigation}, several improvements were proposed including linear array \cite{hart2012active}, and RMT-based methods \cite{elliott2015modeling}. To design larger ZoQ, a method was proposed by introducing one and two secondary sources, which aimed at maximizing the area of the 10dB quiet zones \cite{tseng2011local}. By introducing feedforward adaptive ANC structure to control shape and extension of ZoQ, a new avenue was proposed in \cite{ardekani2014adaptive}, which has capability of creating a controllable spherical quiet zone via a compensation factor. Moreover, in \cite{wrona2018shaping}, a ZoQ method equipped with memetic algorithm (MA) and distributed adaptation has been proposed, which holds high potential for practical applications.

The ZoQ based methods have been integrated into ANC for speech enhancement in hearing aids. In \cite{serizel2009zone}, the ANC with noise reduction strategy was proposed for a specific ZoQ. Following a different direction, the method in \cite{serizel2012zone} developed two schemes. The former method is based on the MSE criterion at a remote point, and the latter method is based on an average MSE criterion over a desired ZoQ. Moreover, several approaches on integrated ANC and noise reduction for hearing aids have been proposed, such as those proposed in \cite{serizel2009integrated,pradhan2017speech,albu2018efficient} and references therein.

\subsection{Spatial ANC approach}

Spatial ANC is designed to attenuate noise within a specific space. Since there is not only one space point but the space must be controlled, many transducers are needed. A direct method to achieve spatial ANC is to employ the multi-channel FxLMS algorithm. However, such method cannot guarantee that the spatial noise is also reduced, and the calculation cost is high. To address this drawback, several mode-domain methods were proposed in \cite{zhangg2016sparse,zhang2018active,zhang2018spatial,maeno2018mode}. Such methods can significantly reduce computational complexity since the mode-domain conversions can reduce the cross-correlations between transducers. In addition, the filter update can be computed by a mode-independent manner. Very recently, a new spatial ANC method was developed, where the higher-order source (HOS) scheme is used to improve controllability of efficiency and error \cite{murata2019global}.

\subsection{fMRI acoustic noise control}

The problem of fMRI acoustic noise control is a research focus in the application of ANC technology. Studies showed that when magnetic resonance imaging (MRI) scans, it can generate acoustic noise up to 130dB \cite{reddy2011hybrid}. So far, a large number of ANC algorithms were proposed for fMRI problem, see \cite{kumamoto2011active,kannan2010efficient} and references therein. The followed articles implemented by head-mounted ANC system with piezoelectric loudspeakers and optical microphones \cite{kida2009head} or RMT \cite{miyazaki2015head}. The effectiveness of the FxRLS-FxNLMS algorithms on fMRI noise control was validated that the hybrid algorithm can gain a very fast convergence \cite{reddy2011hybrid}. This group then proposed an improved parallel feedback ANC system \cite{ganguly2015improved}, which parallels two FxNLMS algorithms with the NLMS decomposition to achieve exceeding 40dB noise attenuation level.

The adaptive comb filters (ACFs) that perform a coordinated search, have been applied to fMRI noise control. Based on the previous works, i.e., SONIC (self-optimizing narrowband interference canceller) \cite{niedzwiecki2009new,niedzwiecki2009new2,niedzwiecki2009self}, the ACF incorporated SONIC algorithm was suggested for active fMRI noise control \cite{niedzwiecki2013estimation}. This algorithm can track the fundamental frequency and the amplitude of different frequency components of non-stationary harmonic signals embedded in white measurement noise. However, its reduction level is limited to 20dB. For performance improvement, a discrete Fourier transformation (DFT)-based frequency-domain ANC algorithm for fMRI noise control was proposed in \cite{lee2017frequency}, which can achieve approximately 35dB overall reduction in the frequency range of interest.

\subsection{Periodic noise control}

Based on different characteristic distributions with time, the exhaust noise can be divided into three types: steady, periodic, and intermittent \cite{zhou2013optimal}. Among the three types, the periodic noise is the most widely encountered and has great research significance. The repetitive impulsive noise belongs to a category of periodic noise and has many practical backgrounds. For example, when compressed air escapes regularly through the exhaust pipe, it generates repetitive impulsive noise, up to 120dB. In the past decade, many studies have focused on and analyzed this kind of noise. Similarly, the FxLMS algorithm and its variants can also be used to attenuate this noise \cite{sun2015convergence}.

The iterative learning control (ILC) algorithm is another method that can effectively control repetitive impulsive noise, which can be regarded as the adaptive feedforward algorithm. In ILC controllers, the control signal is adapted to the specific repetitive of transient noise. It runs repetitively in an open-loop over a finite time interval to reduce the same discrete interference at each time instant. Related to the ILC method is the well-known repetitive control (RC) algorithm, which can also be interpreted as the adaptive feedforward algorithm. The ILC and RC algorithms share the same core idea, that is, it should be possible to improve the performance of a system that repeatedly performs the same repetitive task. The difference is that the RC algorithm has continuously active property.

It demonstrated in \cite{pinte2010novel,stallaert2010novel} that the drawbacks of conventional ILC and RC methods are high computational complexity and low robustness of the controllers. To solve these problems, a novel strategy to design non-causal ILC and RC filters was described, with emphasis on the application of ANC. It can achieve 4dB noise reduction with the frequency-domain ILC filter. Owing to the fact that the ILC and RC algorithms are closely related, many design methods can be used in both algorithms. For example, in \cite{zhou2013active}, a non-causal transversal FIR filter was combined with the optimal ILC algorithm for active repetitive impulsive noise control, and then it was applied to the optimal RC algorithm \cite{zhou2013optimal}.

\subsection{Transformer noise control}
Transformer noise is typically characterized by even and multiples of the line frequency of the fundamental and harmonic frequencies. Generally, the fundamental frequency does not vary by more than $\pm$0.5\% \cite{zhang2012active}. Therefore, using the internally synthesized preset frequency signal as a reference signal for transformer noise control may be a reasonable way. In \cite{zhang2012active}, the classical FxLMS algorithm with an internally synthesized reference signal was proposed, which can obtain a maximum of 15dB noise reduction. Motivated by the online and offline secondary path modeling strategy, a novel FxLMS algorithm was presented for 50000KVA transformer noise control \cite{zhao2017new}. Simulations demonstrated that an 84.10$\sim$96.86\% decrease in average sound energy density can be obtained in a certain area.

\subsection{ANC in vehicles}
ANC in vehicles requires a good spatial and temporal matching between the fields of the actuator or secondary sources \cite{samarasinghe2016recent}. Spatial matching requires that the wavelength of the acoustic wave for active noise control must be small, whereas time matching demands the system can adapt to changes in vehicle speed and load \cite{elliott2010active}. Owing to the fact that ANC has a significant improvement, particularly at low-frequency and on vehicles with specific problems, such technique has been widely used in vehicles \cite{cheer2012active}. In the past ten years, many prominent companies and corporations have proposed effective solutions to this problem, see references in \cite{samarasinghe2016recent}. In \cite{cheer2014design}, three feedback ANC systems were investigated based on internal model control (IMC) architecture and frequency discretization method. Moreover, the multi-channel  algorithm has addressed ANC of road noise in vehicles \cite{cheer2015multichannel,jung2019local}. In \cite{jung2019local}, eight reference sensors located around the four wheels were devised to obtain a broadband reduction of 4dB up to 1kHz. In \cite{guo2018active}, a variable step-size (VSS) median-least mean square algorithm was presented for active rail vehicle interior noise control. In \cite{cheer2016active} and \cite{peretti2013adaptive}, the leaky FxLMS and multi-channel FxLMS algorithms have been respectively applied to ANC for yacht.

The ANC technique is quite general and its application is not limited to the above mentioned works. The ANC principle has already been used for mobile phones \cite{cheer2018application}, earmuff systems \cite{moon2015performance,zhu2019design}, high-speed elevators \cite{yang201423}, mufflers \cite{chen2009simulations}, and washing machines \cite{mazur2019active}. In \cite{rohlfing2014ventilation,gardonio2014modular}, the single and multi-channel ANC systems have been employed for ventilation duct. Furthermore, new problems in the ANC system, such as broadband noise fields in rooms \cite{ma2019active}, electroencephalogram (EEG) signal analysis \cite{bagha2018understanding} have also been studied. These applications show the great application potential of ANC.

We note that a number of papers on NLANC assume the presence of nonlinearities that are difficult to encounter in real-world applications. Therefore, most of above mentioned works are based on linear assumptions and do not consider nonlinear characteristics. To clearly show the benefits offered by a nonlinear approach with respect to the linear one, the implementations and applications of NLANC are summarized in Table \ref{Table003}.
\begin{table}[]
    \scriptsize
    \centering \doublerulesep=0.05pt
    \caption{Implementations and applications of NLANC in Section \ref{sec:5}.}
    \begin{tabular}{l|l|l}
        \hline
        \textbf{Contributions}                & \textbf{Benefits}                                          & \textbf{References} \\ \hline
        FsLMS implementation           &\begin{tabular}[c]{@{}l@{}}The proposed delayed block FsLMS (DBFsLMS)\\ structure involves 58\% less area-delay product\\ and 41\% less energy per sample than the FxLMS\\ structure. \end{tabular}         &\cite{mohanty2017hardware}                     \\ \hline
        \begin{tabular}[c]{@{}l@{}}Generalized filter bank based\\nonlinear ANC algorithm\end{tabular}           &\begin{tabular}[c]{@{}l@{}}The new algorithm is suitable for a headrest\\ application with head movement, and the effective\\ area of active headrest is achieved as high as\\ 45 cm $\times$ 45 cm. \end{tabular}         &\cite{behera2017adaptive}                     \\ \hline
        \begin{tabular}[c]{@{}l@{}}NLANC with virtual\\ microphone control\end{tabular}
        &\begin{tabular}[c]{@{}l@{}}These methods can control noise at a remote\\ location when the system is nonlinear.\end{tabular}
        &\cite{das2012nonlinear,das2012adjoint,das2013nonlinear}                     \\ \hline
    \end{tabular}
    \label{Table003}
\end{table}

\section{Future research challenges}
\label{sec:6}

Some future challenges of ANC have been discussed in \cite{kuo2010active,krish}. In this section, we present a number of promising directions for future ANC research from other perspectives.

\subsection{Theoretical analysis of algorithms}

As we discussed in Part I of this work, the fractional lower order moment (FLOM)-based algorithms can be effectively used for AINC. The problem of interest is to analyze the theoretical performance of the FLOM-based algorithms in the presence of $\alpha$-stable noise. However, such a topic is extremely difficult since the variance of $\alpha$-stable noise is infinite, leading to some assumptions does not satisfy the conditions. As far as we know, researchers did not consider the case in which the error in $\alpha$-stable noise environment is changing across iterations. Moreover, the steady-state and transient performance of conventional least mean $p$th power (LMP) algorithm has not been investigated in the presence of $\alpha$-stable noise. The lack of these analysis foundations causes the analysis of the FLOM-based algorithms exactly difficult.

The similar problem also exists in the analysis of LIP algorithms, such as the FN and EMFN. Moreover, since the input vectors of these algorithms are expanded by different basis functions, the behavior of algorithms is more complicated.

\subsection{Sparsification for KAF-ANC}

Despite the above KAF algorithms are well applied in NLANC, there are still problems to be solved. For system identification problem, several KAF algorithms with curb of network size were proposed. Since the network size of the KAF linearly increases with the number of training data, which means that the filter length is very high during adaptation. The problem of NLANC is the same as that of KAF for system identification. This increase in complexity can be harnessed through the application of some sparsification techniques, for example, the quantized scheme in \cite{chen2011quantized} and the set-membership scheme in \cite{flores2019set-membership}. These strategies can curb the number of input data into the kernel function at each iteration. If it is detected that not enough new information is carried in the input signal, it will be abandoned as the redundant data.

\subsection{The application of Internet of Things to ANC}

IoT attracts significant attention in the fields of location, smart transportation, and signal processing. Based on the concept of IoT, the ANC technique has been applied to IoT. By utilizing wireless communication with acoustics, IoT was brought into the ANC context, which provided a potential solution to urban noise pollution \cite{shen2018mute}. Through DSP implementation verification, the proof of concept (PoC) system was stated, which allows the concept of IoT to be applied in the field of ANC \cite{galambos2015active}. In this work, the analog cables were replaced by digital ethernet transmission, resulting in a more flexible system that is easier to be rearranged and modified.

It is noted that existing ANC-IoT schemes rely on the single-node device, and the impact of practical constraints, such as impulsive noise, the nonlinear distortions, and distributed implementation, on the performance of ANC-IoT has not been investigated yet. Furthermore, as discussed before, ANC can effectively suppress the noise in vehicles. Hence, ANC may be specifically applied to the Internet of Vehicles (IoV).

\subsection{Discussion}
Currently, NLANC has made much progress. In the following, some future directions that are peculiar to NLANC are presented.

1) For NLANC algorithms, most of the algorithms are still at the theoretical and laboratory stage, and they are not fully utilized in practical applications. For the future directions of nonlinear algorithms in NLANC, it needs to curb its computational complexity. Since the nonlinear filtering problem is inherently infinite dimensional, the substantial increase of data samples can prohibit practical applications of NLANC algorithms. Therefore, low complexity NLANC algorithms are an important research direction in the future.

2) Selecting the optimal memory length of the NLANC algorithm for the corresponding nonlinear systems is also a research direction. In the field of echo cancellation, the memory control scheme of Volterra filter has been studied  \cite{zeller2010adaptive}. However, it has not been investigated in the NLANC adaptive filtering algorithms.  Such a problem is also applicable to the heuristic-based algorithms in NLANC. The optimal size of populations for heuristic-based algorithms is of great significance to the reduced computational complexity and improved performance.

3) Further exploration of NLANC technology in specific practical scenarios are required. Currently, NLANC does not seem to find its most suitable breakthrough point for large-scale commercial applications. Thus, finding a cost-effective application scenario is also an important future research direction of NLANC.

\section{Conclusion}
\label{sec:7}
In this survey paper, we have discussed state-of-the-art methods related to NLANC technique, including FsLMS-based algorithms, distributed extension and selected applications. We have also summarized FLANN approaches as well as the Chebyshev, EMFN, and LN filters. The spline adaptive filter, kernel adaptive filter, and nonlinear distributed ANC algorithms are the new nonlinear modeling approaches considered for NLANC. Due to technical limitations, these open areas need further investigations of sophisticated NLANC techniques. The review of implementations for both linear and nonlinear ANC are merged together. We have also discussed the main challenges for future ANC techniques.

\section*{Acknowledgment}
\label{sec:Ack}
The authors would like to thank the associate editor and the anonymous referees for their
valuable comments.

\section*{References}









\section{Introduction}

In 2014, the United Nations Joint Programme on HIV and AIDS (UNAIDS) set the 90-90-90 targets as the primary strategy to end the HIV/AIDS epidemic \citep{joint201490}, which includes identifying 90\% of people living with HIV through expanded testing, placing 90\% of positively identified individuals on antiretroviral therapy, and ensuring that 90\% of those on therapy can achieve undetectable viral loads by 2020.
Considerable progress has been made towards Tri90 \citep{Gaolathe2016, gisslen2017sweden, joint201690, joint2017ending, marukutira2018tale}.
Yet, as of 2019, there was a significant gap; instead of the targets of 90-90-90 (Tri90), it was 81-67-59 globally \citep{joint2020good}; and now the COVID-19 crisis has the potential to undermine existing efforts towards the HIV/AIDS epidemic \citep{joint2020seizing}.



To end HIV/AIDS epidemic, we shall learn from the past efforts towards Tri90 and identify important features that could guide more targeted and effective health policies. Novel datasets and sophisticated modeling tools are needed to enhance our understanding of Tri90 achievements. The Population-based HIV Impact Assessment (PHIA) survey is a nationally representative HIV-focused survey that started data collection in 2015. It is designed to measure the reach and impact of HIV programs. We believe it offers a unique opportunity to better understand the key factors that drive the HIV epidemics in the most affected countries in sub-Saharan Africa. In this article, we propose a novel causal structural learning algorithm to discover important covariates and potential causal pathways for the Tri90. Causal structural learning aims to build a directed acyclic graph (DAG) that shows direct causal relations among variables of interest in a given domain. The resulting DAG helps us to understand the mechanism behind data.


Many classical structural learning algorithms are constrained-based, such as the PC algorithm \citep{spirtes1991algorithm}, the PC-stable algorithm \citep{colombo2014order}, and the MMPC algorithm \citep{tsamardinos2003mmpc}. 
The constrained-based algorithms learn graphical structures by d-separation set searching.
\textcolor{black}{
D-separation is an important graphical concept for causal structural learning \citep{geiger1990logic, geiger1990d}, where the ``d" stands for dependence or directed. 
Roughly speaking, in a directed acyclic graph (DAG), two vertices $X$ and $Y$ are d-separated by a set of vertices  $\tf{Z}=\{Z_1,\cdots,Z_d\}$ if any only if all the paths/information between $X$ and $Y$ are blocked by vertices in $\tf{Z}$ (see Section~\ref{section:Preliminaries} for a rigorous definition of d-separation).
The d-separation relationship in a DAG can be related to the conditional independence through the Markov condition and the faithfulness assumption \citep{spirtes2010introduction} (see definitions in Section~\ref{section:Preliminaries}).
That is to say, vertices $X$ and $Y$ are d-separated by a vertex set $\tf{Z}=\{Z_1,\cdots,Z_d\}$ is equivalent to
the conditional independence of the corresponding variables $X$ and $Y$ given the corresponding set of variables $\tf{Z}=\{Z_1,\cdots,Z_d\}$ under the Markov condition and the faithfulness assumption.}

However, existing constrained-based structural learning algorithms can be quite aggressive in edge removal: 
{if a type II error is made such that two connected covariates $X$ and $Y$ are thought to be conditionally independent given some $\tf{Z}$,
the edge between $X$ and $Y$ is removed mistakenly,
and useful information about important features and possible causal pathways get lost during this edge removal process.}
\textcolor{black}{Especially in the case of categorical variables and relatively small sample sizes, the conditional independence tests used by structural learning algorithms can have high type II error rates that lead to many false edge-removals and a severe information loss.}
As illustrated in the PHIA data analysis, very few features are connected by using the constrained-based structural learning algorithms. Results and discussions in detail can be found later in Section~\ref{section:Results}.

\textcolor{black}{In literature dealing with the unreliable conditional independence tests, most literature focuses on their negative effects in the orientation procedure but much less on those in the skeleton learning procedure, with a few exceptions. \citet{bromberg2009improving} proposed a method to resolve the inconsistencies in the conditional independence tests in the skeleton learning procedure. Their idea is to deduce a ``preference" score on the conditional independence test results from a certain set of axioms, such as Pearl's axioms \citep{dawid1979conditional, pearl1988probabilistic}, which all true conditional independence relationships should follow. These axioms can be seen as integrity constraints that can avoid certain inconsistent test outcomes.
However, the computational complexity of the algorithm proposed by \citet{bromberg2009improving} significantly increases as the number of vertices increases. Thus, the algorithm cannot be employed on the MPHIA dataset for which the graph is not sparse enough.
}

In this article, we propose a new causal structural learning algorithm that aims to preserve more information on important features and potential causal pathways. We apply the proposed algorithm on the Malawi PHIA \citep{ministry2017malawi} data set, which we refer to as MPHIA, and obtain interesting findings related to Tri90 pathways. We further compare and validate the proposed algorithm with some classical structural learning algorithms using information criteria and simulations

The remaining part of the paper is organized as follows. In Section~\ref{section:Preliminaries}, we provide definitions for important concepts in causal structural learning, such as d-separation, and also a basic summary for the graphical notations used in the paper.
In Section~\ref{section:Motivation}, we use a simple example to illustrate the problem of aggressive edge removal of existing structural learning algorithms. In Sections~\ref{section:MethodEquation} $\sim$ \ref{section:Overall}, we propose a new causal structural learning algorithm to overcome the aggressive edge removal issue.
In Section~\ref{section:Results}, we apply the proposed algorithm on the MPHIA data set, compare the results obtained from the proposed algorithm with those of the existing algorithms in Section~\ref{section:Validation}, and discuss the discovered Tri90 pathways in detail in Section~\ref{section:Discussions}.
In Section~\ref{section:Simulations},
we compare the numerical performance of the proposed algorithm with classical structural learning algorithms in simulation studies.
Section~\ref{section:Conclusions} provides a summary and discussion for the paper.
To save space, technical details, additional numerical results, and the relevant codebooks are provided in the Supplement.

\section{Method}
\label{section:Method}

In this section, we propose our causal structural learning algorithm.
We first provide definitions for important concepts in causal structural learning, such as d-separation, and also a basic summary for the graphical notations used in the paper in Section~\ref{section:Preliminaries}.
We then discuss the problem of aggressive edge removal of existing structural learning algorithms with a simplified example in Section~\ref{section:Motivation}.
To deal with this problem, we propose our new algorithm and provide a high-level overview of the algorithm in Section~\ref{section:MethodEquation}.
The proposed algorithm consists of two main steps, the forward step and the maximization step, which are explained in detail in Sections~\ref{section:Step1} and \ref{section:Step2}, respectively.
Section~\ref{section:Overall} provides the orientation procedure of our proposed algorithm and also summarizes the proposed method.


\subsection{Preliminaries}
\label{section:Preliminaries}


\textcolor{black}{
Suppose $\mathcal{G}$ is a directed acyclic graphical (DAG) model which represents the joint probability distribution over the vertex set $\tf{V}$ with directed edges and no directed loop. Each vertex in the graph represents a variable. We use variables $X$, $Y$, etc., to refer to the variables corresponding to the vertices $X$, $Y$, etc., and use the edge $X \rightarrow Y$ or $Y \leftarrow X$ to refer to the directed edge from $X$ to $Y$, where $X$ is a parent vertex of $Y$, and $Y$ is a child of $X$.
$X - Y$ denotes an undirected edge that could be either $X \rightarrow Y$ or $Y \rightarrow X$. We use a path to refer to an acyclic sequence of adjacent vertices and a causal path from $X$ to $Y$ to refer to a path that all arrows are pointing away from $X$ and into $Y$. If there is a causal path from $X$ to $Y$, we say that $X$ is an ancestor of $Y$ and that $Y$ is a descendant of $X$.}

\textcolor{black}{Next, we provide the formal definition of d-separation \citep{geiger1990logic, geiger1990d}.
A collider is a vertex on a path with two incoming arrows.
More specifically, a vertex $Z$ is a collider (v-structure) on a path $U$ if and only if the path $U$ contains a subpath $X \rightarrow Z \leftarrow Y$. 
For vertices $X$, $Y$ and a vertex set  $\tf{Z}$ which does not contain $X$ and $Y$, $X$ is d-connected to $Y$ given $\tf{Z}$ if and only if there is an acyclic path $U$ between $X$ and $Y$ such that every collider on $U$ is either a member of $\tf{Z}$ or an ancestor of a member of $\tf{Z}$, and no non-collider on $U$ is in $\tf{Z}$.
$X$ is d-separated from $Y$ given $\tf{Z}$ if and only if $X$ is not d-connected to $Y$ given $\tf{Z}$. As a simple illustration, suppose that $Z$ is the only vertex in $\tf{Z}$ and $U$. $X$ and $Y$ are d-connected given $Z$ for $X \rightarrow Z \leftarrow Y$; $X$ and $Y$ are d-separated given $Z$ for $X \rightarrow Z \rightarrow Y$, $X \leftarrow Z \leftarrow Y$, and $X \leftarrow Z \rightarrow Y$.} 


\textcolor{black}{Moreover, a set of variables $\tf{V}$ is causally sufficient if and only if no variable outside $\tf{V}$ is a direct cause of more than one variable in $\tf{V}$. For a causally sufficient set of variables $\tf{V}$ with probability distribution $P(\tf{V})$, the Markov condition assumes that the d-separation in the DAG $\mathcal{G}$ implies conditional independence in $P(\tf{V})$, i.e. if $X$ is d-separated from $Y$ by $\tf{Z}$ in $\mathcal{G}$, then $X$ is independent of $Y$ conditional on $\tf{Z}$ in $P(\tf{V})$; and the faithfulness assumption assumes that every conditional independence relationship in $P(\tf{V})$ is entailed by the d-separation relationship for the causal DAG $\mathcal{G}$, i.e. if $X$ is independent of $Y$ conditional on $\tf{Z}$, then $X$ is d-separated from $Y$ by $\tf{Z}$.
Therefore, the d-separation in the graph is equivalent to conditional independence in the distribution under the Markov condition and the faithfulness assumption. See \citet{spirtes2010introduction} for a more detailed introduction and discussion. In this paper, we always assume the causal sufficiency, the Markov condition, and the faithfulness assumption.}


\begin{comment}
\textcolor{black}{
We also make the assumption of
causal sufficiency. A domain is causally sufficient if it does not contain any hidden or latent variables.
Under the Markov condition, d-separation in the graph implies conditional independence in the distribution, i.e. if $X$ is d-separated from $Y$ by $\tf{Z}$, then $X$ is independent of $Y$ conditional on $\tf{Z}$.
Under the faithfulness assumption, the conditional independence implies d-separation in the graph, i.e. if $X$ is independent of $Y$ conditional on $\tf{Z}$, then $X$ is d-separated from $Y$ by $\tf{Z}$.
Therefore, d-separation is equivalent to conditional independence under the Markov condition and the faithfulness assumption. See \citet{spirtes2010introduction} for a more detailed introduction and discussion.
}
\end{comment}


\subsection{A Motivation Example}
\label{section:Motivation}


In the MPHIA data set, classical constrained-based structural learning algorithms discover very few important features for Tri90, as shown later in Section~\ref{section:Results}.
In this section, we explain the main cause of the problem and illustrate the motivation of our method with a simplified example from the MPHIA data. 
A complete analysis of the MPHIA data can be seen in Section~\ref{section:Results}.

\begin{example}
    \label{ex:ex1}
    \textcolor{black}{Suppose our graph contains only three vertices $X$, $Y$ and $Z$, and $X$ is Tri90Aware, the indicator variable for awareness of the HIV positive status, which is one of the variables in which we are mainly interested in our Tri90 goal study. 
   
   
    We want to check whether $Y$ (AlcoholFrequency, an ordinal variable for alcohol drinking frequency) and $Z$ (WealthQuintile, an ordinal variable for the wealthiness) are neighbors of $X$ (Tri90Aware). 
    That is to say, we are interested in whether the wealthiness and alcohol drinking frequency have direct causal relationships with the awareness of the HIV positive status.}
    We have the following four (conditional) independence test results among $X$, $Y$, and $Z$ for a confidence level $\alpha=0.05$ in the MPHIA dataset of males who are included in the Tri90 study.
   
    \begin{equation}
    \label{eq:testing}
    \begin{aligned}
    & X\perp Y \text{ is rejected, and } X\perp Z \text{ is rejected}, \\
    & X\perp Y| Z \text{ is not rejected, and } X\perp Z| Y \text{ is not rejected}.
    \end{aligned}
    \end{equation}
    
    \textcolor{black}{Since both $Y$ (AlcoholFrequency) and $Z$ (WealthQuintile) are not independent with $X$ (Tri90Aware) marginally from the testing results,
    they should be connected to $X$ (Tri90Aware) either directly or indirectly, according to the Markov condition.
    However, based on the faithfulness assumption,
    neither $Y$ (AlcoholFrequency) nor $Z$ (WealthQuintile) should be a neighbor of $X$ (Tri90Aware) since $X$ (Tri90Aware) is conditionally independent with $Y$ (AlcoholFrequency) given $Z$ (WealthQuintile) and $X$ (Tri90Aware) is conditionally independent with $Z$ (WealthQuintile) given $Y$ (AlcoholFrequency).}
    Those four testing results lead to contradicted conclusions under the Markov condition and the faithfulness assumption.
    Existing structural learning algorithms, such as the PC-stable algorithm, remove both the edge between Tri90Aware and AlcoholFrequency ($X - Y$) and the edge between Tri90Aware and WealthQuintile ($X - Z$) from the conditional independence testing results and conclude that there are no neighbors of $X$ (Tri90Aware), which might be too strict in detecting edges.

   
   
   
   
    
\end{example}

\textcolor{black}{Many covariates in the MPHIA data are categorical. When used as the conditional set, those categorical covariates lead to relatively high type II error rates for the conditional independence tests. In addition, triples $(X,Y,Z)$ with contradictory/inconsistent testing results such as Equation~(\ref{eq:testing}) are quite common in the MPHIA data. Such contradiction/inconsistency also leads to aggressive edge-removal for the constrained-based causal structural learning algorithms.}
To solve the false edge-removal issue, we propose a new graphical structural learning algorithm,
and we illustrate how the new algorithm successfully finds more edges later in Example~\ref{ex:ex2}.

\subsection{Method}
\label{section:MethodEquation}

In a directed acyclic graph (DAG) $\mathcal{G}$ with vertex set $\tf{V}$ and $X \in \tf{V}$.
Suppose $\tf{N}_X=(N_1,N_2,\cdots,N_q)$ is the parents and children set of $X$.
Under the faithfulness assumption and the Markov condition, respectively,
we have
\begin{eqnarray}
\label{eq:1}
&&X\not\perp N_i | \tf{S} \text{ for any } \tf{S}\subseteq(\tf{N}\backslash\{N_i\}), i=1,\cdots,q,\\
\label{eq:2}
&&X\perp M_i | \tf{S}_i \text{ for some } \tf{S}_i\subseteq\tf{N}_X, i=1,\cdots,r,
\end{eqnarray}
where
$
\tf{M} = \{M_1,M_2,\cdots,M_r\}
= \tf{V} \backslash (\tf{N}_X \cup \{X\})
$
is the set of vertices not connected to $X$.
Furthermore, under the Markov condition and the faithfulness assumption,
suppose the variable set, $\tf{N}$ satisfies the conditions \eqref{eq:1} and \eqref{eq:2},
then it is easy to verify that $\tf{N}$ is the set of parents and children of $X$ and thus $\tf{N} = \tf{N}_X$.

For each vertex $X$, we want to get the best subset of $\tf{V}$ that fulfills the conditions \eqref{eq:1} and \eqref{eq:2}.
Our procedure consists of two steps: a forward step and a maximization step.
The forward step finds all sets $\tf{N}$ that satisfy \eqref{eq:1} for vertex $X$ and will be illustrated in Algorithm~\ref{alg:step1} in Section~\ref{section:Step1}.
The maximization step picks the best set that fulfills \eqref{eq:2} among those sets found by Algorithm~\ref{alg:step1} and will be provided in Section~\ref{section:Step2}.


Before deriving the details of the algorithms, let us come back to Example~\ref{ex:ex1} and see how algorithms based on \eqref{eq:1} and \eqref{eq:2} can solve the problem of aggressive edge removal that presents in the existing classical structural learning algorithms.

\begin{example}
    \label{ex:ex2}
    \textcolor{black}{Same as Example~\ref{ex:ex1}, suppose our graph only contains three vertices $X$ (Tri90Aware), $Y$ (AlcoholFrequency), and $Z$ (WealthQuintile).}
   
    Further, assume that we have the same (conditional) independence test results as in Example~\ref{ex:ex1}.
    Although the four testing results in \eqref{eq:testing} are incompatible under the Markov condition and the faithfulness assumption,
    it is possible to have a valid graphical structure with only three of them.
    For example,
    \begin{equation}
    \label{eq:ex2(1)}
    \tikzfig{Example2a}
    \end{equation}
    and
    \begin{equation}
    \label{eq:ex2(2)}
    \tikzfig{Example2b}
    \end{equation}
    
    \textcolor{black}{Note that the neighborhood of $X$ (Tri90Aware) satisfies \eqref{eq:1} and \eqref{eq:2} for $X$.
    In the situations of \eqref{eq:ex2(1)}, the neighbor of $X$ (Tri90Aware) is $Z$ (WealthQuintile);
    and in the situations of \eqref{eq:ex2(2)}, the neighbor of $X$ (Tri90Aware) is $Y$ (AlcoholFrequency).
    However, $\{Y,Z\}$ does not satisfy \eqref{eq:1} and is not the neighborhood of $X$.
    Algorithm~\ref{alg:step1} finds both $\{Z\}$ and $\{Y\}$ as potential neighborhoods for $X$ (Tri90Aware) since both of them satisfy \eqref{eq:1} for $X$.
    Algorithm~\ref{alg:step2} compares between $\{Z\}$ and $\{Y\}$ and chooses the set which fulfills \eqref{eq:2} better as the neighborhood of $X$.
    In sum, the new algorithm chooses either $\{Z\}$ or $\{Y\}$ as the neighborhood of $X$, and concludes that either AlcoholFrequency ($Y$) or WealthQuintile ($Z$) is the neighbor of Tri90Aware ($X$) but not both of them.}
    So, the proposed algorithm is less aggressive in edge removal than the existing classical structural learning algorithms such as the PC-stable algorithm as discussed in Example~\ref{ex:ex1}.

\end{example}

\subsection{Forward Step}
\label{section:Step1}

Before stating Algorithm~\ref{alg:step1}, we first provide some useful definitions and their properties. Let $p$ be the total number of vertices.
For each vertex $X$, let $\tf{T} = \tf{V}\backslash \{X\} = \{T_i, i=1,2,\cdots,p-1\}$
and vertices are considered/added sequentially in the order of $T_1,T_2,\cdots,T_{p-1}$ to form a candidate neighborhood of $X$. That is to say, for any already formed non-empty candidate neighborhood $\tf{S} = \{T_{s_1},\cdots, T_{s_q}\} \subseteq \tf{T}$, where $1\leq s_1 < s_2 < \cdots < s_q \leq p-1$, we consider whether an additional vertex from $\tf{L}_X(\tf{S}):=\{T_{{s_q}+1},T_{s_q+2},\cdots,T_{p-1}\}$ can be added into the variable set, $\tf{S}$. We also define $\tf{L}_X(\emptyset)=\tf{T}$, which means that we need to consider all the vertices in $\tf{T}$ when $\tf{S}$ starts from an empty set.
Furthermore, let $\tf{C}_{X}(\tf{S})$ be the vertices in $\tf{L}_X(\tf{S})$ that can be added into $\tf{S}$ while still satisfying equation \eqref{eq:1}. That is to say,

\begin{comment}
\end{comment}
\begin{comment}
\begin{equation}
\tf{L}_X(\tf{S}) :=
\{T_{E(S)+1},T_{E(S)+2},\cdots,T_{p-1}\} =
\begin{cases}
\tf{T}, & q = 0 \\
\{T_{{s_q}+1},T_{s_q+2},\cdots,T_{p-1}\}, & q \geq 1
\end{cases}
\end{equation}
where $E_X(\tf{S})=0$ for $q=0$ and $E_X(\tf{S})=s_q$ for $q \geq 1$.
\end{comment}

\begin{equation}
\tf{C}_{X}(\tf{S}) := \{
T \in \tf{L}_X(\tf{S}) | \tf{N} = (\tf{S} \cup \{T\}) \text{ satisfies equation \eqref{eq:1}}
\}.
\end{equation}

Proposition~\ref{prop:cxs} establishes some useful properties of $\tf{C}_X(\tf{S})$, which are used in Algorithm~\ref{alg:step1} to facilitate the calculation of $\tf{C}_X(\tf{S})$.

\begin{proposition}
    \label{prop:cxs}

    $\tf{C}_X(\tf{S})$ has the following properties.
    \begin{enumerate}[topsep=0pt,itemsep=-1ex,partopsep=0ex,parsep=0ex]
        \item Let $\emptyset$ denote the empty set. If $\tf{S}=\emptyset$, then 
        \begin{equation}
        \label{eq:initialC}
        \tf{C}_X(\emptyset) = \{C|C\not\perp X\}.
        \end{equation}

        \item If $\tf{S}$ does not satisfy \eqref{eq:1}, then
        \begin{equation}
        \label{eq:emptyC}
        \tf{C}_X(\tf{S}) = \emptyset.
        \end{equation}

        \item If $\tf{S}_1\subseteq\tf{S}$, then
        \begin{equation}
        \label{eq:inclusion}
        \tf{C}_X(\tf{S}) \subseteq \tf{C}_X(\tf{S}_1).
        \end{equation}

        \item Let $\tf{S} = \{S_1,S_2,\cdots,S_n\}$, $n\geq 1$,
        $\tf{C}_X^{*}(\tf{S}) := \bigcap_{i=1}^n \tf{C}_X(\tf{S}_{-i})$,
        and $\tf{S}_{-i} = \tf{S} \backslash \{S_i\}$, for $i=1,2,\cdots,n$. It follows that
        \begin{equation}
        \label{eq:upperbound}
        \tf{C}_X(\tf{S}) \subseteq \tf{C}_X^{*}(\tf{S}).
        \end{equation}

        \item Let $\tf{S} = \{S_1,S_2,\cdots,S_n\}$, $n\geq 1$. If
        $\tf{S}$ satisfies \eqref{eq:1}, then
        \begin{equation}
        \label{eq:incremental_test}
        \tf{C}_X(\tf{S}) =
        \{C \in \tf{C}_X^{*}(\tf{S}) \cap \tf{L}_X(\tf{S}) :
        C\not\perp X|\tf{S},
        S_i\not\perp X| (\tf{S}_{-i} \cup \{C\}), i=1,\cdots,n
        \}.
        \end{equation}

    \end{enumerate}

\end{proposition}

The properties of $\tf{C}_X(\tf{S})$ in Proposition~\ref{prop:cxs} can be shown by using its definition,
and the proof can be found in Supplement~\ref{section:proof_cxs}.
Proposition~\ref{prop:cxs} establishes a recursive structure for $\tf{C}_X(\tf{S})$.
\eqref{eq:initialC} shows that $\tf{C}_X(\tf{S})$ consists of all vertices that are marginally correlated with $X$ when $\tf{S}$ is an empty set.
\eqref{eq:upperbound} states that $\tf{C}_X^{*}(\tf{S})$ is an upper-bound (concerning the partial order of inclusion) for $\tf{C}_X(\tf{S})$ when $\tf{S}$ is non-empty.
Furthermore, \eqref{eq:incremental_test} tells us
that we only need to examine the conditional independence relationship with conditional set of size $|\tf{S}|$ to get $\tf{C}_X(\tf{S})$ from $\tf{C}_X^{*}(\tf{S}) \cap \tf{L}_X(\tf{S})$ for non-empty $\tf{S}$.
\eqref{eq:incremental_test} is used in Algorithm~\ref{alg:step1} to reduce the number of conditional independence tests and speed up the finding process.
We can show that using \eqref{eq:incremental_test},
no conditional independence test is repeated in Algorithm~\ref{alg:step1}. Indeed,
let $s_{\max} = \max\{s_0,s_1,\cdots,s_q\}$, $S_{\max} = T_{s_{\max}}$, 
$\tf{s}_{-\max} = \{s_0,s_1,\cdots,s_q\}\backslash\{s_{\max}\}$,
and $\tf{S}_{-\max} = \{T_s| s\in\tf{s}_{-\max}\}$, 
then the conditional independence test of $X$ and $T_{s_0}$ given $\tf{S}=\{T_{s_1},T_{s_2}\cdots,T_{s_q}\}$ can only happen when we check whether $S_{\max}$ can be added into candidate set $\tf{S}_{-\max}$.

\begin{singlespace}
\begin{algorithm}[h]
    \caption{Forward Step: Find all potential parent and children sets of $X$ that satisfy  \eqref{eq:1}. The algorithm sequentially adds the candidate set satisfying \eqref{eq:1} into the preliminary result set, $\mathcal{R}$.}
    \label{alg:step1}
    \begin{algorithmic}
        \REQUIRE{
            (1) a vertex set $\textbf{V}$,
            (2) a vertex ${X}$ in $\tf{V}$, and
            (3) size $\alpha$ for conditional independence tests.
           
        }
        \ENSURE{The set $\mathcal{N}_X$ including all possible $\tf{N}$ that (approximately) satisfies equation \eqref{eq:1}.}
        \STATE{Let $\mathcal{S} = \{\emptyset\}$ and $\mathcal{R} = \emptyset$.} 
        \WHILE{$\mathcal{S}\neq\emptyset$}
        \STATE{Calculate $\tf{C}_X(\tf{S})$ for each $\tf{S}\in\mathcal{S}$ of cardinal from low to high as follows:}
        \IF{$\tf{S} = \emptyset$}
        \STATE{$\tf{C}_X(\emptyset) = \{C|C\not\perp X \text{ under size } \alpha\}$.}
        \ELSE
        \IF{$|\tf{S}| > M_{\text{CI}}$}
        \STATE{Set $\tf{C}_X(\tf{S}) = \tf{C}_X^{*}(\tf{S}) \cap \tf{L}_X(\tf{S})$. }
        \ELSE
        \STATE{Calculate $\tf{C}_X(\tf{S})$ from equation \eqref{eq:incremental_test}.}
        \ENDIF
        \ENDIF

        \STATE{Delete $\tf{S}$ from $\mathcal{S}$.}
        \IF{$\tf{C}_X(\tf{S})=\emptyset$}
        \STATE{Add $\tf{S}$ into the preliminary result set, $\mathcal{R}$.}
        \ELSE
        \STATE{Add every $\tf{S}\cup\{C\}$ for $C$ in $\tf{C}_X(\tf{S})$ into $\mathcal{S}$.}
        \ENDIF
        \ENDWHILE
        \RETURN{$\mathcal{N}_X = \{\tf{R}|\tf{R}\in\mathcal{R}, \text{ and \tf{R} is not a proper subset of any other sets in } \mathcal{R}\}$.}
    \end{algorithmic}
\end{algorithm}
\end{singlespace}

\noindent\tf{Remark 1}:
Note that the final output of Algorithm~\ref{alg:step1}, $\mathcal{N}_X$, only contains the candidate sets which are not proper subsets of other candidate sets for the following reasons. If $\tf{N}_X$ is a candidate set that satisfies \eqref{eq:1} for $X$, then all proper subsets of $\tf{N}_X$ satisfy \eqref{eq:1} for $X$, but no proper subset of $\tf{N}_X$ will satisfy \eqref{eq:2}. This is because
for any $\tf{N}_1\subsetneq \tf{N}_X$, take $M\in\tf{N}_X\backslash\tf{N}_1$, and we have $X\not\perp M|\tf{S}$ for any $\tf{S}\subseteq\tf{N}_1\subsetneq\tf{N}_X$ since $\tf{N}_X$ satisfies \eqref{eq:1}. Hence $\tf{N}_1$ does not satisfy \eqref{eq:2}.

\noindent\tf{Remark 2}: In the calculation of $\tf{C}_X(\tf{S})$ and $\tf{C}_X^{*}(\tf{S})$, the size of conditional sets in the conditional independence test is not restricted.
However, conditional independence tests with large conditional sets are not accurate, and the number of all possible conditional independence tests also grows exponentially with the size of conditional sets.
Hence we set the upper-bound of sizes of conditional sets in Algorithm~\ref{alg:step1},
which is also discussed by \citet{tsamardinos2006mmhc}
and is commonly implemented in causal structural learning software such as \citet{scutari2010bnlearn}.
In addition, as shown by \eqref{eq:incremental_test}, if we set the upper-bound of sizes of conditional sets to be $M_{\text{CI}}$, then we can approximate $\tf{C}_X(\tf{S})$ by
\begin{equation}
\tf{C}_X(\tf{S}) \approx
\tf{C}_X^{*}(\tf{S}) \cap \tf{L}_X(\tf{S}), \text{ if } |\tf{S}| > M_{\text{CI}}.
\end{equation}

\textcolor{black}{\noindent{\bf Remark 3}:
$M_{\text{CI}}$ is a tuning parameter in Algorithm 1.
Theoretically, one should choose $M_{\text{CI}}$ that is no less than the largest degree of vertices in the graph. However, there are concerns about using a large $M_{\text{CI}}$, and the true degrees are also unknown in practice.
First, a large $M_{\text{CI}}$ corresponds to big conditional sets, which will increase the computing cost. Second, when the sample size is not large enough, a large $M_{\text{CI}}$ leads to less reliable results due to the increased type II errors. 
Generally, it is recommended to set $M_{\text{CI}}=3$ when the true graph is expected to be sparse or moderately sparse \citep{yan2020effective}. We also carry on a simulation study on different choices of $M_{\text{CI}}$, and more details and discussions can be seen in Supplement~S.5.1.
}


\subsection{Maximization Step}
\label{section:Step2}

To choose the best neighbor set from the set of candidate sets $\mathcal{N}_X$, we check how well each candidate $\tf{N} \in \mathcal{N}_X$ satisfies equation \eqref{eq:2}.
Define
\begin{equation}
\label{eq:snxy}
S_N(X, Y) := \max_{S\subseteq N} \operatorname{CI}(X,Y|\tf{S}),
\end{equation}
where $\operatorname{CI}(X,Y|\tf{S})$ is the $p$-value of some chosen conditional independence test for $X$ and $Y$ given $\tf{S}$.
\textcolor{black}{
To see why we use $\max$ in \eqref{eq:snxy}, note that conditional independence relationship implies big $p$-values, and $S_N(X,Y)$ measures whether any subset of $\tf{N}$ makes $X$ and $Y$ conditional independent.
The maximum in \eqref{eq:snxy} is from the fact that \eqref{eq:2} only requires one $\tf{S}$ that makes $X$ and $Y$ conditionally independent given $\tf{S}$ as suggested by the large $p$-value.
The idea of using conditional independence set of the largest $p$-value in the causal structural learning algorithm can also be found in \citet{ramsey2016improving}.}

Define
\begin{equation}
\label{eq:qxn}
Q_X(\tf{N}) := \min_{M_i\in \tf{V} \backslash (\tf{N}\cup\{X\})} S_N(M_i,X),
\end{equation}
which measures how well it is for subsets of $\tf{N}$ to ``separate" $X$ from any vertices not in $\tf{N}\cup\{X\}$. \textcolor{black}{Here we are concerned about whether there is any $M_i\in \tf{V} \backslash (\tf{N}\cup\{X\})$ violates the conditional independence {between $M_i$ and $X$ given $\tf{N}$}, and hence the minimum of $p$-values across $M_i$'s is used. Equations~\eqref{eq:snxy} and \eqref{eq:qxn} together can be seen as a minimax procedure. If $\tf{N}$ satisfies Equation \eqref{eq:2}, then for any $M_i$ not in the neighborhood of $X$, there should exist $\tf{S}\subseteq\tf{N}$ such that $M_i\perp X|\tf{S}$ leading to a large value of $S_N(M_i, X)$. If $Q_X(\tf{N})$ is large, then for any $M_i$ not in the neighborhood of $X$, $S_N(M_i, X)$ is large and it is likely that $M_i\perp X|\tf{S}$. So $\tf{N}$ with larger $Q_X(\tf{N})$ is likely to be the true neighbor set of $X$.} Furthermore, we have the following useful properties of $S_N(X, Y)$ in Proposition~\ref{prop:snxy}, which are used in
Algorithm~\ref{alg:step2} to facilitate the calculation of $S_N(X, Y)$.

\begin{proposition}
    \label{prop:snxy}

    $S_N(X, Y)$ has the following properties.
    \begin{enumerate}[topsep=0pt,itemsep=-1ex,partopsep=0ex,parsep=0ex]
        \item Let $\emptyset$ be the empty set. Then
        \begin{equation}
        \label{eq:initialS}
        S_{\emptyset}(X, Y) = \operatorname{CI}(X,Y,\emptyset).
        \end{equation}

        \item If  $\tf{N}_1\subseteq\tf{N}$, then
        \begin{equation}
        \label{eq:inclusionS}
        S_N(X, Y) \geq S_{N_1}(X, Y).
        \end{equation}

        \item  Let $\tf{N} = \{N_1,N_2,\cdots,N_n\}$, $n\geq 1$,
        $S_N^{*}(X, Y) := \max_{i=1}^n S_{N_{-i}}(X, Y)$,
        and $\tf{N}_{-i} = \tf{N} \backslash \{N_i\}$, $i=1,2,\cdots,n$. Then
        \begin{equation}
        \label{eq:lowerboundS}
        S_N(X, Y) \geq S_N^{*}(X, Y).
        \end{equation}

        \item If $\tf{N}\neq\emptyset$, then
        \begin{equation}
        \label{eq:incremental_calcS}
        S_N(X, Y) = \max\{S_N^{*}(X, Y), \operatorname{CI}(X,Y|\tf{N})\}.
        \end{equation}

    \end{enumerate}

\end{proposition}

The properties of $S_N(X, Y)$ in Proposition~\ref{prop:snxy} can be shown by using its definition,
and the proof can be found in Supplement~\ref{section:proof_snxy}.
Similar to Proposition~\ref{prop:cxs}, Proposition~\ref{prop:snxy} shows a recursive structure in $S_N(X, Y)$.
\eqref{eq:initialS} shows how to calculate $S_N(X, Y)$ for an empty set.
\eqref{eq:lowerboundS} states that $S_N^{*}(X, Y)$ is a lower-bound for $S_N(X, Y)$ for a non-empty set $\tf{N}$.
Furthermore, \eqref{eq:incremental_calcS} tells us
that we only need to calculate the $p$-value of conditional independence test $ \operatorname{CI}(X,Y|\tf{N})$ to get $S_N(X, Y)$ from $S_N^{*}(X, Y)$ when $\tf{N}$ is non-empty.
\eqref{eq:incremental_calcS} is used in
Algorithm~\ref{alg:step2} to reduce the number of conditional independence tests and speed up the searching process.
It is easy to prove that using \eqref{eq:incremental_calcS},
no conditional independence test is repeated in
Algorithm~\ref{alg:step2} because the conditional independence test of $X$ and $Y$ given $\tf{N}$ only happens when we calculate $S_{N}(X,Y)$.

\noindent\tf{Remark}:
Similar to Algorithm~\ref{alg:step1}, we propose to set upper-bound for sizes of conditional sets to be $M_{\text{CI}}$,
then $S_N(X, Y)$ can be approximated by
\begin{equation}
\label{eq:appox_S}
S_N(X, Y) \approx \max_{S\subseteq N, |S| \leq M_{\text{CI}}} \operatorname{CI}(X,Y|\tf{S}).
\end{equation}
Furthermore, from \eqref{eq:incremental_calcS}, we can approximate $S_N(X, Y)$ by
\begin{equation}
\label{eq:incremental_calcS_approx}
S_N(X, Y) \approx
S_N^{*}(X, Y), \text{ if } |\tf{N}| > M_{\text{CI}}.
\end{equation}

Based on \eqref{eq:incremental_calcS} and \eqref{eq:incremental_calcS_approx}, we have Algorithm~\ref{alg:step2} which selects the $\tf{N}$ with the largest $Q_X(\tf{N})$ to be the neighbor set of $X$.

\begin{comment}
\begin{singlespace}
\begin{algorithm}[H]
    \caption{Calculation of $S_N(X, Y)$.}
    \label{alg:step2(1)}
    \begin{algorithmic}
        \REQUIRE{(1)
            a vertex set $\textbf{V}$,
            (2) a vertex ${X}$ in $\tf{V}$,
            (3) a vertex ${Y}$ in $\tf{V}$, and
            (4) vertex set $\tf{N} \subseteq \tf{V} \backslash \{X,Y\}$. \\ }
        \ENSURE{$S_N(X, Y)$.}
        \IF{$\tf{N} = \emptyset$}
        \STATE{$S_N(X, Y) = \operatorname{CI}(X,Y,\emptyset)$.}
        \ELSE
        \IF{$|\tf{N}| > M_{\text{CI}}$}
        \STATE{Set $S_N(X,Y) = S_N^{*}(X, Y)$.}
        \ELSE
        \STATE{Calculate $S_N(X, Y)$ from equation \eqref{eq:incremental_calcS}.}
        \ENDIF
        \ENDIF
        \RETURN{$S_N(X, Y)$.}
    \end{algorithmic}
\end{algorithm}
\end{singlespace}
\end{comment}

\begin{singlespace}
\begin{algorithm}[!h]
    \caption{Maximization Step: Find $\tf{N} \in \mathcal{N}_X$ with the largest $Q_X(\tf{N})$.}
    \label{alg:step2}
    \begin{algorithmic}
        \REQUIRE{
            (1) a vertex set $\textbf{V}$,
            (2) a vertex ${X}$ in $\tf{V}$, and
            (3) candidate neighbor set $\mathcal{N}_X$ of $X$.}
        \ENSURE{$\tf{N} \in \mathcal{N}_X$ with the largest $Q_X(\tf{N})$.}
        \STATE{Let $r = 0$ and $\tf{N}_X = \emptyset$.}
        \FOR{$\tf{N} \in \mathcal{N}_X$}
        \STATE{Let $\tf{M} = \tf{V} \backslash (\tf{N}\cup\{X\})$.}
        \STATE{Set $Q_X(\tf{N}) = \infty$.}
        \FOR{each vertex $M_i \in \tf{M}$}
       
       
        \IF{$\tf{N} = \emptyset$}
        \STATE{$S_N(M_i, X) = \operatorname{CI}(M_i,X,\emptyset)$.}
        \ELSE
        \STATE{Calculate $S_N(M_i, X)$ from equation \eqref{eq:incremental_calcS} or approximate $S_N(M_i, X)$ using \eqref{eq:incremental_calcS_approx}.}
       
       
       
       
       
        \ENDIF
       
        \IF{$S_N(M_i, X) \leq r$}
        \STATE{Continue the outer loop for next $\tf{N} \in \mathcal{N}_X$.}
        \ENDIF
        \STATE{Set $Q_X(\tf{N}) = \min\{Q_X(\tf{N}), S_N(M_i, X)\}$.}
        \ENDFOR
        \IF{$Q_X(\tf{N}) > r $}
        \STATE{Set $\tf{N}_X = \tf{N}$.}
        \STATE{Set $r = Q_X(\tf{N})$.}
        \ENDIF
        \ENDFOR
        \RETURN{$\tf{N} = \tf{N}_X$.}
    \end{algorithmic}
\end{algorithm}
\end{singlespace}

\subsection{Overall Structural Learning with Orientations}
\label{section:Overall}

With Algorithms~\ref{alg:step1} and \ref{alg:step2}, we can learn the neighbor set of every vertex $X$ in the DAG $\mathcal{G}$ and hence the skeleton of the DAG.
We can further orient edges according to the conditional independence relationship.
The overall structural learning algorithm is summarized in Algorithm~\ref{alg:structure_learn}. The orientation procedure is quite similar to those of classic constraint-based causal structural learning algorithms, such as the PC algorithm \citep{spirtes1991algorithm}. We also incorporate some prior knowledge into the orientations.
For example, we know covariates such as AgeGroup cannot be affected by other covariates such as Education, so if there is an edge between AgeGroup and Education, then we orient the edge as AgeGroup to Education.
If there are multiple d-separation sets $S(X, Y)$ for a non-adjacent pair $(X, Y)$, we shall use the d-separation set $S(X,Y)$ with the largest $p$-value to make the orientation results stable.
For $Z$ adjacent to both $X$ and $Y$, we call $X-Z-Y$ a v-structure and make the orientation $X\rightarrow Z\leftarrow Y$ if $Z\not\in S(X,Y)$.
Finally, note that there could be conflicting v-structures.
For example, if there is a chain $X-Z-Y-W$, and both $X-Z-Y$ and $Z-Y-W$ are v-structures, then the edge $Z-Y$ should be oriented as $Z\leftarrow Y$ from the v-structure of $X-Z-Y$ but $Z\rightarrow Y$ from the v-structure of $Z-Y-W$.
In Algorithm~\ref{alg:structure_learn}, we resolve the conflict between v-structures by comparing the $p$-value of the v-structure.
Going back to the previous example, if $S(X,Y)$ has a larger $p$-value than $S(Z,W)$, then we orient $Z-Y$ as $Z\leftarrow Y$ following v-structure $X-Z-Y$ and vice-versa.




\begin{singlespace}
\begin{algorithm}[t]
    \caption{Structural Learning Algorithm}
    \label{alg:structure_learn}
    \begin{algorithmic}
        \REQUIRE{
            (1) A vertex set $\textbf{V}$, and
           
            (2) size $\alpha$ for conditional independence tests.
           
        }
        \ENSURE{CPDAG $\mathcal{G}$.}
        \FOR{every vertex $X \in \tf{V}$}
        \STATE{Calculate the neighbor set $\tf{N}_X$ using Algorithms~\ref{alg:step1} and \ref{alg:step2}.}
        \ENDFOR
        \STATE{Start from a complete undirected graph $\mathcal{G}$ with the vertex set $\tf{V}$.}
        \FOR{every pair of vertices $(X,Y)$ connected in $\mathcal{G}$}
        \IF{$X\not\in \tf{N}_Y$ and $Y\not\in \tf{N}_X$}
        \STATE{Delete the edge between $X$ and $Y$ in $\mathcal{G}$.}
        \ENDIF
        \ENDFOR
        \FOR{each edge $X -Y$}
            \IF{there is prior knowledge on the orientation $X - Y$}
            \STATE{Orient $X - Y$ according to the prior knowledge.}
            \ENDIF
        \ENDFOR
        \FOR{each pair of non-adjacent variables $(X, Y)$ with a common neighbor $Z$}
        \STATE{Find d-separation set $S(X,Y)$ with the largest $p$-value in the neighborhood of $X$ and $Y$.}
        \IF{$Z \not\in S(X,Y)$}
        \STATE{Orient $X - Z - Y$ as $X \rightarrow Z \leftarrow Y$.}
        \ENDIF
        \ENDFOR
        \STATE{Form $\mathcal{G}$ by recursive orientation according to the following two rules:\\
            1. If $X - Y$ and there is a directed path from $X$ to $Y$, then orient $X - Y$ as $X \rightarrow Y$;\\
            2. If $X$ and $Y$ are not adjacent and there is a $Z$ such that $X \rightarrow Z$ and $Z - Y$, then orient $Z - Y$ as $Z \rightarrow Y$.}
    \end{algorithmic}
\end{algorithm}
\end{singlespace}

\noindent\textcolor{black}{\tf{Remark 1}: Note that multiple orientations may satisfy the inferred d-separation (conditional independence) structure. Hence, in Algorithm~3, we first use prior knowledge in the orientation process to establish orientations and to enhance the interpretability of the orientation result. For instance, we presume that {age may lead to education status} but not the other way around.
The prior knowledge is provided in the sample code of the supplemental material.}

\noindent\textcolor{black}{\tf{Remark 2}: Note that in the d-separation set searching procedure in Algorithm~3, we use the d-separation set $S(X,Y)$ with the largest $p$-value to stabilize the orientation results. 
It agrees with equation~\eqref{eq:snxy} used in the maximization step for skeleton learning.
The approach of using the conditional independence set with the largest $p$-value in the orientation process has also been used in \citet{ramsey2016improving}.
Sections~2.5 and 2.6 both use the largest $p$-value among the conditional independence tests but for different purposes. Section~2.5 is about skeleton learning, and Section~2.6 is about orientation. \citet{ramsey2016improving} only concerned the orientation phase but not the skeleton learning.}

\noindent\tf{Remark 3}:
To make the graph learned by the proposed algorithm more interpretable, we calculate $p$-values for the significance of the (undirectional) connections for all edges in the graph learned by Algorithm~\ref{alg:structure_learn}.
More specifically, for edge $X-Y$ between the covariates $X$ and $Y$, define $P(X,Y)$ to measure the significance of $X - Y$ as follows
\begin{equation}
\label{eq:PXY}
P(X,Y) = \min\{S_{N_X}(X,Y), S_{N_Y}(X,Y)\},
\end{equation}
where $\tf{N}_X$ and $\tf{N}_Y$ are the neighbors of $X$ and $Y$, respectively.
Note that the measure $P(X,Y)$ is undirectional, i.e.,
the orientation of the edge $X - Y$ has no effect on $P(X,Y)$ from its definition,
and $P(X,Y)$ also has no information on the edge orientation.



\section{Application to MPHIA Data}
\label{section:Results}

The MPHIA survey is a new HIV-focused, cross-sectional, household-based, nationally representative survey of adults and adolescents aged 15 years and older as well as children aged 0-14 years. In addition to HIV testing results, the survey contains demographic questions, such as age group, gender (preference to the gender collected in the adult or adolescent questionnaire), ethnic group, and HIV-related questions, such as access to preventive care and treatment services. \textcolor{black}{There were 26,871 survey participants and 1,407 covariates in total. We will focus on adults and adolescents aged 15 years and older with complete Tri90 related information. There are 2,217 such individuals in the MPHIA survey data, including 712 males and 1,505 females. That is to say, these 2,217 individuals are HIV-positive with known status of HIV awareness (Aware), antiretroviral therapy (ART), and viral load suppression (VLS).}

\textcolor{black}{In PHIA surveys, the 1,407 covariates include 34 continuous variables such as age and time, 32 discrete variables such as the number of partners, and 1,341 categorical variables (nominal ones such as gender and ethnic group, ordinal ones such as alcohol frequency). Some covariates are applicable to males or females only, i.e., the questions related to pregnancy are only applicable to females. Therefore, the numbers of vertices are quite different between the female graph and the male graph. 
After a data-preprocessing procedure (illustrated later), we will apply the causal structural learning algorithm to six datasets for each combination of gender and 90-90-90 goal separately.}



\input{tables/Table0}


\noindent{\bf Overview of the data-preprocessing procedure}: 
We first drop all the covariates with a dominant level (one level has $>99\%$ samples), because the MPHIA sample size is not large enough to reject any null hypothesis of the conditional independence involving those variables. Also, we ``merge" all the closely related covariates into a single covariate in the MPHIA dataset. For example, some multi-option questions are dummy-coded by many two-option ones in the MPHIA dataset,
and we combine them together to create multi-level categorical covariates. By ``merging" these kinds of covariates, we can reduce the number of covariates and improve the interpretability of our results.
Also, some covariates appear in the MPHIA dataset multiple times with the same meaning but different names, such as EthnicGroup and EthnicCode, and we keep one covariate and drop the others in such situations.
We further remove some covariates that are direct indicators of the Tri90 goals from the MPHIA data. For example, there are questions like whether the subject takes a certain ART medicine or not in the MPHIA survey. These covariates are strongly correlated to the ART status, but not helpful for the purpose of building the causal pathways because they may block the connection between the ART status and other meaningful covariates.
Furthermore, some categorical variables have many levels which can complicate the analysis. To reduce the number of levels while keeping the main information of each covariate, we keep the biggest levels of each covariate which cover at least 95\% of the individuals, and combine the remaining levels into the ``Others" category.
\textcolor{black}{We sub-sample the whole dataset by gender and by Tri90 goal to create six datasets, and within each of the six datasets, we further drop the covariates with a dominant level (one level has $>99\%$ samples). The sample sizes and the numbers of vertices in each dataset are summarized in Table~\ref{tab:Table0}. There is one covariate included in the female awareness dataset but not in the female ART and VLS datasets: LiveHere (whether the individual lives here or not).
The covariate is categorical and has a dominating level, which is below the threshold of 99\% in the female awareness dataset but is above the threshold in the female ART and VLS datasets. More details about the six datasets are provided in Supplement~\ref{section:6datasets}.}



\begin{comment}
\textcolor{black}{We notice that females and males can have different important covariates and causal pathways for Tri90. For example, pregnancy-related covariates play an important role in Tri90 pathways among females but not males. So, we apply the causal structural learning algorithm to females and males separately and identify the important covariates for each Tri90 goal. Within each of the six datasets (by gender and by Tri90 goal), we further drop the covariates with a dominant level (one level has $>99\%$ samples). Details about the six datasets are provided in Supplement~\ref{section:6datasets}.}
\end{comment}

In the next section, we compare the results obtained from the proposed algorithm with those of the existing algorithms in Section~\ref{section:Validation}. We further discuss the potential Tri90 pathways discovered by the proposed algorithm in detail in Section~\ref{section:Discussions}.


\subsection{Model Comparison}
\label{section:Validation}

We also apply the following existing structural learning algorithms to each dataset: 
\textbf{PC-stable} algorithm proposed by \citet{colombo2014order}, which is a stable/order-independent variant of the original PC (initials of the first names) algorithm proposed by \citet{spirtes1991algorithm};
\textbf{MMPC} (Max Min Parents and Children) algorithm proposed by \citet{tsamardinos2003mmpc}; 
\textbf{IAMB} (Incremental Association Markov Blanket) algorithm proposed by \citet{tsamardinos2003iamb}; 
\textbf{GS} (Grow-Shrink) algorithm proposed by \citet{margaritis2003gs}.
Notice that we use the PC-stable algorithm instead of the PC algorithm since the PC-stable algorithm is order-independent.
Order independence means that a random reordering of the variables does not affect the graphical learning result, which enhances the stability of the algorithm and also makes the results more interpretable.

\textcolor{black}{We use the existing graphical learning algorithms implemented in the \textbf{R} package \textbf{bnlearn} \citep{scutari2010bnlearn} in Sections~3 and 4.
The default conditional independence test statistic in the \textbf{bnlearn} package is the mutual information for categorical variables and the linear correlation for continuous variables.
We set the upper bound of sizes of conditional sets $M_{\text{CI}} = 3$ and the size of the conditional independence test $\alpha=0.05$ by default for all the causal structural learning algorithms used in the comparison.}


We first summarize the structural learning results by different graphical learning methods. 
\textcolor{black}{The results in Table~\ref{tab:Table1} shows that the proposed algorithm makes new discoveries regarding Tri90 goals.}
More specifically, Table~\ref{tab:Table1} show that the numbers of edges (NE) and the number of directed edges (NDE) are both much smaller than the number of vertices (NV) in the graphs learned by the existing algorithms including the PC-stable, the MMPC, the IAMB, and the GS algorithms. These existing algorithms lean toward fractured graphs and do not have much conditional independence information for the orientation. On the contrary, the proposed algorithm produces larger numbers of edges and well-connected graphs and has more potential to infer the directions of edges.

\begin{singlespace}
\input{tables/Table1.tex}
\end{singlespace}

\textcolor{black}{\noindent\tf{Remark}: It is important to note that the goal of the causal graphical algorithms is not to produce as many edges and directed edges as possible. Later in this section, we will use the Bayesian Information Criterion (BIC) to show that these discoveries made by the proposed algorithm provide useful information about the MPHIA data.
Furthermore, in Section~\ref{section:Discussions}, we discuss the discovered pathways in detail, which are reasonable, also confirmed in other HIV Tri90 literature, and can provide useful insight for the three Tri90 goals.}

Let $\mathbf{d}$ be the distance (defined as the length of shortest path regardless of direction) from a particular 90-90-90 goal (awareness of HIV, ART, or VLS) to a covariate and N($\mathbf{d} \leq k$) be the number of covariates whose distances to a 90-90-90 goal are smaller than or equal to $k$. Very few covariates are close to 90-90-90 goals in the graphs learned by the existing algorithms with only two covariates whose distance to the 90-90-90 goals are smaller or equal to three. On the contrary, the proposed algorithm discovers many covariates that are of a small distance to the 90-90-90 goals including several direct neighbors ($\mathbf{d} \leq 1$), which includes the ones discovered by the existing algorithms. See details in Supplement Table~\ref{tab:Table2}.




The Bayesian information criterion (BIC, \citep{BIC}) is a classical statistical tool for model selection.
We compare our proposed graphical learning algorithm with the aforementioned classical PC-stable, MMPC, IAMB, and GS algorithms using BIC criterion and summarize the results in Table~\ref{tab:IC}. Let $\tf{D}_k$, $k=1,2,\cdots,6$, be the data sets corresponding to the three 90-90-90 goals of each gender, respectively.
For $k=1,2,\cdots,6$,
we use each of the graphical learning algorithm $A_i$, for $i=1,2,\cdots,5$,
to learn a DAG $\mathcal{G}_{i,k}$ on $\tf{D}_{k}$,
where $A_i$, $i=1,2,\cdots,5$, stand for the PC-stable, MMPC, IAMB, GS, and the proposed algorithm, respectively.

The log-likelihood of a DAG $\mathcal{G}$ can be decomposed as follows:
\[
\ell(\theta|\mathcal{G})
=\frac{1}{n} \sum_{i=1}^n \log(p(x_{i1},\cdots,x_{ip}|\mathcal{G};\theta))
\]

\begin{equation}
= \sum_{j=1}^p \left[
    \frac{1}{n} \sum_{i=1}^n  \log(p(x_{ij}|\pi_{ij}(\mathcal{G});\theta_j)\right] 
= \sum_{j=1}^p \ell_j(\theta_j|\pi_{j}(\mathcal{G})),
\label{eqn:loglikelihood}
\end{equation}
where $\theta_j$'s are parameters of the model, $n$ is the sample size, $p$ is the number of covariates, $x_{ij}$, $1\leq i\leq n$, $1\leq j\leq p$, is the $i$-th observation of the $j$-th covariate $X_j$, $\pi_{j}(\mathcal{G})$ is the set of parents of $X_j$ in the DAG $\mathcal{G}$, and $\pi_{ij}(\mathcal{G})$ is the $i$-th observation of $\pi_{j}(\mathcal{G})$.
It shows that the log-likelihood of $\mathcal{G}$ can be decomposed as the sum of the log-likelihood of local structures of a covariate given its parents.
For the estimation of the local structures,
existing literature on causal graphical models often assumes a linear model of a covariate on its parents \citep{spirtes2010introduction, Valente2010searching, bolla2019graphical}.
Note that in the MPHIA data set, many covariates are categorical, so we assume generalized linear models (GLM) of a covariate given its parents instead of linear models, and we fit the GLM of local structures by MLE.
The degree of freedom (DF) and the log-likelihood of a DAG are the sums of DFs and log-likelihoods of the local structures, respectively.

Note that all the \textit{Log-likelihood*} in Table~\ref{tab:IC} are positive
since they are the differences between the log-likelihood and the log-likelihood of the null model (model with intercept only).
Furthermore, in Table~\ref{tab:IC},
BIC score is calculated by $-2\textit{Log-likelihood*} + \textit{DF} \log(n)$,
where $n$ is the sample size. Hence lower BIC scores correspond to better models.

\begin{singlespace}
\input{tables/IC.tex}
\end{singlespace}

Our proposed algorithm learns a much larger number of edges in all the six graphs compared with the existing algorithms in Table~\ref{tab:Table1}, and thus it has much larger degrees of freedom for the log-likelihood defined in Equation \eqref{eqn:loglikelihood}. Table~\ref{tab:IC} shows that the proposed algorithm has the largest degree of freedom
and the largest log-likelihood.
The much larger log-likelihoods imply that the proposed algorithm discovers a lot more useful information in MPHIA, which is further confirmed by the best (smallest) BIC scores in Table~\ref{tab:IC}.

\subsection{90-90-90 Pathways}
\label{section:Discussions}


Table~\ref{tab:Neighbors} lists the neighbors of the 90-90-90 goals discovered by the proposed causal structural learning algorithm.
Details of the graphs learned by the proposed algorithm are provided in Figures~\ref{fig:tri90awarefemale} and \ref{fig:tri90awaremale}|\ref{fig:tri90vlsmale}, which render parts of the graphs surrounding each 90-90-90 goal for each gender.
\textcolor{black}{The partial graphs we present in Figures~\ref{fig:tri90awarefemale} and \ref{fig:tri90awaremale}|\ref{fig:tri90vlsmale} are uniquely determined using our proposed algorithm with prior knowledge in the orientation process.
Hence the problem of multiple orientations does not affect the interpretation of the results of the MPHIA data. See Remark 1 of Algorithm~\ref{alg:structure_learn}.}
Also note that the covariates are renamed to have more intuitive meaning than their original names in MPHIA data codebooks. The MPHIA data codebook and the meaning of important covariates are provided in Supplement~\ref{section:Codebook}. Those covariates discovered by the proposed algorithm are consistent with findings in the existing literature, and we discuss them by Tri90 goals and by genders.
\textcolor{black}{
However, we should be cautious about the interpretation of the results of the proposed algorithm. Notably, the proposed algorithm assumes causal sufficiency, implying no hidden confounding (or common causes) covariates.
Although the MPHIA survey has included 1,407 related covariates, and most of the confounding covariates are likely observed, there may still be some unobserved confounding factors that can affect the interpretation of the proposed algorithm results.
}

\begin{singlespace}
\input{tables/Neighbors.tex}
\end{singlespace}

\subsubsection{HIV awareness among females}
For HIV awareness among females, we use Figure~\ref{fig:tri90awarefemale} to visually illustrate the covariates that are closely connected to HIV awareness (distances $\leq 2$).
The learned graphs for other Tri90 goals and genders are provided in Supplement~\ref{section:TableFigure} Figures~\ref{fig:tri90awaremale}|\ref{fig:tri90vlsmale}.


We find that HIV-positive females with more pregnant times are more likely to know their HIV status. We think this is because pregnant women are more likely to be tested for HIV during antenatal clinic visits or laboring. \citet{Peltzer2009} reported that age, condom usage, and education are associated with HIV status awareness. The proposed algorithm finds that AgeGroup and Education are potential ``reasons" for female's HIV awareness, and that EasyCondom and PLWHSupportGroup may be results of female's HIV awareness: (1) older individuals are more likely to know their HIV status; (2) individuals who are unaware of their HIV status often do not know whether it is easy to get condoms or have difficulties in getting condoms; (3) females unaware of their status never answer the PLWHSupportGroup question; (4) individuals with higher education are less likely to be aware of their HIV positive status, which is opposite of the marginal association reported in \citet{Peltzer2009}. It is because (4) is the conditional association given AgeGroup (which affects both Education and HIV awareness): while females with higher education levels in the age group 35 to 44 are more likely to be aware of HIV positive status, females with higher education levels in other age groups are less likely to be aware of the HIV status. Further studies will be needed to establish the relationship between education and HIV status awareness.



Figure~\ref{fig:tri90awarefemale} also presents other covariates that indirectly connect to the HIV awareness among female HIV patients. Among those, employment status (WorkLast12Mo) and marital status are impacted by both urban residence and education; and education is impacted by ethnic group, urban residence, and age group.
AgeGroup is a potential ``reason" for SellSexEver and SyphilisTestInPreg: 
younger generations are less likely to have SyphilisTestInPreg; older generations are less likely to have sold sex ever.
CircumcisedHIVRisk is a score measuring whether the individual agrees that men who are circumcised are not at risk of HIV at all, do not need to use condoms, and can have multiple sexual partners without the risk of HIV. The graph suggests 
females that are more cautious toward HIV risk (lower CircumcisedHIVRisk score) have smaller numbers of pregnancies on average. 
Finally, more sexual partners (PartNumber12Mo), being a partner to the head of house (RelationToHeadOfHouse), and lower tolerance towards violence (ViolenceOK) lead to easiness in getting condoms.


\begin{figure}[!ht]
    \begin{center}
    \begin{adjustbox}{center}
        \includegraphics[width=10in,height=6.6in]{"./DAGs/Tri90aware_female_bigfont.pdf"}

        \end{adjustbox}
       
       
        \begin{singlespace}
        \caption{90-90-90 Awareness graph in female.
        Vertices representing the Tri90 goals are biggest and marked by orange; vertices closer to goals have bigger sizes and darker colors than those farther away from goals. Widths of edges reflect the significance of the non-directional connection (conditional dependence) between vertices.
        Red and blue edges represent positive and negative relationships with Tri90 goals, respectively.
        Codebook can be found in Supplement~\ref{Codebook:tri90awarefemale}.}
        \label{fig:tri90awarefemale}
        \end{singlespace}
    \end{center}
\end{figure}


\subsubsection{HIV awareness among males}
We find males whose partners are older are more likely to be aware of their HIV status, both marginally, or conditioning on their age groups. Comparing the male awareness pathway and the female awareness pathway, we suspect that older females contribute more to HIV awareness because someone's age group is associated with his/her partner's age.

\begin{comment}
\citet{Peltzer2009} showed that the number of sexual partners is significantly associated with HIV status awareness.
\CC{More specifically, individuals with more sexual partners are more likely to be aware of their HIV status.}
\CC{\citet{Peltzer2009} did not make conclusions about the orientation, but I think it implies that more sexual partners lead to HIV status awareness.}
\lb{What is the direction of the association reported in the literature?}
\CC{WifeNumLiveElseWhere and PartnerNum1Mo share a common ``reason" WealthQuintile in Figure~\ref{fig:tri90awaremale} and are strongly positively correlated measured by Kendall's $\tau$.}
The proposed causal structural learning algorithm suggests that WifeNumLiveElsewhere and PartnerNumber12Mo be potential ``reasons" for male awareness and ART: (1) individuals with more than one partner in the last 12 months are less likely to be aware of their HIV status among HIV positive males and also less likely to be on ART treatment among HIV aware males \CC{conditioning on whether they have partners/wives live elsewhere}; and (2) individuals with one or more partners/wives live elsewhere are more likely to be on ART treatment among HIV aware males \CC{conditioning on whether they have more than one partners in the last 12 months}. \lb{I do not quite follow WifeNumLiveElsewhere.}
\CC{WifeNumLiveElsewhere means the number of wives or partners live elsewhere, so maybe we can rename it to PartnerNumLiveElsewhere?}
\end{comment}

In Figure~\ref{fig:tri90awaremale}, the proposed algorithm finds that PartnerNumber12Mo (sexual partners in the last 12 months) and WifeNumLiveElsewhere (number of wive/partners that live elsewhere) to be important covariates for male HIV awareness. These discoveries are supported by \citet{Peltzer2009} which showed that HIV awareness is significantly associated with the number of sexual partners. 
WifeNumLiveElsewhere and PartnerNum12Mo are strongly correlated and have closely related meanings.
Instead of studying their impacts on HIV awareness separately,
we study the relationship between HIV awareness and WealthQuintile which is the parent of both WifeNumLiveElsewhere and PartnerNum12Mo.
We find that ``poor" HIV positive males (with wealth in lower 20\%) are
more likely to have no partner in the last 12 months, more likely to have wives/partners live elsewhere, and less likely to be aware of their HIV status than the wealthier (in the upper 80\%) HIV positive males. This finding is consistent with \citet{dokubo2014awareness}.

We also find that ViolenceOK and AlcoholFrequency may be possible ``reasons" of HIV awareness in males, where ViolenceOK is a score to measure whether the individual believes it is right for a man to beat his wife/partner under various scenarios.
Further investigation shows that among males who are HIV positive: (1) those with a higher violence score are more likely to be unaware of their HIV status; (2) those who never drink are most likely to be aware of their HIV status while those with a high frequency of alcohol use (more than four drinks a week) are least likely to be aware of their HIV positive status, and such an effect is more significant in males who have not worked in the last 12 months. Alcohol usage is an important factor associated with sexual risk behavior \citep{Kalichman2007, Hahn2011}, and violence score can also be an important factor for sexual behavior. These factors deserve more attention when advocating HIV awareness.

The proposed causal structural learning algorithm also finds that WorkLast12Mo has a small distance to HIV awareness which is supported by \citet{Peltzer2009}. PLWHSupportGroup is a result of HIV awareness because males unaware of their status never answer the PLWHSupportGroup question.


\subsubsection{ART pathways}
Next, we discuss the ART pathways for females and males in one subsection because they are relatively simpler compared with the HIV awareness pathways. For the same reason, we discuss the VLS pathways for females and males in the same subsection.

\citet{Hodgson2014}, \citet{Tomori2014}, and \citet{Kebaabetswe2019} showed the connection between TravelTime to the receiving of ART. The proposed algorithm indicates that among people living with HIV and being aware of their positive HIV status, those who do not know the travel time to HIV care are less likely to receive ART treatment than the individuals who know the travel time.

\citet{Friedman2015, Hatcher2015, Sullivan2015} showed that violence and other kinds of abuse are significantly associated with ART initiation. We find females with smaller violence scores are more likely to be on ART treatment

\begin{comment}
The proposed algorithm suggests (1) having a large number of wives in the family (WifeNumOfHusband) may lead to ART initiation for females; (2) having a large number of sexual partners (PartenerNumber12Mo) likely causes ART initiation for males.
\CC{Individuals with more than one partner in the last 12 months are less likely to be on ART treatment among HIV aware males conditioning on whether they have partners/wives live elsewhere; and individuals with one or more partners/wives live elsewhere are more likely to be on ART treatment among HIV aware males conditioning on whether they have more than one partners in the last 12 months.}
\lb{Again, I am 100\% not sure how WifeNumLiveElsewhere is related to the number of sexual partners.}
\CC{
WifeNumLiveElsewhere actually means the number of wives or partners who live elsewhere, so maybe we can rename it to PartnerNumLiveElsewhere later?
Males with more sexual partners/wives live else are more likely to be on ART.}
\end{comment}


We find that wealthier males are more likely to initiate the ART treatment, which is in agreement with \citet{gebru2018perceived}. On the other hand, the males with ``poor" wealth are likely to have no partners in the last 12 months, and more likely to have wives or partners live elsewhere.
The proposed algorithm also discovers that females with known status of SyphilisTestInPreg are more likely to be on ART treatment. \textcolor{black}{The learned female ART pathway reveals that those  received ART treatments had easy access to the HIV care (short travel time) and easy access to antenatal care where they could be offered the syphilis test during pregnancy. They were both indicators of receiving good health care services which increased the chance of receiving ART treatment.}



Finally, we find AbnormPenisDischarge to be a potential ``reason" for male ART initiation. Further investigation shows that males with the missing response on abnormal penis discharge problems are more likely to be on ART treatment among the HIV-positive males who are aware of their HIV status. Unfortunately, it was unclear why people did not respond AbnormPenisDischarge question.



\subsubsection{VLS pathways}

We find females who attend the support groups more frequently are more likely to have their viral load suppressed. It is supported by \citet{Roberts2000} \citet{Tomori2014}, and \citet{Rangarajan2016}.


\citet{Friedman2015, Hatcher2015, Sullivan2015} showed that violence and other kinds of abuse are significantly associated with viral load suppression (VLS), especially among females. It supports the connection between VLS and ForceSexTimes found in the female VLS pathway. 



We also find SeekMedicalHelp and WifeNum to be potential ``reasons" of VLS among males who are on ART. Males who seek help from doctors or nurses because of health issues such as abnormal penis discharge and painful urination are more likely to have viral load unsuppressed; and males with more wives are more likely to have viral load suppressed.














\section{Simulation Studies}
\label{section:Simulations}

In Section~3, we see that the proposed causal structural learning algorithm discovers many more new edges than the existing ones, and we validate the results of the proposed algorithm through the Bayesian information criterion (BIC).
However, since the true causal graph for MPHIA data is unknown, we cannot validate the edge discoveries directly.
To get more insights on the Type I and Type II error rates and their trade-off for the proposed algorithm against existing ones, we carry out simulation studies with settings mimicking the MPHIA data in this section.

\input{newsimu1}

\textcolor{black}{We also present additional simulation studies in Supplement~\ref{section:AdditionalSimulation} to save space.}

\section{Conclusions}
\label{section:Conclusions}




UNAIDS 90-90-90 goals are important milestones to end AIDS.
To understand the progress on the 90-90-90 goals better, we analyze the Malawi PHIA (MPHIA) data set to discover important covariates and potential causal pathways for the 90-90-90 goals through causal structural learning in the paper.
Existing classical constrained-based causal structural learning algorithms are quite aggressive in edge removal and can lead to information losses while building directed graphical models, \textcolor{black}{especially in the case of categorical variables and relatively small sample sizes}.
To deal with the problem, we propose a new causal structural learning algorithm.
The proposed algorithm can preserve more information about important features and potential causal pathways as shown by various numerical studies \textcolor{black}{when many covariates in the domain are  categorical}.
More specifically, the proposed algorithm improves true positive rates over the existing classical algorithms while having a comparable true negative rate. It shows that our proposed algorithm has a great potential to discover important features and potential causal pathways, \textcolor{black}{especially in a domain with many categorical variables}.

Carrying out the causal graphical analysis on the MPHIA data set, we obtain interesting results on important covariates and possible causal pathways related to the UNAIDS 90-90-90 goals.
For example, the proposed algorithm discovers age and  condom usage to be important covariates for female HIV awareness and number of sexual partners to be important for male HIV awareness, which agrees with literature, such as \citet{Peltzer2009}.
The proposed algorithm also discovers similarities as well as differences between female and male pathways.
For example, travel time is discovered to be an important covariate for both female and male ART.
However, there are also different important covariates for female and male ART, such as attitude towards violence for female ART and partner numbers for male ART.

\textcolor{black}{It is also important to pay attention to the assumptions of the proposed causal structural learning algorithm when we interpret the results.
One important assumption behind the proposed algorithm is causal sufficiency which is critical for our algorithm and many other constrained-based causal structural learning algorithms. Although the MPHIA survey provides many related covariates, the causal sufficiency assumption may still not hold perfectly. }

\textcolor{black}{In the paper, we stratify the MPHIA dataset by sex and learn the graphical models for each sex and each Tri90 goal.
It is also possible to learn the DAG with further stratified data by other covariates such as age, but the sample size would be too small for each stratum to make reliable inferences. }


The discoveries on causality are important extensions for existing literature where only correlation instead of causation is established.
More studies can be carried out to further validate the potential causal discoveries, and other statistical inference tools such as mediation analysis can be applied to further the understanding of the causal relationship.
The discoveries on causality can help develop better HIV response strategies and related policies.





\begin{singlespace}












