\mathbf{e}gin{document}

\title{A Sparse Bayesian Approach to the Identification of Nonlinear State-Space Systems
\author{Wei Pan, Ye Yuan, Jorge Gon\c{c}alves and Guy-Bart Stan
	\thanks{W.~Pan and G.-B.~Stan are with the Centre for Synthetic Biology and
		Innovation and the Department of Bioengineering, Imperial College London, United Kingdom. Email: {\tt\small  \{w.pan11,g.stan\}@imperial.ac.uk.}
	Y.~Yuan and J. Gon\c{c}alves are with the Control Group, Department of Engineering, University of Cambridge, United Kingdom. J. Gon\c{c}alves is also with the Luxembourg Centre for Systems Biomedicine, Luxembourg. Email: {\tt\small  \{yy311,jmg77\}@cam.ac.uk.}
	Corresponding Author: Y.~Yuan and G.-B.~Stan.}
}
}
\maketitle

\mathbf{e}gin{abstract}
This technical note considers the identification of nonlinear discrete-time systems with additive process noise but without measurement noise. In particular, we propose a method and its associated algorithm to identify the system nonlinear functional forms and their associated parameters from a limited number of time-series data points. For this, we cast this identification problem as a sparse linear regression problem and take a Bayesian viewpoint to solve it. As such, this approach typically leads to nonconvex optimisations. We propose a convexification procedure relying on an efficient iterative re-weighted $\ell_1$-minimisation algorithm that uses general sparsity inducing priors on the parameters of the system and marginal likelihood maximisation. Using this approach, we also show how convex constraints on the parameters can be easily added to the proposed iterative re-weighted $\ell_1$-minimisation algorithm. In the supplementary material \cite{tacappendix}, we illustrate the effectiveness of the proposed identification method on two classical systems in biology and physics, namely, a genetic repressilator network and a large scale network of interconnected Kuramoto oscillators.
\end{abstract}
\mathbf{e}gin{IEEEkeywords}
Nonlinear System Identification, Sparse Bayesian Learning, Re-weighted $\ell_1$-Minimisation
\end{IEEEkeywords}

\section{INTRODUCTION}\label{sec:introduction}

Identification from time-series data of nonlinear discrete-time state-space systems with additive process noise is relevant to many different fields such as systems/synthetic biology, econometrics, finance, chemical engineering, social networks, etc. Yet, the development of general identification techniques remains challenging, especially due to the difficulty of adequately identifying nonlinear systems~\cite{ljung1999system,billings2013nonlinear}. Nonlinear dynamical system identification aims at recovering the set of nonlinear equations associated with the system from time-series observations. The importance of nonlinear dynamical system identification and its associated difficulties have been widely recognised \cite{billings2013nonlinear,sjoberg1995nonlinear}.

Since, typically, nonlinear functional forms can be expanded as sums of terms belonging to a family of parameterised functions (see \cite[Sec.~5.4]{ljung1999system} and \cite{billings2013nonlinear}), an usual approach to identify nonlinear state-space models is to search amongst a set of possible nonlinear terms (e.g., basis functions) for a parsimonious description coherent with the available data \cite{haber1990structure}. A few choices of basis functions are provided by classical functional decomposition methods such as Volterra expansion, Taylor polynomial expansion or Fourier series~\cite{ljung1999system,billings2013nonlinear,barahona1996detection}. This is typically used to model systems such as those described by Wiener and Volterra series \cite{barahona1996detection,wiener1966nonlinear}, neural networks \cite{narendra1990identification}, nonlinear auto-regressive with exogenous inputs (NARX) models \cite{leontaritis1985input}, and Hammerstein-Wiener \cite{bai1998optimal} structures, to name just a few examples.

Recently, graphical models have been proposed to capture the structure of nonlinear dynamical networks. In the standard graphical models where each state variable represents a node in the graph and is treated as a random variable, the nonlinear relations among nodes can be characterised by factorising the joint probability distribution according to a certain directed graph \cite{kollar2009probabilistic,pearl1988probabilistic,spirtes2000causation}. However, standard graphical models are often not adequate for dealing with times series directly. This is mainly due to two aspects inherent to the construction of graphical models. The first aspect pertains to the efficiency of graphical models built using time series data. In this case, the building of graphical models requires the estimation of conditional distributions with a large number of random variables \cite{barber2010graphical} (each time series is modelled as a finite sequence of random variables), which is typically not efficient. The second aspect pertains to the estimation of the moments of conditional distribution, which is very hard to do with a limited amount of data especially when the system to reconstruct is nonlinear. In the case of linear dynamical systems, the first two moments can sometimes be estimated from limited amount of data \cite{bach2004learning,materassi2012problem}. However, higher moments typically need to be estimated if the system under consideration is nonlinear.

In this technical note, we propose a method to alleviate the problems mentioned above. This method relies on the assumption that there exits a finite set of candidate dictionary functions whose linear combination allows to describe the dynamics of the system of interest. In particular, we focus on discrete-time nonlinear systems with additive noise represented in a general state-space form. Based on this, we develop an identification framework that uses time series data and \emph{a priori} knowledge of the type of system from which these time series data have been collected, e.g.,\ biological, biochemical, mechanical or electrical systems. For example in Genetic Regulatory Network (GRN), only polynomial or rational nonlinear functional forms typically need to be considered in the identification process.

To identify the network efficiently given the available time series data, we cast this nonlinear system identification problem as a sparse linear regression problem \cite{pan2012cdc,Candes2005decoding,donoho2006compressed}. Although such problems have been widely applied in the context of sparse coding, dictionary learning or image processing \cite{mairal2008discriminative,mairal2010online}, they have received little attention in nonlinear dynamical system identification. Besides the work presented here, one of the rare example of sparse estimation technique used for dynamical system identification is the multiple kernel-based regularisation method, which has been used to estimate finite impulse response models \cite{TianshiTAC}.

Furthermore, very few contributions are available in the literature that address the identification problem with \emph{a priori} information or constraints on the parameters of the system \cite{cerone2011enforcing,zavlanos2011inferring}. In contrast, our proposed framework allows us to incorporate convex constraints on the associated model parameters, e.g.,\ equality or inequality constraints imposed among parameters, or \emph{a priori} required stability conditions.

In sparse linear regression problems, finding the sparsest solution is desirable but typically NP-hard. The classic ``Lasso'' or $\ell_{1}$-minimisation algorithm are typically used as a relaxation to alleviate this numerical difficulty \cite{tibshirani1996regression}. However, these algorithms usually only work well or have performance guarantees when the considered dictionary matrix has certain properties such as the \textit{restricted isometry property} (RIP) \cite{Candes2005decoding,dai2009subspace} or the \emph{incoherence} property \cite{donoho2003optimally}. Loosely speaking, these properties require that the columns of the dictionary matrix are orthogonal, or nearly so. Unfortunately, such properties are hardly guaranteed for nonlinear identification problems and, as a consequence, $\ell_{1}$-relaxation based algorithms typically do not work well when these conditions are not satisfied.

In this technical note, we shall explain, from a probabilistic viewpoint, how a Bayesian approach can attenuate problems arising in the case of high correlations between columns of the dictionary matrix. In particular, the main contributions of this technical note are:
\mathbf{e}gin{itemize}
\item To formulate the problem of reconstructing discrete-time nonlinear systems with additive noise into a sparse linear regression problem. The model class in this technical note covers a large range of systems, e.g.,\ systems with multiple inputs and multiple outputs, systems with memory in their states and inputs, and autoregressive models.
\item To derive a sparse Bayesian formulation of the nonlinear system identification problem, which is casted into a nonconvex optimisation problem.
\item To develop an iterative re-weighted $\ell_1$-minimisation algorithm to convexify the nonconvex optimisation problem and solve it efficiently.
This formulation can also take into account additional convex constraints on the parameters of the model.
\end{itemize}

The generality of our framework allows it to be applied on a broad class of nonlinear system identification problems. In particular, to illustrate our results, we applied our approach to two examples: (1) the Genetic Repressilator Network, where we identify nonlinear regulation relationships between genes, transcriptional and translational strengths and degradation rates, and (2) a network of Kuramoto Oscillators, where we identify the network topology and nonlinear coupling functions. Details about these examples can be found in the supplementary material \cite{tacappendix}.

This technical note is organised as follows. Section \ref{sec:model} introduces the class of nonlinear models considered. Section \ref{sec:probstate} formulates the nonlinear identification problem into a sparse linear regression problem. Section~\ref{sec:bayesian} re-interprets the sparse problem from a Bayesian point of view, while Section~\ref{sec:cost1} shows how the resulting nonconvex optimisation problem can be convexified and solved efficiently using an iterative re-weighted $\ell_1$-minimisation algorithm. Finally, we conclude and discuss several future open questions.

\section{FORMULATION OF THE NONLINEAR IDENTIFICATION PROBLEM}\label{sec:identification}

\subsection{Considered Nonlinear Dynamical Model Class}\label{sec:model}

We consider dynamical systems described by discrete-time nonlinear state-space equations driven by additive Gaussian noise. The discrete-time dynamics of the $i$-th state variable $x_{i}$, $i=1,\ldots,n_{\mathbf{x}}$ is assumed to be described by:
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
x_i(t_{k+1})
& = \mathbf{F}_i(\mathbf{x}(t_{k}),\mathbf{u}(t_{k}))+\xi_i(t_{k})  \\
&=\sum_{i=1}^nolimits_{s=1}^{N_{i}}v_{is}f_{is}(\mathbf{x}(t_{k}),\mathbf{u}(t_{k}))+\xi_i(t_{k})  \\
&=\mathbf{f}^\top_i(\mathbf{x}(t_{k}),\mathbf{u}(t_{k}))\mathbf{v}_i+\xi_i(t_{k}),
\label{eq:expansion}
\end{aligned}
\end{equation}
where $\mathbf{x}=[x_{1},\ldots ,x_{n_\mathbf{x}}]^\top\in {\mathbb{R}}^{n_\mathbf{x}}$ denotes the state vector, $\mathbf{u}=[u_{1},\ldots ,u_{n_\mathbf{u}}]^\top\in {\mathbb{R}}^{n_\mathbf{u}}$ denotes the input vector, and $\mathbf{F}_i(\cdot): \mathbb{R}^{n_\mathbf{x}+n_\mathbf{u}}\rightarrow \mathbb{R}$ is a smooth nonlinear function which is assumed to be represented as a linear combination of several dictionary functions $f_{is}(\mathbf{x}(t_{k}),\mathbf{u}(t_{k})): \mathbb{R}^{n_\mathbf{x}+n_\mathbf{u}}\rightarrow \mathbb{R}$ (see Sec.~5.4 in \cite{ljung1999system}). These constituent dictionary functions can be monomial, polynomial, constant or any other functional form such as rational, exponential, trigonometric etc. $\mathbf{f}_i(\mathbf{x}(t_{k}),\mathbf{u}(t_{k}))$ is the vector of considered dictionary functions (which does not contain unknown parameters) while $\mathbf{v}_i \in {\mathbb{R}^{N_i}}$ appearing in \eqref{eq:expansion} is the weight vector associated with the dictionary functions vector. The additive noise $\xi_i(t_{k})$ is assumed to be i.i.d. Gaussian distributed with zero mean: $\xi_i(t_{k})\thicksim\mathcal{N}(0, Noise_i)$, with
$
\mathbb{E}(\xi_i (t_{p}))=0, \ \mathbb{E}(\xi_i(t_{p})\xi_i(t_{q}))=Noise_i\delta_{pq},
\label{model:noise}
$
where
$
\delta _{pq}=\left\{
\mathbf{e}gin{array}{ll}
1, & p=q, \\
0, & pNeq q
\end{array}
\right.
$.
$\xi_i(\cdot)$ and $\xi_j(\cdot)$ are assumed independent $\forall i Neq j$.

\mathbf{e}gin{remark}
The class of systems considered in (\ref{eq:expansion}) can be extended to the more general dynamics class
$
\mathbf{x}_i(t_{k+1})=\mathbf{F}_i(\mathbf{x}(t_{k}),\ldots, \mathbf{x}(t_{k-m_\mathbf{x}}), \mathbf{u}(t_{k}),\ldots,\mathbf{u}(t_{k-m_\mathbf{u}}) )+\mathbf{x}i(t_{k}),
$
where the ``orders'' $m_\mathbf{x}$ and $m_\mathbf{u}$ are assumed to be known \emph{a priori}, and $\mathbf{F}_i(\cdot): \mathbb{R}^{(m_{\mathbf{x}}+1)n_\mathbf{x}+(m_{\mathbf{u}}+1)n_\mathbf{u}}\rightarrow \mathbb{R}$. An example of such system can be found in the supplementary material \cite{tacappendix} (see Example 1). In particular, MIMO nonlinear autoregressive models belong to such descriptions.
\end{remark}

\subsection{Identification Problem Statement}
\label{sec:probstate}

If $M$ data samples satisfying (\ref{eq:expansion}) can be obtained from the system of interest, the system in (\ref{eq:expansion}) can be written as
$
\mathbf{y}_i=\mathbf{\Phi}c_i\mathbf{v}_i+\boldsymbol{\xi}_i, \ i=1,\ldots,n_{\mathbf{x}},
$
where
$\mathbf{y}_i \triangleq\left[x_i(t_{1}),\ldots,x_i(t_{M})\right]^\top\in {\mathbb{R}}^{M\times 1}$,
$\mathbf{v}_i \triangleq \left[v_{i1},\ldots,v_{iN_i}\right]^\top \in {\mathbb{R}}^{N_i\times 1}$,
$\boldsymbol{\xi}_i \triangleq
\left[\xi_i(t_{0}),\ldots,\xi_i(t_{M-1})\right]^\top\in {\mathbb{R}}^{M\times 1}$, and $\mathbf{\mathbf{\Phi}c}_i\in {\mathbb{R}}^{M\times N_i}$ represents the dictionary matrix with its $j$-th column being $[f_{ij}(\mathbf{x}(t_{0}),\mathbf{u}(t_{0})), \ldots, f_{ij}(\mathbf{x}(t_{M-1}),\mathbf{u}(t_{M-1}))]^{\top}$.

In this framework, the identification problem amounts to finding $\mathbf{v}_i \in \mathbb{R}^{N_i \times 1}$ given the measured data stored in $\mathbf{y}_i$. This, in turn, amounts to solving a linear regression problem, which can be done using standard least square approaches, provided that the structure of the nonlinearities in the model are known, i.e., provided that $\mathbf{\Phi}c_i$ is known. In what follows, we make the following assumption on the measurements contained in $\mathbf{y}_i$.
\mathbf{e}gin{assumption}
\label{assumption:observable}
The system (\ref{eq:expansion}) is fully measurable, i.e., time series data of all the state variables $x_i$ can be obtained.
\end{assumption}

Depending on the field for which the dynamical model needs to be built, only a few typical nonlinearities specific to this field need to be considered. In what follows we gather in a matrix $\mathbf{\Phi}_i$ similar to $\mathbf{\Phi}c_i$ the set of \emph{all} candidate/possible  dictionary functions that we want to consider for identification:
\mathbf{e}gin{equation}
\mathbf{y}_i=\mathbf{\Phi}_i\mathbf{w}_i+\boldsymbol{\xi}_i, \ i=1,\ldots,n_{\mathbf{x}}.
\label{problem:expand}
\end{equation}
The solution $\mathbf{w}_i$ to \eqref{problem:expand} is typically going to be sparse,  which is mainly due to the potential introduction of non-relevant and/or non-independent dictionary functions in $\mathbf{\Phi}_i$.

Since the $n_{\mathbf{x}}$ linear regression problems in (\ref{problem:expand}) are independent, for simplicity of notation, we omit the subscript $i$ used to index the state variable and simply write:
\mathbf{e}gin{equation}
\mathbf{y}=\mathbf{\Phi} \mathbf{w}+\boldsymbol{\xi}.
\label{problem}
\end{equation}
It should be noted that $N$, the number of dictionary functions or number of columns of the dictionary matrix $\mathbf{\Phi} \in \mathbb{R}^{M \times N}$, can be very large, at least larger than the number of observations $M$. Moreover, since $\mathbf{y}$ is constructed from time series data, typically two or more of the columns of the $\mathbf{\Phi}$ matrix are highly correlated. In this case standard methods, which involve some form of $\ell_1$-regularised minimisation, often yield poor performance on system identification~\cite{Candes2006robust}. 

\section{BAYESIAN VIEWPOINT ON THE RECONSTRUCTION PROBLEM}
\label{sec:bayesian}

\subsection{Sparsity Inducing Priors }
Bayesian modelling treats all unknowns as stochastic variables with certain probability distributions \cite{bishop2006pattern}. For $\mathbf{y}=\mathbf{\Phi} \mathbf{w}+\boldsymbol{\xi}$, it is assumed that the stochastic variables in the vector $\boldsymbol{\xi}$ are Gaussian i.i.d. with $\boldsymbol{\xi}\thicksim\mathcal{N}(\mathbf{0}, Noise\mathbf{I})$.
In such case, the likelihood of the data given $\mathbf{w}$ is
$
\mathcal{P}(\mathbf{y}|\mathbf{w})
=\mathcal{N}(\mathbf{y}|{\mathbf{\Phi}} {\mathbf{w}},Noise\mathbf{I})
\propto\exp \left[ -\frac{1}{2 Noise}\| \mathbf{y}-\mathbf{\Phi}\mathbf{w}\|_2^{2}\right].
$
We define a prior distribution $\mathcal{P}(\mathbf{w})$ as
$
\mathcal{P}(\mathbf{w})\propto\exp \left[-\frac{1}{2}\sum_{j}g(w_j)\right]=\prod_{j}\exp \left[-\frac{1}{2}g(w_j)\right]=\prod_{j}\mathcal{P}(w_j),
$
where $g(w_j)$ is a given function of $w_j$.
To enforce sparsity on $\mathbf{w}$, the function $g(\cdot)$ is usually chosen as a concave, non-decreasing function of $|w_j|$. Examples of such functions $g(\cdot)$ include Generalised Gaussian priors and Student's \emph{t} priors (see \cite{palmer2006variational} for details).

Computing the posterior mean $\mathbb{E}(\mathbf{w}|\mathbf{y})$ is typically intractable because the posterior $\mathcal{P}(\mathbf{w}|\mathbf{y})$ is highly coupled and non-Gaussian.
To alleviate this problem, ideally one would like to approximate $\mathcal{P}(\mathbf{w}|\mathbf{y})$ as a Gaussian distribution for which efficient algorithms to compute the posterior exist \cite{bishop2006pattern}. Another approach consists in considering \emph{super-Gaussian} priors, which yield a lower bound for the priors $\mathcal{P}(w_j)$ \cite{palmer2006variational}.
The sparsity inducing priors mentioned above are \emph{super-Gaussian}. More specifically, if we define $\boldsymbol{\gamma} \triangleq \left[\gamma_1, \ldots, \gamma_N\right]^\top \in \mathbb{R}^N_{+}$, we can represent the priors in the following relaxed (variational) form:
$
\mathcal{P}(\mathbf{w})=\prod_{j=1}^{n}\mathcal{P}(w_j)$, $\mathcal{P}(w_j)=\max_{\gamma_j>0}\mathcal{N}(w_j|0,\gamma_j)\varphi(\gamma_j)
$,
where $\varphi(\gamma_j)$ is a nonnegative function which is treated as a hyperprior with $\gamma_j$ being its associated hyperparameters. Throughout, we call $\varphi(\gamma_j)$ the ``\emph{potential function}''. This Gaussian relaxation is possible if and only if $\log \mathcal{P}(\sqrt{w_j})$ is concave on $(0,\infty)$.
The following proposition provides a justification for the above:

\mathbf{e}gin{proposition}\cite{palmer2006variational}
\label{super-Gaussian}
A probability density $\mathcal{P}(w_j)\equiv \exp(-g(w_j^2))$ can be represented in the convex variational form:
$
\mathcal{P}(w_j)=\max_{\gamma_j>0}\mathcal{N}(w_j|0,\gamma_j)\varphi(\gamma_j)
$
if and only if $-\log \mathcal{P}(\sqrt{w_j})=g(w_j)$ is concave on $(0,\infty)$. In this case  the potential function takes the following expression:
$
\varphi(\gamma_j)=\sqrt{{2\pi}/{\gamma_j}}\exp\left(g^{*}\left({\gamma_j}/{2}\right)\right)
$
where $g^{*}(\cdot)$ is the concave conjugate of $g(\cdot)$.
A symmetric probability density $\mathcal{P}(w_j)$ is said to be super-Gaussian if  $\mathcal{P}(\sqrt{w_j})$ is log-convex on $(0,\infty)$.
\end{proposition}

\subsection{Marginal Likelihood Maximisation}
For a fixed $\mathbf{g}amma=\left[{\gamma}_{1},\ldots, {\gamma}_{N}\right]$,
we define a relaxed prior, which is a joint probability distribution over $\mathbf{w}$ and ${\mathbf{g}amma}$, as
$
\mathcal{P}(\mathbf{w}; \mathbf{g}amma)
=\prod_{j}\mathcal{N}(w_j|0,{\gamma}_j)\varphi({\gamma}_j)
=\mathcal{P}(\mathbf{w}|\mathbf{g}amma)\mathcal{P}(\mathbf{g}amma) \le \mathcal{P}(\mathbf{w})
$,
where $\mathcal{P}(\mathbf{w}|\mathbf{g}amma)\triangleq \prod_{j}\mathcal{N}(w_j|0,{\gamma}_j),
\mathcal{P}(\mathbf{g}amma) \triangleq \prod_{j} \varphi({\gamma}_j)$.
Since the likelihood is  $\mathcal{P}(\mathbf{y}|\mathbf{w})$ is Gaussian, we can get a relaxed posterior which is also Gaussian
$
\mathcal{P}(\mathbf{w}|\mathbf{y},\mathbf{g}amma)
=\frac{\mathcal{P}(\mathbf{y}|\mathbf{w})\mathcal{P}(\mathbf{w};\mathbf{g}amma)}{\int \mathcal{P}(\mathbf{y}|\mathbf{w})\mathcal{P}(\mathbf{w};\mathbf{g}amma)d\mathbf{w}} 
=\mathcal{N}(\mathbf{m}_{\bw},\mathbf{\Sigma}_{\bw}).
$
Defining $\boldsymbol{\Gamma}amma\triangleq \operatorname{diag}[\mathbf{g}amma]$, the posterior mean and covariance are given by:
\mathbf{e}gin{eqnarray}
\mathbf{m}_{\bw}&=& \boldsymbol{\Gamma}amma \mathbf{\Phi}^\top (\noise\mathbf{I}+\dic\bGamma\dic^\top)^{-1} \mathbf{y}, \label{mean} \\
\mathbf{\Sigma}_{\bw}&=& \boldsymbol{\Gamma}amma- \boldsymbol{\Gamma}amma \mathbf{\Phi}^\top (\noise\mathbf{I}+\dic\bGamma\dic^\top)^{-1} \mathbf{\Phi}. \label{variance}
\end{eqnarray}
Now the key question is how to choose the most appropriate $\mathbf{g}amma=\hat{\mathbf{g}amma}=\left[\hat{\gamma}_{1},\ldots, \hat{\gamma}_{N}\right]$ to maximise $\prod_{j}\mathcal{N}(w_j|0,{\gamma}_j)\varphi({\gamma}_j)$ such that $\mathcal{P}(\mathbf{w}|\mathbf{y},\hat{\mathbf{g}amma})$ can be a ``good'' relaxation to $\mathcal{P}(\mathbf{w}|\mathbf{y})$.
Using the product rule for probabilities, we can write the full posterior as:
$
\mathcal{P}(\mathbf{w}, \mathbf{g}amma|\mathbf{y})
\propto \mathcal{P}(\mathbf{w}|\mathbf{y},\mathbf{g}amma)\mathcal{P}(\mathbf{g}amma|\mathbf{y})
= \mathcal{N}(\mathbf{m}_{\mathbf{w}},\mathbf{\Sigma}_{\mathbf{w}}) \times {\mathcal{P}(\mathbf{y}|\mathbf{g}amma)\mathcal{P}(\mathbf{g}amma)}/{\mathcal{P}(\mathbf{y})}.
$
Since $\mathcal{P}(\mathbf{y})$ is independent of $\mathbf{g}amma$, the quantity $\mathcal{P}(\mathbf{y}|\mathbf{g}amma)\mathcal{P}(\mathbf{g}amma)=\int \mathcal{P}(\mathbf{y}|\mathbf{w})\mathcal{P}(\mathbf{w}|\mathbf{g}amma)\mathcal{P}(\mathbf{g}amma)d\mathbf{w}$ is the prime target for variational methods \cite{wainwright2008graphical}. This quantity is  known as evidence or marginal likelihood. A good way of selecting $\hat{\mathbf{g}amma}$ is to choose it as the minimiser of the sum of the misaligned probability mass, e.g.,
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
\hat{\mathbf{g}amma}&=  \argmin\limits_{\mathbf{g}amma\geq\mathbf{0}} \int \mathcal{P}(\mathbf{y}|\mathbf{w})\left|\mathcal{P}(\mathbf{w})-\mathcal{P}(\mathbf{w};\mathbf{g}amma)\right|d\mathbf{w} \\
&= \argmax\limits_{\mathbf{g}amma\geq\mathbf{0}} \int \mathcal{P}(\mathbf{y}|\mathbf{w})\prod_{j=1}^{n}\mathcal{N}(w_j|0,\gamma_j)\varphi(\gamma_j)d\mathbf{w}.
\label{mass}
\end{aligned}
\end{equation}
The second equality is  a consequence of $\mathcal{P}(\mathbf{w};\mathbf{g}amma)\le\mathcal{P}(\mathbf{w})$. The procedure in (\ref{mass}) is referred to as  evidence maximisation or type-II maximum likelihood \cite{tipping2001sparse}. It means that the marginal likelihood can be maximised by selecting the most probable hyperparameters able to explain the observed data. Once $\hat{\mathbf{g}amma}$ is computed, an estimate of the unknown weights can be obtained by setting $\hat{\mathbf{w}}$ to the posterior mean (\ref{mean}) as $\hat{\mathbf{w}}=\mathbb{E}(\mathbf{w}|\mathbf{y};\hat{\mathbf{g}amma})=\hat{\boldsymbol{\Gamma}amma} \mathbf{\Phi}^\top (Noise\mathbf{I+{\mathbf{\Phi}}\hat{\boldsymbol{\Gamma}amma}}{{\mathbf{\Phi}}}^\top)^{-1} \mathbf{y}$, with $\hat{\boldsymbol{\Gamma}amma}\triangleq \operatorname{diag}[\hat{\mathbf{g}amma}]$. If an algorithm can be proposed to compute $\hat{\mathbf{g}amma}$ in (\ref{mass}), we can, based on it, obtain an estimation of the posterior mean $\hat{\mathbf{w}}$.

\subsection{Enforcing Additional Constraints on $\mathbf{w}$}

It is often important to be able to impose constraints on $\hat{\mathbf{w}}$ when formulating the optimisation problem (\ref{mass}) used to compute $\hat{\mathbf{w}}$ from $\hat{\mathbf{g}amma}$. In physical and biological systems, positivity of the parameters $\mathbf{w}$ of the system is an example of such constraints. Another example of constrained optimisation comes from stability considerations, which emerge naturally when the underlying system is known \emph{a priori} to be stable\footnote{Many stability conditions can be formulated as convex optimisation problems (see for example \cite{boyd1987linear,horn1990matrix}).}. Yet, only a few contributions in the literature address the problem of how to take into account \emph{a priori} information on system stability in the context of system identification \cite{cerone2011enforcing,zavlanos2011inferring}. To be able to integrate constraints on $\mathbf{w}$ into the problem formulation, we consider the following assumption on $\mathbf{w}$.

\mathbf{e}gin{assumption}
\label{assumption-constraints}
Constraints on the weights $\mathbf{w}$ can be described by a set of convex functions:
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
H^{[I]}_{i}(\mathbf{w})&\leq0,\ \ i=1,\ldots, m_I, \\
H^{[E]}_{j}(\mathbf{w})&=0,\ \ j=1,\ldots, m_E.
\label{convexconstraints}
\end{aligned}
\end{equation}
where the convex functions $H^{[I]}_{i}: \mathbb{R}^{N}\rightarrow \mathbb{R}$ are used to define inequality constraints, whereas the convex functions $H^{[E]}_{j}: \mathbb{R}^{N}\rightarrow \mathbb{R}$ are used to define equality constraints.

\end{assumption}

\section{NONCONVEX OPTIMISATION FOR IDENTIFICATION PROBLEMS}\label{sec:cost1}

In this section, we derive a sparse Bayesian formulation of the problem of system identification with convex constraints, which is casted into a nonconvex optimisation problem. The nonconvex optimisation problem can be dealt by an iterative re-weighted $\ell_1$-minimisation algorithm.

\subsection{Nonconvex Objective Function in Hyperparameter}
\label{sec:hyper}

\mathbf{e}gin{theorem}
\label{theorem:hyper}
The optimal hyperparameters $\hat{\mathbf{g}amma}$ in (\ref{mass}) can be obtained by minimising the following objective function
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
{\mathcal{L}_{\mathbf{g}amma}}\left( \mathbf{g}amma \right) &= \log \left\vert \noise\mathbf{I}+\dic\bGamma\dic^\top \right\vert \\
&+\mathbf{y}^\top(\noise\mathbf{I}+\dic\bGamma\dic^\top)^{-1}\mathbf{y}+\sum_{i=1}^nolimits_{j=1}^{N}p(\gamma_j),
\label{hypercostfunction}
\end{aligned}
\end{equation}
where $p(\gamma_j)=-2\log \varphi(\gamma_j)$.
The posterior mean is then given by
$\hat{\mathbf{w}}=\hat{\boldsymbol{\Gamma}amma} \mathbf{\Phi}^\top (Noise\mathbf{I+{\mathbf{\Phi}}\hat{\boldsymbol{\Gamma}amma}}{{\mathbf{\Phi}}}^\top)^{-1} \mathbf{y}$,
where $\hat{\boldsymbol{\Gamma}amma}=\operatorname{diag}[\hat{\mathbf{g}amma}]$.

\end{theorem}
\mathbf{e}gin{proof}
See Section A in the Appendix \cite{tacappendix}.
\end{proof}

\mathbf{e}gin{lemma}\label{lemma:concave}
The objective function in the hyperparameter $\mathbf{g}amma$-space, ${\mathcal{L}_{\mathbf{g}amma}}\left( \mathbf{g}amma \right)$ in (\ref{hypercostfunction}), is nonconvex.
\end{lemma}

\mathbf{e}gin{proof}
See Section B in the Appendix \cite{tacappendix}.
\end{proof}

\subsection{Nonconvex Objective Function in $\mathbf{w}$ with Convex Constraints }
\label{sec:with-constraint}

Based on the analysis in Section \ref{sec:hyper}, we first derive a dual objective function in the $\mathbf{w}$-space with convex constraints by considering the equivalent objective function of (\ref{hypercostfunction}) in the $\mathbf{g}amma$-space. We then show that this equivalent objective function is also nonconvex.

\mathbf{e}gin{theorem}
\label{theorem:MAP}
The estimate for $\mathbf{w}$ with constraints can be obtained by solving the optimisation problem
\mathbf{e}gin{equation}
\label{w-nonconvex}
\mathbf{e}gin{split}
\min_{\mathbf{w}} \|\mathbf{y}-\mathbf{\Phi}\mathbf{w}\|_{2}^{2}+Noise g_{\text{SB}}(\mathbf{w}), \,\,\,\,\, \mathrm{subject}\,\,\mathrm{to}  \,\,\,\,\
(\ref{convexconstraints})
\end{split}
\end{equation}
where
$
g_{\text{SB}}(\mathbf{w})=\min\limits_{\mathbf{g}amma\geq \mathbf{0}}\{\mathbf{w}^\top\boldsymbol{\Gamma}amma^{-1}\mathbf{w}+\log|\noise\mathbf{I}+\dic\bGamma\dic^\top|+\sum_{i=1}^nolimits_{j=1}^{N}p(\gamma_j)\} \label{penalty}
$
and the estimate of the stochastic variable $\mathbf{w}$ is given by the poseterior mean $\mathbf{m}_{\bw}$ defined in (\ref{mean}).
\end{theorem}

\mathbf{e}gin{proof}
See Section C in the Appendix \cite{tacappendix}.
\end{proof}

Although all the constraint functions are convex in Theorem \ref{theorem:MAP}, we show in the following Lemma that the objective function in \eqref{w-nonconvex} is nonconvex since it is the sum of convex and concave functions.
\mathbf{e}gin{lemma}
\label{lemma:gsb}
The penalty function $g_{\text{SB}}(\mathbf{w})$ in Theorem \ref{theorem:MAP} is a non-decreasing, concave function of $|\mathbf{w}|$ which promotes sparsity on the weights $\mathbf{w}$.
\end{lemma}
\mathbf{e}gin{proof}
The proof uses the duality lemma (see Sec. 4.2 in \cite{jordan1999introduction}).
See Section D in the Appendix \cite{tacappendix}.
\end{proof}

\subsection{Lasso Type Algorithm}
\label{sec:cost-constraint-convex}

We define the terms excluding $h^{*}(\mathbf{g}amma^{*})$ as
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
\mathcal{L}_{\mathbf{g}amma^{*}}(\mathbf{g}amma,\mathbf{w})\triangleq   Noiseinv\| \mathbf{y}-\mathbf{\Phi}\mathbf{w}\|_{2}^{2}+\sum_{i=1}^nolimits_j\left(\frac{w_j^2}{\gamma_j}+\gamma_j^{*}\gamma_j\right).
\end{aligned}
\end{equation}
For a fixed $\mathbf{g}amma^{*}$, we notice that $\mathcal{L}_{\mathbf{g}amma^{*}}(\mathbf{g}amma,\mathbf{w})$  is jointly convex in $\mathbf{w}$ and $\mathbf{g}amma$ and can be globally minimised by solving over $\mathbf{g}amma$ and then $\mathbf{w}$.
Since ${w_j^2}/{\gamma_j}+\gamma_j^{*}\gamma_j \geq 2w_j\sqrt{\mathbf{g}amma^{*}_{j}}$, for any $\mathbf{w}$,
$\gamma_j={|{w}_j|}/{\sqrt{\gamma^{*}_{j}}}$ minimises $\mathcal{L}_{\mathbf{g}amma^{*}}(\mathbf{g}amma,\mathbf{w})$.
When  $\gamma_j={|{w}_j|}/{\sqrt{\gamma^{*}_{j}}}$ is substituted into $\mathcal{L}_{\mathbf{g}amma^{*}}(\mathbf{g}amma,\mathbf{w})$, $\hat{\mathbf{w}}$ can be obtained by solving the following weighted convex $\ell_1$-minimisation procedure
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
\mathbf{\hat{w}}
=\argmin\limits_{\mathbf{w}}\left\{\| \mathbf{y}-\mathbf{\Phi}
\mathbf{w}\|_{2}^{2}+2Noise \sum_{i=1}^nolimits_{j=1}^{N}\sqrt{\mathbf{g}amma_{j}^{*}}|w_j| \right\}.
\label{wupdate}
\end{aligned}
\end{equation}
We can then set
$
\gamma_j=|\hat{w}_j|/\sqrt{\gamma^{*}_{j}}
\label{gammaupdate}
$, $\forall j$.
As a consequence, $\mathcal{L}_{\mathbf{g}amma^{*}}(\mathbf{g}amma,\mathbf{w})$ will be minimised for  any fixed $\mathbf{g}amma^{*}$.
Due to the concavity of $g_{\text{SB}}(\mathbf{w})$, the objective function in (\ref{w-nonconvex}) can be optimised using a re-weighted $\ell_1$-minimisation in a similar way as was considered in (\ref{wupdate}). The updated weight at the $k^{th}$ iteration is then given by
$
u_j^{(k)}\triangleq \left.\frac{\partial g_{\text{SB}}(\mathbf{w})}{2\partial |w_j|}\right|_{\mathbf{w}=\mathbf{w}^{(k)}}=\sqrt{\mathbf{g}amma^{*}_j}.
$

We can now explain how the update of the parameters can be performed based on the above. We start by setting the iteration count $k$ to zero and $u_j^{(0)}=1, \ \forall j$. At this stage, the solution is a typical $\ell_1$-minimisation solution. Then at the $k^{th}$ iteration, we initialise $u_j^{(k)}=\sqrt{\mathbf{g}amma^{*(k)}_j}$, $\forall j$ and then minimise over $\mathbf{g}amma$ using $\gamma_j=|w_j|/\sqrt{\gamma^{*}_{j}},\ \forall j$.
Consider again $\mathcal{L}_{\mathbf{g}amma,\mathbf{w}}(\mathbf{g}amma,\mathbf{w})$. For any fixed $\mathbf{g}amma$ and $\mathbf{w}$, {the tightest bound} can be obtained by minimising over $\mathbf{g}amma^{*}$. The tightest value of $\mathbf{g}amma^{*}=\hat{\mathbf{g}amma^{*}}$ equals the gradient of the function $h(\mathbf{g}amma)\triangleq \log|\noise\mathbf{I}+\dic\bGamma\dic^\top|+\sum_{i=1}^nolimits_{j=1}^{N}p(\gamma_j)$ defined in Lemma \ref{lemma:concave} at the current $\mathbf{g}amma$. $\mathbf{g}amma^{*}$ has the following analytical expression:
$
\hat{\mathbf{g}amma^{*}}
=Nabla_{\mathbf{g}amma}\left( \log |\noise\mathbf{I}+\dic\bGamma\dic^\top|+\sum_{i=1}^nolimits_{j=1}^{N}p(\gamma_j) \right)
= \operatorname{diag} \left[\mathbf{\Phi}^\top\left(\noise\mathbf{I}+\dic\bGamma\dic^\top\right)^{-1}{\mathbf{\Phi}} \right]+p'(\mathbf{g}amma),
$
where $p'(\mathbf{g}amma)=\left[p'(\gamma_1),\ldots,p'(\gamma_N)\right]^\top$. The optimal $\mathbf{g}amma^{*(k+1)}$ can then be obtained as
$
\mathbf{g}amma^{*(k+1)}= \operatorname{diag} \left[\mathbf{\Phi}^\top\left(Noise\mathbf{I}+\mathbf{\Phi} \mathbf{\Gamma}^{(k)} \mathbf{\Phi}^\top\right)^{-1}{\mathbf{\Phi}} \right]+p'(\mathbf{g}amma^{(k)}).
$
After computing the estimation of ${\gamma_j}^{(k)}={|w_j^{(k)}|}/{\sqrt{\gamma^{*(k)}_{j}}}$,
we can compute $\mathbf{g}amma^{*(k+1)}$, which gives
$
\gamma_j^{*(k+1)}= \mathbf{\Phi}_j^\top\left(Noise\mathbf{I}+\mathbf{\Phi} \mathbf{U}^{(k)} \mathbf{W}^{(k)} \mathbf{\Phi}^\top\right)^{-1}{\mathbf{\Phi}_j}+p^{\prime}(\gamma_j^{(k)}),
$
where
$
\mathbf{{\Gamma}}^{(k)}\triangleq \operatorname{diag}\left[{\mathbf{g}amma}^{(k)}\right]$, $
\mathbf{U}^{(k)}\triangleq \operatorname{diag}\left[\mathbf{u}^{(k)}\right]^{-1}=\operatorname{diag}\left[\sqrt{\mathbf{g}amma^{*(k)}}\right]^{-1} $,
$
\mathbf{W}^{(k)}\triangleq \operatorname{diag}\left[|\mathbf{w}^{(k)}|\right].
$
We can then define
$
u_j^{(k+1)}\triangleq \sqrt{\gamma_j^{*(k+1)}}
$
for the next iteration of the weighted $\ell_1$-minimisation. The above described procedure is summarised in Algorithm \ref{alg:weight}.
\mathbf{e}gin{algorithm}[!]

\caption{Nonlinear Identification Algorithm}
\label{alg:weight}
\mathbf{e}gin{algorithmic}[1]
	\State Collect time series data from the system of interest (assuming the system can be described by~\eqref{eq:expansion});
	\State Select the candidate dictionary functions that will be used to construct the dictionary matrix described in Section~\ref{sec:probstate};
	\State Initialise $u_j^0=1, \ \forall j$
	\For  {$k=0, \ldots, k_{\max}$}
	\State Solve the weighted $\ell_1$-minimisation problem with convex constraints on $\mathbf{w}$
	\mathbf{e}gin{equation}

	\mathbf{e}gin{split}
	\min_{\mathbf{w}} \| \mathbf{y}-\mathbf{\Phi}\mathbf{w}\|_{2}^{2}+2Noise\sum_{i=1}^nolimits_{j}u_j^{(k)}|w_j|, \,\,\,\,\, \mathrm{subject}\,\,\mathrm{to}  \,\,
	(\ref{convexconstraints}); Notag
	\end{split}
	\end{equation}
	\State Set $\mathbf{U}^{(k)}\triangleq \operatorname{diag}\left[\mathbf{u}^{(k)}\right]^{-1}, \
	\mathbf{W}^{(k)}\triangleq \operatorname{diag}\left[|\mathbf{w}^{(k)}|\right];$
	\State Update weights $u_j^{(k+1)}$ for the next iteration 
	$ u_j^{(k+1)}=\left[\mathbf{\Phi}_j^\top\left(Noise\mathbf{I}+\mathbf{\Phi} \mathbf{U}^{(k)} \mathbf{W}^{(k)} \mathbf{\Phi}^\top\right)^{-1}{\mathbf{\Phi}_j}+p^{\prime}(\gamma_j^{(k)})\right]^{{1}/{2}}$;
	
	\If   {a stopping criterion is satisfied}
	\State Break;
	\mathbb{E}ndIf
	\mathbb{E}ndFor
\end{algorithmic}
\end{algorithm}

\mathbf{e}gin{remark}
\label{remark:threshold}
There are two important aspects of the re-weighted $\ell_1$-minimisation algorithm presented in Algorithm~\ref{alg:weight}. First, for convex optimisation, there will be no exact zeros during the iterations and strictly speaking, we will always get a solution without any zero entry even when the RIP condition holds. However, some of the estimated weights will have very small magnitudes compared to those of other weights, e.g., $\pm 10^{-5}$ compared to $1$, or the ``energy'' some of the estimated weights will be several orders of magnitude lower than the average ``energy'', e.g., $\|w_j\|_2^2 \ll \|\mathbf{w}\|_2^2$. Thus a threshold needs to be defined \emph{a priori} to prune ``small'' weights at each iteration. The second aspect concerns the computational complexity of this approach. The repeated execution of Algorithm \ref{alg:weight} is very cheap computationally since it scales as $\mathcal{O}(MN\|\mathbf{w}^{(k)}\|_0)$ (see \cite{candes2008enhancing, wipf2010iterative}). Since at each iteration certain weights are estimated to be zero, certain dictionary functions spanning the corresponding columns of $\mathbf{\Phi}$ can be pruned out for the next iteration.
\end{remark}

\subsection{Convergence}

It is natural to investigate the convergence properties of this iterative re-weighted $\ell_1$-minimisation procedure. Let $\mathcal{A}(\cdot)$ denote a mapping that assigns to every point in $\mathbb{R}^{N}_{+}$ the subset of  $\mathbb{R}^{N}_{+}$ which satisfies Steps 5 and 6 in Algorithm \ref{alg:weight}. Then the convergence property can be established as follows:

\mathbf{e}gin{theorem}
\label{theorem:convergence}
Given the initial point $\mathbf{g}amma^{(0)}\in \mathbb{R}^{n}_{+}$ a sequence $\{\mathbf{g}amma^{(k)}\}_{k=0}^{\infty}$ is generated such that
$
\mathbf{g}amma^{(k+1)}\in \mathcal{A}(\mathbf{g}amma^{(k)}), \, \forall k.
$
This sequence is guaranteed to converge to a local minimum (or saddle point) of ${\mathcal{L}_{\mathbf{g}amma}}$ in (\ref{hypercostfunction}).
\end{theorem}
\mathbf{e}gin{proof}
The proof is in one-to-one correspondence with that of the Global Convergence Theorem \cite{zangwill1969nonlinear}.
See Section E in the Appendix \cite{tacappendix}.
\end{proof}

\section{Illustrative Numerical Examples}\label{sec:examples}

To implement Algorithm~\ref{alg:weight}, we use CVX, a popular package for specifying and solving convex programs \cite{grant2008cvx}.
To illustrate our results, the approach is applied to two classic examples: (1) the Genetic Repressilator Network, where we identify nonlinear regulation relationships between genes, transcriptional and translational strengths and degradation rates, and (2) a network of Kuramoto Oscillators, where we identify the network topology and nonlinear coupling functions. More details about these two examples and algorithmic comparisons with other algorithms described in \cite{sparselab} in terms of the Root of the Normalised Mean Square Error (RNMSE) and computational running time for different Signal-to-Noise Ratios (SNR) can be found in the supplementary material \cite{tacappendix}.
Importantly, this comparison shows that Algorithm~\ref{alg:weight} outperforms other classical algorithms \cite{sparselab} in terms of RNMSE, when used to identify the nonlinear systems associated with these illustrative examples.

\section{CONCLUSION AND DISCUSSION}\label{sec:conclusion}
This technical note proposed a new method for the identification of nonlinear discrete-time state-space systems with additive process noise. This method only required time-series data and some prior knowledge about the type of system from which these data have been acquired (e.g.,\ biochemical, mechanical or electrical). Based on this prior knowledge, candidate nonlinear functions (dictionary functions) can be selected for the particular type of system to be identified.

Due to the typical sparsity in terms of number of dictionary functions used to describe the dynamics of nonlinear systems and the fact that the number of measurements is typically small (at least smaller than the number of candidate nonlinear functions), the corresponding identification problem falls into the class of sparse linear regression problems. We considered this problem in a Bayesian framework and solved it efficiently using an iterative re-weighted $\ell_1$-minimisation algorithm. This approach also allowed us to easily add convex constraints from prior knowledge of some properties of the system (e.g.,\ positivity of certain variables, stability of the system, etc.). Finally, we illustrated how this approach can be efficiently used to accurately reconstruct discrete-time nonlinear models of the genetic repressilator and of Kuramoto networks.

Several important questions remain currently open for further research. Possibly, the most important is the assumption that the system is fully measurable. Typically, only part of the state is measured~\cite{yuan2011robust,yuan2012decentralised}, and, in particular, the number of hidden/unobservable nodes and their position in the network are usually unknown. We are currently investigating partial-measurement extensions of the method presented in this technical note.  Meanwhile, our algorithm is relatively more computationally expensive than other algorithms such as those in \cite{sparselab} but outperforms them all in terms of the accuracy of the identification as measured by the RNMSE. In future work, we plan to improve our proposed algorithm by exploiting further the structure of the optimisation problem at hand and reducing the associated algorithmic complexity. Another issue is that we assume that only process noise is present, and thus do not directly take into account measurement noise. We are currently working on an extension of the method allowing the incorporation of measurement noise into the presented framework.

\subsection{Proof of Theorem \ref{theorem:hyper}}
\label{proof:1}
We first re-express $\mathbf{m}_{\bw}$ and $\mathbf{\Sigma}_{\bw}$ in (\ref{mean}) and (\ref{variance}) using the Woodbury inversion identity:
\mathbf{e}gin{eqnarray}
\mathbf{m}_{\bw}&=& \boldsymbol{\Gamma}amma \mathbf{\Phi}^\top (\noise\mathbf{I}+\dic\bGamma\dic^\top)^{-1} \mathbf{y}=Noiseinv\mathbf{\Sigma}_{\bw}\mathbf{\Phi}^\top\mathbf{y}, \label{mean2} \\
\mathbf{\Sigma}_{\bw}&=& \boldsymbol{\Gamma}amma- \boldsymbol{\Gamma}amma \mathbf{\Phi}^\top (\noise\mathbf{I}+\dic\bGamma\dic^\top)^{-1} \mathbf{\Phi} \boldsymbol{\Gamma}amma=(\boldsymbol{\Gamma}amma^{-1}+Noiseinv\mathbf{\Phi}^\top \mathbf{\Phi})^{-1}. \label{variance2}
\end{eqnarray}
Since the data likelihood $\mathcal{P}(\mathbf{y}|\mathbf{w})$ is Gaussian, we can write the integral for the marginal likelihood in (\ref{mass}), as
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
&\int \mathcal{N}(\mathbf{y}|{\mathbf{\Phi}} {\mathbf{w}},Noise\mathbf{I})\prod_{j=1}^{N}\mathcal{N}(w_j|0,\gamma_j)\varphi(\gamma_j)d\mathbf{w} \\
=&\left(\frac{1}{2\piNoise}\right)^{M/2}\left(\frac{1}{2\pi}\right)^{N/2}
\int \exp\{-E(\mathbf{w})\}d\mathbf{w} \prod_{j=1}^{N}\frac{\varphi(\gamma_j)}{\sqrt{\gamma_j}},
\label{integral}
\end{aligned}
\end{equation}
where $E(\mathbf{w})=\frac{1}{2Noise}\|\mathbf{y}-\mathbf{\Phi}\mathbf{w}\|^2+\frac{1}{2}\mathbf{w}^\top\boldsymbol{\Gamma}amma^{-1}\mathbf{w}, \ \boldsymbol{\Gamma}amma=\operatorname{diag}(\mathbf{g}amma)$. Equivalently, we get
\mathbf{e}gin{eqnarray}
E(\mathbf{w})=\frac{1}{2}(\mathbf{w}-\mathbf{m}_{\bw})^\top\mathbf{\Sigma}_{\bw}^{-1}(\mathbf{w}-\mathbf{m}_{\bw})+E(\mathbf{y}),
\label{integral2}
\end{eqnarray}
where $\mathbf{m}_{\bw}$ and $\mathbf{\Sigma}_{\bw}$ are given by (\ref{mean2}) and (\ref{variance2}). Using the Woodbury inversion identity, we obtain:
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
E(\mathbf{y}) =\frac{1}{2}\left(Noiseinv\mathbf{y}^\top\mathbf{y}-Noiseinv\mathbf{y}^\top\mathbf{\Phi}\mathbf{\Sigma}_{\bw}\mathbf{\Sigma}_{\bw}^{-1}\mathbf{\Sigma}_{\bw}\mathbf{\Phi}^\top\mathbf{y}Noiseinv\right) =\frac{1}{2} \mathbf{y}^\top(\noise\mathbf{I}+\dic\bGamma\dic^\top)^{-1}\mathbf{y}.
\label{data-dependent-term0}
\end{aligned}
\end{equation}
Using (\ref{integral2}), we can evaluate the integral in (\ref{integral}) to obtain
$$\int \exp\{-E(\mathbf{w})\}d\mathbf{w}=\exp\{-E(\mathbf{y})\}(2\pi)^{N/2}|\mathbf{\Sigma}_{\bw}|^{1/2}.$$
Exploiting the determinant identity, we have $|\boldsymbol{\Gamma}amma^{-1}||\noise\mathbf{I}+\dic\bGamma\dic^\top | =|Noise\mathbf{I}| | \boldsymbol{\Gamma}amma^{-1}+Noiseinv\mathbf{\Phi}^\top \mathbf{\Phi} |$, 
from which we can compute the first term in (\ref{hypercostfunction}) as
$\log \left\vert \noise\mathbf{I}+\dic\bGamma\dic^\top \right\vert =-\log |\mathbf{\Sigma}_{\bw}|+M\logNoise+\log  |\boldsymbol{\Gamma}amma|$.
Then applying a $-2\log(\cdot)$ transformation to (\ref{integral}), we have
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
&-2\log \int \mathcal{P}(\mathbf{y}|\mathbf{w})\prod_{j=1}^{n}\mathcal{N}(w_j|0,\gamma_j)\varphi(\gamma_j)d\mathbf{w}  Notag \\
= &-\log |\mathbf{\Sigma}_{\bw}|+M\log2\pi Noise+\log  |\boldsymbol{\Gamma}amma|+\mathbf{y}^\top(\noise\mathbf{I}+\dic\bGamma\dic^\top)^{-1}\mathbf{y}+\sum_{i=1}^nolimits_{j=1}^{N}p(\gamma_j) Notag \\
= &  \log \left\vert \noise\mathbf{I}+\dic\bGamma\dic^\top \right\vert+M\log 2\pi Noise +\mathbf{y}^\top(\noise\mathbf{I}+\dic\bGamma\dic^\top)^{-1}\mathbf{y}+\sum_{i=1}^nolimits_{j=1}^{N}p(\gamma_j).
\label{integral3}
\end{aligned}
\end{equation}
From (\ref{mass}), we then obtain $\hat{\mathbf{g}amma} = \arg \min\limits_{\mathbf{g}amma\geq\mathbf{0}} \{\log \left\vert \noise\mathbf{I}+\dic\bGamma\dic^\top\right\vert +\mathbf{y}^\top(\noise\mathbf{I}+\dic\bGamma\dic^\top)^{-1}\mathbf{y}+\sum_{i=1}^nolimits_{j=1}^{N}p(\gamma_j)\}.$ We compute the posterior mean to get an estimate of $\mathbf{w}$: $\hat{\mathbf{w}}=\mathbb{E}(\mathbf{w}|\mathbf{y};\hat{\mathbf{g}amma})=\hat{\boldsymbol{\Gamma}amma} \mathbf{\Phi}^\top (Noise\mathbf{I+{\mathbf{\Phi}}\hat{\boldsymbol{\Gamma}amma}}{{\mathbf{\Phi}}}^\top)^{-1} \mathbf{y}$ where $\hat{\boldsymbol{\Gamma}amma}=\operatorname{diag}[\hat{\mathbf{g}amma}]$.

\subsection{Proof of Lemma \ref{lemma:concave}}
\label{proof:2}
We first show that the data-dependent term in (\ref{hypercostfunction}) is convex in $\mathbf{w}$ and $\mathbf{g}amma$. From (\ref{mean2}), (\ref{variance2}) and (\ref{data-dependent-term0}), the data-dependent term can be re-expressed as
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
&\mathbf{y}^\top\left(\noise\mathbf{I}+\dic\bGamma\dic^\top\right)^{-1}\mathbf{y}  \\
=& Noiseinv\mathbf{y}^\top\mathbf{y}-Noiseinv\mathbf{y}^\top\mathbf{\Phi}\mathbf{\Sigma}_{\bw} \mathbf{\Phi}^\topNoiseinv\mathbf{y} \\
=&  Noiseinv\|\mathbf{y}-\mathbf{\Phi}\mathbf{m}_{\bw}\|_2^2+ \mathbf{m}_{\bw}^\top\boldsymbol{\Gamma}amma^{-1}\mathbf{m}_{\bw}\\
=&\min_{\mathbf{x}} \{Noiseinv\| \mathbf{y}-\mathbf{\Phi}\mathbf{x}\|_{2}^{2}+\mathbf{x}^\top\boldsymbol{\Gamma}amma^{-1}\mathbf{x}\},
\label{data-dependent-term}
\end{aligned}
\end{equation}
where $\mathbf{m}_{\bw}$ is the posterior mean defined in (\ref{mean}). It can easily be shown that the minimisation problem is convex in $\mathbf{w}$ and $\mathbf{g}amma$, where  $\boldsymbol{\Gamma}amma\triangleq \operatorname{diag}[\mathbf{g}amma]$.

Next we define $h(\mathbf{g}amma)\triangleq \log|\noise\mathbf{I}+\dic\bGamma\dic^\top|+\sum_{i=1}^nolimits_{j=1}^{N}p(\gamma_j) \label{hr}$, and show $h(\mathbf{g}amma)$  is a concave function with respect to~$\mathbf{g}amma$. $\log|\cdot|$ is concave in the space of positive semi-definite matrices. Moreover, $\noise\mathbf{I}+\dic\bGamma\dic^\top$ is an affine function of $\mathbf{g}amma$ and is positive semidefinite for any $\mathbf{g}amma\geq 0$. This implies that $\log \left\vert \noise\mathbf{I}+\dic\bGamma\dic^\top \right\vert$ is a concave, nondecreasing function of $\mathbf{g}amma$.  Since we adopt a super-Gaussian prior with potential function $\varphi(\gamma_j), \forall j,$ as described in Proposition \ref{super-Gaussian}, a direct consequence is that $p(\gamma_j)=-\log\varphi(\gamma_j)$ is concave.

\subsection{Proof of Theorem \ref{theorem:MAP}}
\label{proof:3}
Using the data-dependent term in (\ref{data-dependent-term}), together with $\mathcal{L}_{\mathbf{g}amma}(\mathbf{g}amma)$  in (\ref{hypercostfunction}), we can create a strict upper bounding auxiliary function on $\mathcal{L}_{\mathbf{g}amma}(\mathbf{g}amma)$ as $\mathcal{L}_{\mathbf{g}amma,\mathbf{w}}(\mathbf{g}amma,\mathbf{w})= Noiseinv\| \mathbf{y}-\mathbf{\Phi}\mathbf{w}\|_{2}^{2}+\mathbf{w}^\top\boldsymbol{\Gamma}amma^{-1}\mathbf{w}+\log|\noise\mathbf{I}+\dic\bGamma\dic^\top|+\sum_{i=1}^nolimits_{j=1}^{N}p(\gamma_j)$. When we minimise over $\mathbf{g}amma$ instead of $\mathbf{w}$, we obtain
\mathbf{e}gin{eqnarray}
\mathcal{L}_{\mathbf{w}}(\mathbf{w})
&\triangleq& \min\limits_{\mathbf{g}amma\geq\mathbf{0}}\mathcal{L}_{\mathbf{g}amma,\mathbf{w}}(\mathbf{g}amma,\mathbf{w}) Notag \\
&=& Noiseinv\| \mathbf{y}-\mathbf{\Phi}\mathbf{w}\|_{2}^{2}+\min\limits_{\mathbf{g}amma\geq \mathbf{0}}\{\mathbf{w}^\top\boldsymbol{\Gamma}amma^{-1}\mathbf{w}+\log|\noise\mathbf{I}+\dic\bGamma\dic^\top|+\sum_{i=1}^nolimits_{j=1}^{N}p(\gamma_j)\}.
\end{eqnarray}
Then for $\mathbf{w}$ with convex constraints as described in Assumption \ref{assumption-constraints}, we obtain the formulation in Theorem~\ref{theorem:MAP}.

From the derivations in (\ref{data-dependent-term}), we can clearly see that the estimate of the stochastic variable $\mathbf{w}$ is the poseterior mean $\mathbf{m}_{\bw}$ defined in (\ref{mean}).

\subsection{Proof of Lemma \ref{lemma:gsb}}
\label{proof:4}

It is shown in Lemma \ref{lemma:concave} that $h(\mathbf{g}amma)$ is concave with respect to $\mathbf{g}amma\geq 0$. According to the duality lemma (see Sec. 4.2 in \cite{jordan1999introduction}), we can express the concave function $h(\mathbf{g}amma)$ as
$
h(\mathbf{g}amma)
=\min_{\mathbf{g}amma^{*}\geq 0}
\left<\mathbf{g}amma^{*}, \mathbf{g}amma\right>-h^{*}(\mathbf{g}amma^{*}),
$
where $h^{*}(\mathbf{g}amma^{*})$ is defined as the concave conjugate of
$h(\mathbf{g}amma)$ and is given by
$
h^{*}(\mathbf{g}amma^{*})
=\min_{\mathbf{g}amma\geq 0} \left<\mathbf{g}amma^{*}, \mathbf{g}amma\right>-h(\mathbf{g}amma).
$

From the proof of Lemma \ref{lemma:concave}, the data-dependent term $\mathbf{y}^\top\left(\noise\mathbf{I}+\dic\bGamma\dic^\top\right)^{-1}\mathbf{y}$ can be re-expressed as
$\min_{\mathbf{w}} \{Noiseinv\| \mathbf{y}-\mathbf{\Phi}\mathbf{w}\|_{2}^{2}+\mathbf{w}^\top\boldsymbol{\Gamma}amma^{-1}\mathbf{w}\}$.
Therefore we can create a strict upper bounding auxiliary function $\mathcal{L}_{\mathbf{g}amma,\mathbf{w}}(\mathbf{g}amma,\mathbf{w})$ on $\mathcal{L}_{\mathbf{g}amma}(\mathbf{g}amma)$ in (\ref{hypercostfunction}) by considering the fact that, in the dual expression,
$
\mathcal{L}_{\mathbf{g}amma,\mathbf{w}}(\mathbf{g}amma,\mathbf{w})
\triangleq \left<\mathbf{g}amma^{*}, \mathbf{g}amma\right>-h^{*}(\mathbf{g}amma^{*})+\mathbf{y}^\top\left(\noise\mathbf{I}+\dic\bGamma\dic^\top\right)^{-1}\mathbf{y}
= Noiseinv\| \mathbf{y}-\mathbf{\Phi}\mathbf{w}\|_{2}^{2}+\sum_{i=1}^nolimits_j\left({w_j^2}/{\gamma_j}+\gamma_j^{*}\gamma_j\right)-h^{*}(\mathbf{g}amma^{*}).
$
We can then re-express $g_{\text{SB}}(\mathbf{w})$ as
\mathbf{e}gin{equation}
g_{\text{SB}}(\mathbf{w})=\min\limits_{\mathbf{g}amma,\mathbf{g}amma{*}\geq \mathbf{0}}\left\{\sum_{i=1}^nolimits_j\left({w_j^2}/{\gamma_j}+\gamma_j^{*}\gamma_j\right)-h^{*}(\mathbf{g}amma^{*})\right\}.
\label{gsbl-dual1}
\end{equation}
$g_{\text{SB}}(\mathbf{w})$ is minimised over $\mathbf{g}amma$ when
$
\gamma_j=|w_j|/\sqrt{\gamma^{*}_{j}},\ \forall j.
$
Substituting this expression into $g_{\text{SB}}(\mathbf{w})$, we get
\mathbf{e}gin{equation}
g_{\text{SB}}(\mathbf{w})=\min\limits_{\mathbf{g}amma^{*}\geq \mathbf{0}}\left\{\sum_{i=1}^nolimits_j2\sqrt{\gamma^{*}_{j}}|w_j|-h^{*}(\mathbf{g}amma^{*})\right\}.
\label{gsbl-dual2}
\end{equation}
This indicates that $g_{\text{SB}}(\mathbf{w})$ can be represented as a minimum over upper-bounding hyperplanes in $\|\mathbf{w}\|_1$, and thus must be concave. $g_{\text{SB}}(\mathbf{w})$ thus promotes sparsity. Moreover, $g_{\text{SB}}(\mathbf{w})$ must be non-decreasing since $\mathbf{g}amma^{*}\geq \mathbf{0}$.

\subsection{Proof of Theorem \ref{theorem:convergence}}
\label{proof:convergence}

The proof is in one-to-one correspondence with that of the Global Convergence Theorem \cite{zangwill1969nonlinear}.
\mathbf{e}gin{enumerate}
\item The mapping $\mathcal{A}(\cdot)$ is compact. Since any element of $\mathbf{g}amma$ is bounded, ${\mathcal{L}}\left( \mathbf{g}amma\right)$ will not diverge to infinity.
In fact, for any fixed $\mathbf{y}$, $\mathbf{\Phi}$ and $\mathbf{g}amma$, there will always exist a radius $r$ such that for any $\|\mathbf{g}amma^{(0)}\|\leq 0$, $\|\mathbf{g}amma^{(k)}\|\leq 0$.
\item We denote $\mathbf{g}amma'$ as the non-minimising point of $\mathcal{L}(\mathbf{g}amma'')<\mathcal{L}(\mathbf{g}amma')$, $\forall$ $\mathbf{g}amma'' \in \mathcal{A}(\mathbf{g}amma')$. At any non-minimising $\mathbf{g}amma'$ the auxiliary objective function $\mathcal{L}_{(\mathbf{g}amma^{*})'}$ obtained from $\mathbf{g}amma^{*}_{\text{SB}}$ will be strictly tangent to $\mathcal{L}(\mathbf{g}amma)$ at $\mathbf{g}amma'$. It will therefore necessarily have a minimum elsewhere since the slope at $\mathbf{g}amma'$ is nonzero by definition. Moreover, because the $\log |\cdot|$ function is strictly concave, at this minimum the actual cost function will be reduced still further. Consequently, the proposed updates represent a valid descent function \cite{zangwill1969nonlinear}.
\item $\mathcal{A}(\cdot)$ is closed at all non-stationary points.
\end{enumerate}

Newpage

\onecolumn

\mathbf{e}gin{center}
	{\Large \textbf{Supplementary Material}}
\end{center}

\makeatletter

\renewcommand{S\arabic{section}}{S\arabic{section}}   
\renewcommand{S\arabic{subsection}}{S\arabic{subsection}}   
\renewcommand{S\arabic{table}}{S\arabic{table}}   
\renewcommand{S\arabic{figure}}{S\arabic{figure}}

\section{A Motivating Example}
\label{example:motivate}
In this example of a system with states and inputs memories, we will show how to construct the expanded dictionary matrix by adding candidate nonlinear functions.
\mathbf{e}gin{example}
\label{example:NARX3}
As an illustrative example, we consider  the following model of polynomial terms for a single-input single-output (SISO) nonlinear autoregressive system with exogenous input (NARX model) \cite{leontaritis1985input}:
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
x(t_{k+1})=0.7x^5(t_{k})x(t_{k-1})-0.5x(t_{k-2})
+0.6u^4(t_{k-2})-0.7x(t_{k-2})u^2(t_{k-1})
+\xi(t_k),
\label{NARX1}
\end{aligned}
\end{equation}
with $x, u, \xi \in {\mathbb{R}}$.
We can write \eqref{NARX1} in extended form as:
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
x(t_{k+1})&=w_1+w_2x(t_k)+\ldots+w_{m_x+2}x(t_{k-m_x})+\ldots+w_{N}x^{d_x}(t_{k-m_x})u^{d_u}(t_{k-m_u})+\xi(t_k) \\
&= \mathbf{w}^\top\mathbf{f}(x(t_{k}),\ldots, x(t_{k-m_x}), u(t_{k}),\ldots,u(t_{k-m_u}))+\xi(t_{k}),
\label{NARX3}
\end{aligned}
\end{equation}
where $d_x$ (resp. $d_u$) is the degree of the output (resp. input); $m_x$ (resp. $m_u$) is the maximal memory order of the output (resp. input); $\mathbf{w}^\top=[w_{1},\ldots,w_{N}]\in \mathbb{R}^{N}$ is the weight vector; and $\mathbf{f}(\mathbf{x}(t_{k}),\ldots, \mathbf{x}(t_{k-m_\mathbf{x}}), \mathbf{u}(t_{k}),\ldots,\mathbf{u}(t_{k-m_\mathbf{u}}))=[f_{1}(\cdot),
\ldots,f_{N}(\cdot)]^\top\in \mathbb{R}^{N}$ is the dictionary functions vector.
By identification of \eqref{NARX3} with the NARX model (\ref{NARX1}), we can easily see that $d_{x}=5$, $d_{u}=4$,  $m_{x}=2$, $m_{u}=2$.
To define the dictionary matrix, we consider all possible monomials up to degree $d_{x}=5$ (resp. $d_u=4$) and up to memory order $m_x=5$ (resp. $m_u=2$) in $x$ (resp. $u$).
This yields $\mathbf{f}(\cdot) \in \mathbb{R}^{1960}$ and thus $\mathbf{w} \in \mathbb{R}^{1960}$. Since $\mathbf{v} \in \mathbb{R}^{4}$, only 4 out of the 1960 associated weights $w_i$ are nonzero. 
\end{example}

\section{Illustrative Numerical Examples}\label{sec:examples}
\subsection{Experiment setup}

We hereafter present two classic examples that we use to illustrate our proposed method and on which we apply Algorithm~\ref{alg:weight}. 
To implement Algorithm~\ref{alg:weight}, we use CVX, a popular package for specifying and solving convex programs \cite{grant2008cvx}. The algorithm is implemented in MATLAB R2013a. The calculations were performed on a standard laptop computer (Intel Core i5 2.5GHz with 8GB RAM). In the examples that follow, we set the pruning threshold (mentioned in Remark~\ref{remark:threshold}) to $10^{-4}$, i.e., $\|w_j\|_2^2 / \|\mathbf{w}\|_2^2<10^{-4}$.

\subsection{Two classic examples}
\mathbf{e}gin{example}
\label{example:o}
In this example, we consider a classical dynamical system in systems/synthetic biology, the repressilator, which we use to illustrate the reconstruction problem at hand.
The repressilator is a synthetic three-gene regulatory network where the dynamics of mRNAs and proteins follow an oscillatory behaviour \cite{elowitz2000son}.
A discrete-time mathematical description of the repressilator, which includes both
transcription and translation dynamics, is given by the following set of discrete-time equations:
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
x_{1}(t_{k+1}) &=x_{1}(t_{k})+(t_{k+1}-t_{k})\left[-\gamma _{1}x_{1}(t_{k})+\frac{\alpha _{1}}{(1+x_{6}^{n_{1}}(t_{k}))}\right]+\xi_1(t_{k}),  Notag \\
x_{2}(t_{k+1}) &=x_{2}(t_{k}) +(t_{k+1}-t_{k})\left[-\gamma_{2}x_{2}(t_{k})+\frac{\alpha _{2}}{(1+x_{4}^{n_{2}}(t_{k}))}\right]+\xi_2(t_{k}),  Notag \\
x_{3}(t_{k+1})&=x_{3}(t_{k})+(t_{k+1}-t_{k})\left[-\gamma _{3}x_{3}(t_{k})+\frac{\alpha _{3}}{(1+x_{5}^{n_{3}}(t_{k}))}\right]+\xi_3(t_{k}),
Notag \\
x_{4}(t_{k+1}) &=x_{4}(t_{k}) +(t_{k+1}-t_{k})\left[-\gamma _{4}x_{4}(t_{k})+\mathbf{e}ta _{1}x_{1}(k)\right]+\xi_4(t_{k}),  Notag \\
x_{5}(t_{k+1})&=x_{5}(t_{k}) +(t_{k+1}-t_{k})\left[-\gamma _{5}x_{5}(k)+\mathbf{e}ta _{2}x_{2}\right]+\xi_5(t_{k}),  Notag \\
x_{6}(t_{k+1}) &=x_{6}(t_{k}) +(t_{k+1}-t_{k})\left[-\gamma _{6}x_{6}(t_{k})+\mathbf{e}ta _{3}x_{3}(t_{k})\right]+\xi_6(t_{k}).
\label{oscillator}
\end{aligned}
\end{equation}
Here, $x_{1},x_{2},x_{3}$ (resp. $x_{4},x_{5},x_{6}$) denote the concentrations of the mRNA transcripts (resp. proteins) of genes 1, 2, and 3, respectively.
$\xi_i$, $\forall i$ are i.i.d. Gaussian noise.
$\alpha _{1},\alpha
_{2},\alpha _{3}$ denote the maximum promoter strength for their corresponding gene,
$\gamma_{1},\gamma _{2},\gamma _{3}$ denote the mRNA degradation rates, $\gamma
_{4},\gamma _{5},\gamma _{6}$ denote the protein degradation rates, $\mathbf{e}ta
_{1},\mathbf{e}ta _{2},\mathbf{e}ta _{3}$ denote the protein production rates, and
$n_{1},n _{2},n _{3}$ the Hill coefficients.
The set of equations in (\ref{oscillator}) corresponds to a topology where gene $1$ is repressed by gene $2$,
gene $2$ is repressed by gene $3$, and gene $3$ is repressed by gene $1$.
Take gene 1 for example. The hill coefficient $n_{1}$ will typically have a value within a range from $1$ to $4$ due to biochemical constraints. The core question here is:
how can we determine the topology and kinetic parameters of the set of equations in (\ref{oscillator}) from time series data of $x_{1}, \dots, x_{6}$?

Note that we do not assume \emph{a priori} knowledge of the form of the nonlinear functions appearing on the right-hand side of the equations in \eqref{oscillator}, e.g., whether the degradation obeys first-order or enzymatic catalysed dynamics or whether the proteins are repressors or activators. It should also be noted that many linear and nonlinear functions can be used to describe the dynamics of GRNs in terms of biochemical kinetic laws, e.g., first-order functions $f(\left[ S\right] )=\left[ S\right]$, mass action functions $f(\left[ S_{1}\right] ,\left[ S_{2}\right] )=\left[ S_{1}\right]\cdot \left[ S_{2}\right]$,  Michaelis-Menten functions $f(\left[ S\right])=V_{\max }\left[ S\right] /(K_{M}+\left[ S\right] )$, or Hill functions $f(\left[ S\right] )=V_{\max }\left[ S\right] ^{n}/(K_{M}^{n}+[S]^{n})$.
These kinetic laws typical of biochemistry and GRN models will aid in the definition of the dictionary function matrix.
Next we show how the network construction problem of the repressilator model in \eqref{oscillator} can be formulated in a linear regression form.

Following the procedure described in Section II 
of the main text, we construct a candidate dictionary matrix $\mathbf{\Phi}$, by selecting as candidate basis functions, nonlinear functions typically used to represent terms appearing in biochemical kinetic laws of GRN models.
As a proof of concept, we only consider Hill functions as potential nonlinear candidate functions. The set of Hill functions with Hill coefficient $h$, both in activating and repressing from, for each of the $6$ state variables are:
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
\text{hill}_h(t_{k}) &\triangleq &\left[ \frac{1}{1+x_{1}^{h}(t_{k})},\ldots ,
\frac{1}{1+x_{6}^{h}(t_{k})}, \frac{x_{1}^{h}(t_{k})}{1+x_{1}^{h}(t_{k})},\ldots ,\frac{x_{6}^{h}(t_{k})}{1+x_{6}^{h}(t_{k})}\right]_{1\times 12}, 
\end{aligned}
\end{equation}
where $h$ represents the Hill coefficient. In what follows we consider
that the Hill coefficient can take any of the following integer values: $1$, $2$, $3$ or $4$. Since there are 6 state variables, we can construct the dictionary matrix $\mathbf{\Phi}$ with $6$ (dictionary functions for linear terms) $+(4*12)$ (dictionary functions for Hill functions) $=54$ columns.

\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
\mathbf{\Phi} =\left[
\mathbf{e}gin{array}{ccccccc}
x_1(t_0) & \ldots & x_6(t_0) & \text{hill}_1(t_0) & \ldots &
\text{hill}_4(t_0)  \\
\vdots &  & \vdots & \vdots &  & \vdots \\
x_1(t_{M-1}) & \ldots & x_6(t_{M-1}) & \text{hill}_1(t_{M-1}) & \ldots &
\text{hill}_4(t_{M-1})
\end{array}
\right] \in \mathbb{R}^{M\times (6+48)}.
\label{oscillator-dic}
\end{aligned}
\end{equation}

Then the output can be defined as
$$\mathbf{y}_i\triangleq
\left[\frac{x_i(t_{1})-x_i(t_{0})}{t_{1}-t_0},\ldots,\frac{x_i(t_{M})-x_i(t_{M-1})}{t_{M}-t_{M-1}}\right]^\top\in {\mathbb{R}}^{M\times 1}, i=1,\ldots,6.$$
Considering the dictionary matrix $\mathbf{\Phi}$ given in (\ref{oscillator-dic}), the corresponding target $\mathbf{w}_i$ for the ``correct'' model in (\ref{oscillator}) should be:
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
\mathbf{w}_{true} &= [\mathbf{w}_1, \mathbf{w}_2,\mathbf{w}_3,\mathbf{w}_4,\mathbf{w}_5, \mathbf{w}_6]  \\
&=
\left[
\mathbf{e}gin{array}{cccccc}
-\gamma_{1}(=-0.3) & 0 & 0 & \mathbf{e}ta _{1}(=1.4) & 0 & 0 \\
0 & -\gamma_{2}(=-0.4) & 0 & 0 & \mathbf{e}ta _{2}(=1.5) & 0 \\
0 & 0 & -\gamma_{3}(=-0.5)  & 0 & 0 & \mathbf{e}ta _{3}(=1.6) \\
0 & 0 & 0 & -\gamma_{4} (=-0.2) & 0 & 0 \\
0 & 0 & 0 & 0 & -\gamma_{5} (=-0.4) & 0 \\
0 & 0 & 0 & 0 & 0 & -\gamma_{6}(=-0.6)  \\
\mathbf{0}_{47\times 1} & \mathbf{0}_{45\times 1} & \mathbf{0}_{46\times 1} &  &  &  \\
\alpha _{1}(=4)  & \alpha _{2}(=3)  & \alpha _{3}(=5)  & \mathbf{0}_{48\times 1} & \mathbf{0}_{48\times 1} &
\mathbf{0}_{48\times 1} \\
\mathbf{0}_{0\times 1}  & \mathbf{0}_{2\times 1} & \mathbf{0}_{1\times 1} &  &  &  
\end{array}
\right] . 
\label{truepara}
\end{aligned}
\end{equation}
with values in brackets indicating the correct parameter values.

To generate the time-series data, we took `measurements' every $t_{k+1}-t_{k}=1$ between $t = 0$ and $t = 50$ (arbitrary units) from random initial conditions  which are drawn from a standard uniform distribution on the open interval $(0,1)$.
Thus a total of 51 measurements for each state are collected (including the initial value). It should be noted that the number of rows is less than the number of columns in the dictionary matrix.

\end{example}

\mathbf{e}gin{example}
\label{example:k}
A classical example in physics, engineering and biology is the Kuramoto oscillator network \cite{strogatz2000kuramoto}. We consider a network where the Kuramoto oscillators are nonidentical (each has its own natural oscillation frequency $\omega_i$) and the coupling strengths between nodes are not the same. The corresponding discrete-time dynamics can be described by
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
{\phi_i} (t_{k+1})={\phi_i} (t_{k})+(t_{k+1}-t_{k})\left[\omega_i+\sum\limits_{j=1, jNeq i}^{n} w_{ij}g_{ij}(\phi_j(t_k)-\phi_i(t_k))+\xi_i(t_k)\right], i=1, \ldots, n,
\label{Kuramoto}
\end{aligned}
\end{equation}
where $\phi_i\in[0,2\pi)$ is the phase of oscillator $i$, $\omega_i$ is its natural frequency, and the coupling function $g_{ij}$ is a continuous and smooth function, usually taken as $\sin, \ \forall i, j$. $w_{ij}$ represent the coupling strength between oscillators $i$ and $j$ thus $[w_{ij}]_{n\times n}$ defines the topology of the network.
Here, assuming we don't know the exact form of $g_{ij}$, we reconstruct from time-series data of the individual phases $\phi_i$ a dynamical network consisting of $n$ Kuramoto oscillators, i.e., we identify the coupling functions $g_{ij}(\cdot)$ as well as the model parameters, i.e., $\omega_i$ and $w_{ij}$, $i, j=1,\dots,n$.

To define the dictionary matrix $\mathbf{\Phi}$, we assume that all the dictionary functions are functions of a pair of  state variables only and consider $5$ candidate coupling functions $g_{ij}$:  $\sin(x_j-x_i)$, $\cos(x_j-x_i)$, $x_j-x_i$, $\sin^2 (x_j-x_i)$, and $\cos^2 (x_j-x_i)$. Based on this, we define the dictionary matrix as
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
\mathbf{\Phi}_{ij}(x_j(t_k),x_i(t_k)) \triangleq& [\sin(x_j(t_k)-x_i(t_k)), \cos(x_j(t_k)-x_i(t_k)), x_j(t_k)-x_i(t_k), Notag \\
& \sin^2 (x_j(t_k)-x_i(t_k)), \cos^2 (x_j(t_k)-x_i(t_k))] \in \mathbb{R}^{5}. Notag
\end{aligned}
\end{equation}

To also take into account the natural frequencies, we add to the last column of $\mathbf{\Phi}_i$ a unit vector. This leads to the following dictionary matrix $\mathbf{\Phi}_i$:
\mathbf{e}gin{equation}
\mathbf{e}gin{aligned}
\mathbf{\Phi}_{i} &\triangleq
\left[
\mathbf{e}gin{array}{cccc}
\mathbf{\Phi}_{i1}(x_1(t_0),x_i(t_0))& \ldots  &\mathbf{\Phi}_{in}(x_n(t_0),x_i(t_0)) & 1\\
\vdots  & \vdots  & \vdots  & \vdots \\
\mathbf{\Phi}_{i1}(x_1(t_{M-1}),x_i(t_{M-1})) & \ldots  & \mathbf{\Phi}_{in}(x_n(t_{M-1}),x_i(t_{M-1})) & 1
\end{array}
\right]
\in {\mathbb{R}}^{M\times {(5n+1)}}. Notag
\end{aligned}
\end{equation}
Then the output can be defined as
$$\mathbf{y}_i\triangleq
\left[\frac{\phi_i(t_{1})-\phi_i(t_{0})}{t_{1}-t_0},\ldots,\frac{\phi_i(t_{M})-\phi_i(t_{M-1})}{t_{M}-t_{M-1}}\right]^\top\in {\mathbb{R}}^{M\times 1}, i=1,\ldots,n.$$

To generate the time-series data, we simulated a Kuramoto network with $n=100$ oscillators, for which $10\%$ of the non-diagonal entries of the weight matrix $[w_{ij}]_{n\times n}$ are nonzero (assuming $g_{ii}$ and $w_{ii}$ are zeros), and the non-zero $w_{ij}$ values are drawn from a standard uniform distribution on the interval $[-10,10]$. The natural frequencies $\omega_{i}$ are drawn from a normal distribution with mean $0$ and variance $10$.
In order to create simulated data, we simulated the discrete-time model \eqref{Kuramoto} and took `measurements data points' every $t_{k+1}-t_{k}=0.1$ between $t = 0$ and $t = 45$ (in arbitrary units) from random initial conditions drawn from a standard uniform distribution on the open interval $(0,2\pi)$.
Thus a total of 451 measurements for each oscillator phase $\phi_i \in \mathbb{R}^{451\times 501}$ are collected (including the initial value).
Once again, it should be noted that the the number of rows of the dictionary matrix is less than that of columns.

\end{example}

\subsection{Algorithmic Performance Comparisons in terms of Signal-to-Noise Ratio}
Numberwithin{equation}{subsection}
We here investigate the performance of various algorithms including ours (Algorithm~\ref{alg:weight} in the main text)
for different signal-to-noise ratios of the data generated for Example \ref{example:o} and Example \ref{example:k}. We define the signal-to-noise ratio (SNR) as $\mathrm{SNR(dB)}\triangleq 20 \log_{10} (\|\mathbf{\Phi} \mathbf{w}_{\mathrm{true}}\|_2/\|\boldsymbol{\xi}\|_2)$. We considered SNRs ranging from 0 dB to 25 dB for each generated weight. To compare the reconstruction accuracy of the various algorithms considered, we use the root of normalised mean square error (RNMSE) as a performance index, i.e., $\|\hat{\mathbf{w}}-\mathbf{w}\|_2/\|\mathbf{w}\|_2$, where $\hat{\mathbf{w}}$ is the estimate of the true weight $\mathbf{w}$.
For each SNR, we performed 200 independent experiments and calculated the average RNMSE for each SNR over these 200 experiments.
In each ``experiment'' of Example \ref{example:o}, we simulated  the repressilator model with random initial conditions drawn from a standard uniform distribution on the open interval $(0,1)$. The parameters were drawn from a standard uniform distribution with the true values $\mathbf{w}_{true}$ in (\ref{truepara}) taken as the mean and variations around the mean values no more than $10\%$ of the true values. In MATLAB, one can use \texttt{$\mathbf{w}_{true}$.*(0.9 + 0.2*rand(54,6))} to generate the corresponding parameter matrix for each experiment.
In each ``experiment'' of Example \ref{example:k}, we simulated a Kuramoto network with $n=100$ oscillators, for which $10\%$ of the non-diagonal entries of the weight matrix $[w_{ij}]_{n\times n}$ were nonzero (assuming $g_{ii}$ and $w_{ii}$ are always zero). The non-zero $w_{ij}$ values were drawn from a standard uniform distribution on the interval $[-10,10]$. The natural frequencies $\omega_{i}$ were drawn from a normal distribution with mean $0$ and variance $10$.

Based on these settings, we compared Algorithm~\ref{alg:weight} with nine other state-of-the-art sparse linear regression algorithms available at \cite{sparselab}. \cite{sparselab} provides access to a free MATLAB software package managed by David Donoho and his team and contains various tools for finding sparse solutions of linear systems, least-squares with sparsity, various pursuit algorithms, and more. In Table I, we briefly describe the algorithms provided in \cite{sparselab} and used for this comparison.  In Fig. \ref{fig:nmseO} and Fig. \ref{fig:nmseK}, we plot, for various SNRs, the average RNMSE  obtained using our algorithm and other algorithms in \cite{sparselab} for the problems considered in Example \ref{example:o} and Example \ref{example:k} respectively.
In Fig. \ref{fig:timeO} and Fig. \ref{fig:timeK}, we plot, for the various SNRs considered, the average computational running time required by our algorithm and the other algorithms from \cite{sparselab} for the problems considered in Example \ref{example:o} and Example \ref{example:k} respectively. During this comparison, the inputs for the algorithms listed in Table I are always the same, i.e., the dictionary matrix $\mathbf{\Phi}$ and the data contained in $\mathbf{y}$. The initialisation and pre-specified parameters for these algorithms were set to their default values provided in \cite{sparselab}. Interested readers can download the package from \cite{sparselab} and reproduce the results presented here under the default settings of the solvers therein.

It should be noted that the dictionary matrices in all the experiments are rank deficient, i.e., neither column rank nor row rank are full. As a consequence, both the MP and OMP algorithm fail to converge or yield results with extremely large RNMSE. As these two algorithms cannot satisfactorily be used, they have been removed from the comparison results presented in Figures~\ref{fig:nmseO} to~\ref{fig:timeK}.

\mathbf{e}gin{center}
\mathbf{e}gin{table}
\mathbf{e}gin{tabular}{ | p{2.5cm} |  p{3.5cm} | p{9.5cm} |}
\hline
Abbreviation & Solver name in \cite{sparselab} & Method  \\ \hline
BP         &     \verb"SolveBP.m"          &   Basis Pursuit \\ \hline
IRWLS        &     \verb"SolveIRWLS.m"      &    Iteratively ReWeighted Least Squares\\ \hline

ISTBlock       &     \verb"SolveISTBlock.m"       &   Iterative Soft Thresholding, block variant with least squares projection\\ \hline
LARS      &     \verb"SolveLasso.m"         &   Implements the LARS algorithm\\ \hline
MP      &     \verb"SolveMP.m"             &    Matching Pursuit\\ \hline
OMP      &     \verb"SolveOMP.m"            &   Orthogonal Matching Pursuit\\ \hline
PFP      &     \verb"SolvePFP.m"            &   Polytope Faces Pursuit algorithm\\ \hline
Stepwise      &     \verb"SolveStepwise.m"       &    Forward Stepwise\\ \hline
StOMP       &     \verb"SolveStOMP.m"          &   Stagewise Orthogonal Matching Pursuit  \\ \hline
\end{tabular}
\caption{Description of algorithms  in \cite{sparselab} used for comparisons}
\end{table}
\end{center}

\mathbf{e}gin{figure}[H]
\center
\includegraphics[scale=0.4]{NMSE_O}
\caption{Root of Normalised Mean Square Error averaged over 200 independent experiments for the signal-to-noise ratios 0 dB, 5 dB, 10 dB, 15 dB, 20 dB, and 25 dB in Example \ref{example:o}.
}
\label{fig:nmseO}
\end{figure}

\mathbf{e}gin{figure}[H]
\center
\includegraphics[scale=0.4]{NMSE_K}
\caption{Root of Normalised Mean Square Error averaged over 200 independent experiments for the signal-to-noise ratios 0 dB, 5 dB, 10 dB, 15 dB, 20 dB, and 25 dB in Example \ref{example:k}.
	}
\label{fig:nmseK}
\end{figure}
 
\mathbf{e}gin{figure}[H]
\center
\includegraphics[scale=0.4]{TIME_O} 
\caption{Computational running time averaged over 200 independent experiments for the signal-to-noise ratios 0 dB, 5 dB, 10 dB, 15 dB, 20 dB, and 25 dB in Example \ref{example:o}.}
\label{fig:timeO}
\end{figure}

\mathbf{e}gin{figure}[H]
\center
\includegraphics[scale=0.4]{TIME_K} 
\caption{Computational running time averaged over 200 independent experiments for the signal-to-noise ratios 0 dB, 5 dB, 10 dB, 15 dB, 20 dB, and 25 dB in Example \ref{example:k}.}
\label{fig:timeK}
\end{figure}

\subsection{Discussion}
It can be seen from Fig. \ref{fig:nmseO} and Fig. \ref{fig:nmseK} that our algorithm outperforms all the other algorithms in \cite{sparselab} in terms of RNMSE.
However, our algorithm requires more computational running time compared to the other algorithms. There are potentially two reasons for this. The first one is that our algorithm is implemented using the CVX package as a parser \cite{grant2008cvx}. Parsers similar to CVX include YALMIP \cite{lofberg2004yalmip}. CVX and YALMIP call generic SDP solvers, e.g., SDPT3 \cite{toh1999sdpt3} or SeDuMi \cite{sturm1999using}, to solve the convex optimisation problem at hand (we use SeDuMi). While these solvers are reliable for wide classes of optimisation problems, they are not specifically optimised in terms of algorithmic complexity to exploit the specific structure of particular problems, such as ours. The second reason comes from the $5^{th}$ step of Algorithm \ref{alg:weight} where the  matrix $Noise\mathbf{I}+\mathbf{\Phi} \mathbf{U}^{(k)} \mathbf{W}^{(k)} \mathbf{\Phi}^\top \in \mathbb{R}^{M \times N}$ has to be inverted to update the weights for the next iteration. Though a pruning rule has been discussed in Remark \ref{remark:threshold}, such inversion at each iteration is inevitable compared to the algorithms considered in \cite{sparselab}. 
It should also be noted from Fig. \ref{fig:nmseO} and Fig. \ref{fig:nmseK} that the RNMSE of some algorithms is independent of the SNR. This may be due to the fact that the coherence of the dictionary matrix is close to $1$ regardless of the SNR. In such case, these algorithms cannot identify the model correctly even when the SNR is high. Therefore, these algorithms are not appropriate for nonlinear system identification.
In future work, we plan to improve our proposed algorithm by exploiting further the structure of the optimisation problem at hand and reducing the associated algorithmic complexity.

\mathbf{e}gin{thebibliography}{1}

\bibitem{leontaritis1985input}
I.~Leontaritis and S.~Billings, ``Input-output parametric models for non-linear
systems part i: deterministic nonlinear systems,'' \emph{International
journal of control}, vol.~41, no.~2, pp. 303--328, 1985.

\bibitem{grant2008cvx}
M.~Grant, S.~Boyd, and Y.~Ye, ``{CVX}: {MATLAB} software for disciplined convex
programming,'' \emph{Online accessible: http://cvxr.com/}, 2008.

\bibitem{elowitz2000son}
M.~Elowitz and S.~Leibler, ``{A synthetic oscillatory network of
transcriptional regulators},'' \emph{Nature}, vol. 403, no. 6767, pp. 335--8,
2000.

\bibitem{strogatz2000kuramoto}
S.~Strogatz, ``From kuramoto to crawford: exploring the onset of
synchronisation in populations of coupled oscillators,'' \emph{Physica D:
Nonlinear Phenomena}, vol. 143, no.~1, pp. 1--20, 2000.

\bibitem{sparselab}
D.~L. Donoho, V.~C. Stodden, and Y.~Tsaig, ``About {SparseLab},'' \emph{Online
accessible: http://sparselab.stanford.edu/}, 2007.

\bibitem{lofberg2004yalmip}
J.~Lofberg, ``{YALMIP}: A toolbox for modeling and optimization in MATLAB,'' in
\emph{Computer Aided Control Systems Design, 2004 IEEE International
Symposium on}, pp.
284--289, 2004.

\bibitem{toh1999sdpt3}
K.-C. Toh, M.~J. Todd, and R.~H. T{\"u}t{\"u}nc{\"u}, ``{SDPT3}--a {MATLAB}
software package for semidefinite programming, version 1.3,''
\emph{Optimization Methods and Software}, vol.~11, no. 1-4, pp. 545--581,
1999.

\bibitem{sturm1999using}
J.~F. Sturm, ``Using {SeDuMi} 1.02, a {MATLAB} toolbox for optimization over
symmetric cones,'' \emph{Optimization methods and software}, vol.~11, no.
1-4, pp. 625--653, 1999.

\end{thebibliography}

\end{document}