\section{Introduction}

Graphs offer a way to describe and explain relationships in complex datasets, a central entity of modern data analysis, where data
deluge is prominent~\cite{shuman2013Emerging,sandryhaila2014big,slavakis2014modeling}. In particular, the nodes of the graph denote the entities and the edges
encode the pairwise relationship between these entities. Such entities are referred to as \emph{graph signals}. Examples of such complex-structured data beyond traditional time-series include data residing on brain networks,  gene networks, social networks, recommendation systems, transportation networks, and so on.

Having a good quality graph is central to any graph signal processing or machine learning task. In this paper, we are interested in the problem of learning the hidden graph topology behind the data. Due to the sheer quantity of data, we are motivated to select the simplest graphical models that adequately explain the data. In particular, we are interested in learning a sparse graph, i.e., a graph with a limited number of edges that adequately explains the input (or training) data. To realize this, we make a simple, but widely used assumption~\cite{shuman2013Emerging,dong2014learning} that the data is smooth with respect to the discovered graph.

The contributions in this paper are threefold. First, we model the graph learning problem as an edge selection problem, where we parameterize the graph through a \emph{sparse edge sampling vector}. In particular, the proposed model provides an elegant handle to control the number of edges, thus the graph sparsity. Second, for the case when the true smooth graph signals are given, the graph learning problem can be solved optimally, and the solution is based on simple \emph{rank ordering}. Finally, given the noisy graph signals, i.e., for the joint sparse graph learning and denoising problem, we provide a one-step solution based on convex optimization as well as an algorithm based on alternating minimization.

The problem of learning the graph Laplacian or the weighted adjacency matrix from smooth graph signals has been considered before~\cite{dong2014learning,kalofolias2016learn}. Learning sparse graphs from the true graph signals, which is the problem we consider in Section~\ref{sec:Noiseless}, has been studied in~\cite{kalofolias2016learn}. There the graph learning problem is posed as a constrained optimization problem with the constraint set being the set of valid adjacency matrices, and the optimization problem is solved using an iterative primal-dual algorithm. In contrast, our modelling greatly simplifies the solution to simple \emph{rank ordering}. Such a modelling is inspired from~\cite{ghosh2006growing}, where the problem to design edge weights that maximize the algebraic connectivity of the graph has been addressed. In \cite{dong2014learning}, the joint graph learning and denoising problem has been addressed, i.e., the problem that we study in Section~\ref{sec:Noisy}. An alternating minimization algorithm is proposed, alternating between graph learning and denoising, where the graph learning optimization problem involves a search over the space of all valid graph Laplacians. On the contrary,
we show that this problem can be solved in one-step and it boils down to the design of a sparse edge sampling function.

Graph topology identification is also investigated in \cite{segarra2016network} under the assumption that the eigenvectors of the graph Laplacian are known, which is a much stronger assumption. Although the eigenvectors can be computed from graph data (or the sample covariance matrix) when it is stationary with respect to the graph~\cite{marques2016stationary,chepuri2016subsampling}, the graph signals need not always be vertex stationary. In any case, the estimated eigenvectors are not error free due to limited data records. In most of the existing approaches~\cite{dong2014learning,kalofolias2016learn,segarra2016network}, graph sparsification is (or can be) achieved by penalizing the $\ell_1$-norm of the graph Laplacian matrix, adjacency matrix or the shift operator, however, there is no explicit handle to control the number of edges, unlike the proposed approach. In a related line of research, \cite{spielman2011spectral,batson2012twice} investigate computing sparse graphs that approximate a given graph spectrally, which means that their Laplacian matrices have similar quadratic forms.

\section{Problem setup}

Consider a dataset with $N$ real valued elements, which are defined on the vertices of an undirected graph $\mathcal{G} = (\mathcal{V},\mathcal{E})$, where the vertex set $\mathcal{V} = \{v_1,\cdots,v_N\}$ denotes the set of nodes, and the edge set $\mathcal{E}$ reveals the connection between the nodes. We refer to such datasets as {\it graph signals}.  We assume that the length of the graph signals (thus the number of nodes), i.e., $N$ is known. However, the edge set is not known. Therefore, we assume a \emph{complete graph} $\mathcal{G}(\mathcal{V},\mathcal{E})$ as a candidate graph in which each node is connected to every other node with the number of edges $|\mathcal{E}| = M = N(N-1)/2$, and aim to determine a subgraph of $\mathcal{G}$ by choosing a subset of edges, $\mathcal{E}_s$, from the edge set $\mathcal{E}$ of this candidate graph.

Any undirected graph topology is basically determined by its graph Laplacian matrix, which essentially reveals the connectivity of the graph. Let us denote the graph Laplacian matrix (i.e., a symmetric matrix) of the complete graph by ${\boldsymbol L} \in \mathbb{S}^{N}$, where $[{\boldsymbol L}]_{i,j}$ is nonzero, say equal to 1, only if $i=j$ or $(i,j) \in \mathcal{E}$. The symmetric matrix ${\boldsymbol L}$ can be expressed in terms of the so-called \emph{incidence matrix}, ${\boldsymbol A}  = [{\boldsymbol a}_1,\cdots,{\boldsymbol a}_M] \in \mathbb{R}^{N \times M}$ as
\[
{\boldsymbol L} = {\boldsymbol A}{\boldsymbol A}^T = \sum_{m=1}^M  {\boldsymbol a}_m{\boldsymbol a}_m^T,
\]
where the $m$-th column of ${\boldsymbol A}$, i.e., ${\boldsymbol a}_m$ denotes a length-$N$ edge vector with entries $[{\boldsymbol a}_m]_i = 1$, $[{\boldsymbol a}_m]_j= -1$, and it has zeros elsewhere, for an edge $m$ connecting nodes $i$ with $j$ (more generally, ${\boldsymbol a}_m$ is determined only up to a sign).

Let us now denote the subgraph $\mathcal{G}_s(\mathcal{V},\mathcal{E}_s)$ with the edge set $\mathcal{E}_s \subset \mathcal{E}$ such that $|\mathcal{E}_s| = K \ll M$. We will refer to such a subgraph with $K$ edges as a $K$-sparse graph. We connect such a $K$-sparse graph
$\mathcal{G}_s$  to ${\boldsymbol L}$ through a \emph{sparse edge selection} vector ${\boldsymbol w} = [w_1,w_2,\cdots,w_M]^T \in \{0,1\}^M$, where $w_m = 1$ if an edge belongs to the edge subset $\mathcal{E}_s$, and $w_m = 0$ otherwise. In terms of  ${\boldsymbol w}$, $|\mathcal{E}_s| = K$ means $\|{\boldsymbol w}\|_0 = K$. (The notation $\|{\boldsymbol w}\|_0$ counts the number of non-zero entries in ${\boldsymbol w}$.)
Finally, we can write the Laplacian matrix of the $K$-sparse graph, ${\boldsymbol L}_s$, as a function of ${\boldsymbol w}$ as
\begin{equation}
\label{eq:graph_w}
{\boldsymbol L}_s({\boldsymbol w}) = \sum_{m=1}^Mw_m {\boldsymbol a}_m{\boldsymbol a}_m^T.
\end{equation}

In what follows, we will optimally design the edge sampling function ${\boldsymbol w}$ to recover the graph that sufficiently explains the data.

\section{Learning from Noiseless Graph Signals} \label{sec:Noiseless}

Let ${\boldsymbol x} = [x_1,x_2,\cdots,x_N]^T \in \mathbb{R}^N$ be a graph signal defined on the vertices $\mathcal{V}$ of a graph.
The smoothness and the spectral content of the signal both depend on the underlying graph topology. The Laplacian quadratic form given by ${\boldsymbol x}^T{\boldsymbol L}_s({\boldsymbol w}){\boldsymbol x}$ quantifies how smooth the graph signal ${\boldsymbol x}$ is with respect to the underlying graph~\cite{shuman2013Emerging}. In particular, the signal ${\boldsymbol x}$ is smoothest with respect to the graph with $K$ edges for low values of ${\boldsymbol x}^T{\boldsymbol L}_s({\boldsymbol w}){\boldsymbol x}$.

\subsection{Problem statement: noiseless setting}
Suppose we are given $L$ graph signals denoted by the vectors $\{{\boldsymbol x}_{k}\}_{k=1}^L$, and they are collected in an $N \times L$ matrix
${\boldsymbol X} = [{\boldsymbol x}_1,\cdots,{\boldsymbol x}_L]$. We are interested in recovering the graph Laplacian (in other words, the graph topology) under the prior
information that the graph signals are smooth with respect to a $K$-sparse graph. More formally, we state the following.

\begin{myprob}Given the graph signals $\{{\boldsymbol x}_k\}_{k=1}^{L}$, determine a graph with $K$ edges such that the graph signals have smooth variations on the resulting graph.
\end{myprob}

Mathematically, the above problem can be cast as the following optimization problem:
\begin{equation}
\label{eq:graph_learning_1}
\argmin_{{\boldsymbol w} \in \mathcal{W}}  \quad \frac{1}{L}\sum_{k=1}^L {\boldsymbol x}_k^T{\boldsymbol L}_s({\boldsymbol w}){\boldsymbol x}_k  = \frac{1}{L} {\rm tr}\{{\boldsymbol X}^T{\boldsymbol L}_s({\boldsymbol w}){\boldsymbol X}\},
\end{equation}
where
$
\mathcal{W} = \{{\boldsymbol w} \in \{0,1\}^M \, \vert \, \|{\boldsymbol w}\|_0 = K\}
$ is the constraint set that restricts the number of edges.

\subsection{Solver}
Problem \eqref{eq:graph_learning_1} is a cardinality constrained Boolean optimization problem, hence nonconvex.
 By recalling that ${\boldsymbol L}_s({\boldsymbol w}) = \sum_{m=1}^Mw_m {\boldsymbol a}_m{\boldsymbol a}_m^T$, we can express the cost function in  \eqref{eq:graph_learning_1} as a linear function in ${\boldsymbol w}$, i.e., we have
\begin{equation}
\frac{1}{L} {\rm tr}\left\{{\boldsymbol X}^T{\boldsymbol L}_s({\boldsymbol w}){\boldsymbol X}\right\} =  \sum_{m=1}^M w_m {\rm tr}\left\{{\boldsymbol X}^T ({\boldsymbol a_m}{\boldsymbol a_m}^T){\boldsymbol X}\right\}.
\end{equation}
Introducing the length-$M$ vector ${\boldsymbol c} = [c_1,c_2,\ldots,c_M]^T$ with $c_m = {\rm tr}\left\{{\boldsymbol X}^T ({\boldsymbol a_m}{\boldsymbol a_m}^T){\boldsymbol X}\right\}$, we can write \eqref{eq:graph_learning_1} as
\begin{equation}
\label{eq:graph_sorting_problem}
\argmin_{{\boldsymbol w} \in \{0,1\}^M}  \quad {\boldsymbol c}^T{\boldsymbol w} \quad {\rm s. to}  \quad \|{\boldsymbol w}\|_0 = K.
\end{equation}
The above Boolean linear programming problem admits an explicit solution and computing the optimal solution is straightforward. It is solved by sorting the entries of ${\boldsymbol c}$ in an ascending ordering. More specifically, the solution ${\boldsymbol w}$ will have entries equal to 1 at indices corresponding to the $K$ smallest entries of ${\boldsymbol c}$, and others are set to zero (ties may be broken arbitrarily).
Computationally, the sorting algorithm costs $\mathcal{O}(K \log K)$, and with a parallel implementation (e.g., on different processors), the computational complexity will be as low as $\mathcal{O}(K)$~\cite{blum2008energy,chepuri2015sparse}.
We give another interpretation of this result through the following remark.

\begin{myrem}
Let us suppose the graph signal is stochastic with covariance matrix ${\boldsymbol R}_{\boldsymbol x} = \mathbb{E}\{{\boldsymbol x}{\boldsymbol x}^T\} \in \mathbb{R}^{N \times N}$. Then, the solution to~\eqref{eq:graph_sorting_problem} would select $K$ edges between those nodes having the highest cross-correlation, i.e., it will
add an edge between the $i$th and the $j$th node if the variables $x_i$ and $x_j$ are strongly correlated.  To see this, we express the cost function in \eqref{eq:graph_learning_1} as
\begin{equation*}
\begin{aligned}
 L^{-1}{\rm tr}\{{\boldsymbol X}^T{\boldsymbol L}_s({\boldsymbol w}){\boldsymbol X}\} &=  {\rm tr}\left\{{\boldsymbol L}_s({\boldsymbol w})\widehat{\boldsymbol R}_{\boldsymbol x}\right\} \nonumber\\ &= \sum_{m=1}^M w_m  ({\boldsymbol a_m}^T\widehat{\boldsymbol R}_{\boldsymbol x}{\boldsymbol a}_m)
 \end{aligned}
\label{eq:inverse_cov}
\end{equation*}
where $\widehat{\boldsymbol R}_{\boldsymbol x} =  \frac{1}{L}{\boldsymbol X}{\boldsymbol X}^T \in \mathbb{R}^{N \times N}$ is the sample data covariance matrix.
Recalling the definition of ${\boldsymbol a_m}$, it is easy to see that the term ${\boldsymbol a_m}^T\widehat{\boldsymbol R}_{\boldsymbol x}{\boldsymbol a}_m = [\widehat{\boldsymbol R}_{\boldsymbol x}]_{i,i} + [\widehat{\boldsymbol R}_{\boldsymbol x}]_{j,j} - 2[\widehat{\boldsymbol R}_{\boldsymbol x}]_{i,j}$ is small if the $i$th and $j$th nodes are highly correlated and we have sufficient samples to compute the sample covariance matrix.
\end{myrem}
By modelling the graph topology through an {edge selection vector}, the graph learning problem can be solved optimally using a simple and elegant solution with a controlled sparsity level, whereas optimizing directly the graph Laplacian~\cite{dong2014learning} or the adjacency matrix~\cite{kalofolias2016learn} leads to a more complicated suboptimal solution with no explicit handle to control the graph sparsity.

\section{Learning from Noisy Graph Signals} \label{sec:Noisy}
In many cases, we might not have access to the true graph signals. Suppose we observe a noisy version of the graph signal,  ${\boldsymbol x}_k$, as
\begin{equation}
\label{eq:observationmodel}
{\boldsymbol y}_k = {\boldsymbol x}_k + {\boldsymbol n}_k \in \mathbb{R}^N,
\end{equation}
and we are given $L$ such observations for $k=1,2,\ldots,L$, where we assume that ${\boldsymbol n}_k$ is zero-mean white Gaussian noise of variance $\sigma^2$.  To recover ${\boldsymbol x}_k$ based on the smoothness assumption, typically a least-squares problem is solved with a \emph{Tikhonov regularization}, ${\boldsymbol x}_k^T{\boldsymbol L}{\boldsymbol x}_k$, to enforce the prior information that the noiseless graph signal ${\boldsymbol x}_k$ is smooth with respect to the underlying graph. More specifically, the following optimization problem (assuming, for a moment that the graph, i.e., ${\boldsymbol w}$ is known) is solved~\cite{shuman2013Emerging}:
\begin{equation}
\label{eq:graphdenois}
\argmin_{\{{\boldsymbol x}_k\}_{k=1}^L} \quad  \frac{1}{L}\sum_{k=1}^L \left( \|{\boldsymbol y}_k - {\boldsymbol x}_k\|_2^2 + \gamma {\boldsymbol x}_k^T{\boldsymbol L}_s({\boldsymbol w}){\boldsymbol x}_k \right),
\end{equation}
where the regularization parameter $\gamma > 0$ controls the amount of smoothness. This graph denoising problem has an explicit solution given by
\[
\widehat{\boldsymbol x}_k = [{\boldsymbol I} + \gamma {\boldsymbol L}_s({\boldsymbol w})]^{-1}{\boldsymbol y}_k, \, k=1,\cdots,L.
\]
\vskip-1mm
\subsection{Problem statement: noisy setting}
\vskip-1mm
Having given the above denoising inference problem at hand, we will now formally state the problem of interest.

\begin{myprob}Given the observations $\{{\boldsymbol y}_k\}_{k=1}^{L}$ that is related to the unknown graph signal ${\boldsymbol x}_k$ as in \eqref{eq:observationmodel}, determine the $K$-sparse graph such that the estimate $\widehat{\boldsymbol x}_k$ has the lowest possible estimation error, and it is smooth with respect to the recovered graph.
\end{myprob}

Sparse graph learning for the denoising inference problem can be mathematically formulated as follows:
\begin{equation}
\label{eq:learning_prob}
\argmin_{\{{\boldsymbol x}_k\}_{k=1}^L, {\boldsymbol w} \in \mathcal{W}} \, \frac{1}{L}\sum_{k=1}^L ( \|{\boldsymbol y}_k - {\boldsymbol x}_k\|_2^2 + \gamma \, {\boldsymbol x}_k^T{\boldsymbol L}_s({\boldsymbol w}){\boldsymbol x}_k)
\end{equation}
whose solution is denoted as $(\{\widehat{\boldsymbol x}_k\}_{k=1}^L, \widehat{\boldsymbol w})$.
This formulation is different from~\cite{dong2014learning}, as \cite{dong2014learning}  solves an optimization problem over the
space of all possible graph Laplacians (instead of parameterizing the graph with ${\boldsymbol w} \in \mathcal{W}$) without sparsifying the graph. It can nevertheless be done through an extra $\ell_1$-norm penalty term.

The above problem \eqref{eq:learning_prob} is noncovex due to the Boolean and cardinality constraints on ${\boldsymbol w}$ and the coupling between the optimization variables in the second term of \eqref{eq:learning_prob}. We provide two methods to solve it. The first one is a straightforward approach based on alternating descent, while the second one is based on convex relaxation.
\subsection{Alternating minimization} \label{sec:am}
The optimization problem \eqref{eq:learning_prob} can be solved using alternating minimization with respect to $\{{\boldsymbol x}_k\}_{k=1}^L$ and
${\boldsymbol w}$. That is,  given ${\boldsymbol w}$, the problem in \eqref{eq:learning_prob} reduces to a linear system in the unknown ${\boldsymbol X}$, which admits a closed form solution; while given $\{{\boldsymbol x}_k\}_{k=1}^L$, it reduces to a Boolean linear programming problem, which admits an analytical solution with respect to ${\boldsymbol w}$ based on rank ordering. These observations suggest an iterative alternating minimization algorithm yielding successive estimates of $\{{\boldsymbol x}_k\}_{k=1}^L$ with fixed ${\boldsymbol w}$, and alternately of ${\boldsymbol w}$ with fixed
$\{{\boldsymbol x}_k\}_{k=1}^L$. Specifically, with the iterate of ${\boldsymbol w}$ given per iteration $i \geq 0$,  i.e., ${\boldsymbol w}[i]$, we solve for ${\boldsymbol X}[i]$ using a matrix inversion as
\[
{\boldsymbol X}[i] = {\boldsymbol X}_{\rm min}({\boldsymbol w}[i])
\]
with
\begin{equation}
\label{eq:alt_w}
\begin{aligned}
 {\boldsymbol X}_{\rm min}({\boldsymbol w}) &= \argmin_{{\boldsymbol X}} \,  \|{\boldsymbol Y} - {\boldsymbol X}\|_F^2 + \gamma \, {\rm tr} \{{\boldsymbol X}^T{\boldsymbol L}_s({\boldsymbol w}){\boldsymbol X}\} \\
 &=[{\boldsymbol I} + \gamma {\boldsymbol L}_s({\boldsymbol w})]^{-1}{\boldsymbol Y},
\end{aligned}
\end{equation}
where ${\boldsymbol Y} = [{\boldsymbol y}_1,{\boldsymbol y}_2,\cdots, {\boldsymbol y}_L]$ is the data matrix of size $N \times L$.

Once ${\boldsymbol X}[i]$ is available, ${\boldsymbol w}[i+1]$ can be obtained by solving the Boolean linear program [cf. \eqref{eq:graph_sorting_problem}]
\[
{\boldsymbol w}[i+1] = \argmin_{{\boldsymbol w} \in \{0,1\}^M}  \quad \sum_{m=1}^M w_mc_m[i+1] \quad {\rm s. to}  \quad \|{\boldsymbol w}\|_0 = K,
\]
where $c_m[i+1] = {\rm tr}\left\{{\boldsymbol X}^T[i+1] ({\boldsymbol a_m}{\boldsymbol a_m}^T){\boldsymbol X}[i+1]\right\}$. In spite of the Boolean and cardinality constraints in the above problem, there exists a simple analytical solution for ${\boldsymbol w}[i+1]$ based on sorting $\{c_m[i+1]\}_{m=1}^M$, i.e., the solution ${\boldsymbol w}[i+1]$ will have entries equal to 1 at indices corresponding to the $K$ smallest entries in $\{c_m[i+1]\}_{m=1}^M$ and zeros otherwise. The iterations are initialized at $i = 0$ by randomly generating ${\boldsymbol w}[i+1]$ from a uniform distribution over $\mathcal{W}$. The above alternating minimization method is computationally very attractive, and consists of two simple known solutions per iteration. However, the algorithm converges only to a stationary point of \eqref{eq:learning_prob}, and it suffers from the choice of the initial estimate.

The algorithm proposed in \cite{dong2014learning} is also along the lines of alternating minimization, except that the graph learning step involves a complicated optimization over the space of all possible valid Laplacian matrices.
\vskip-3mm
\subsection{Convex relaxation} \label{sec:cvx}
\vskip-1mm
To avoid the issues related to the initialization of the alternating minimization algorithm, in what follows we propose a one-step solution based on convex relaxation. We can rewrite the formulation in \eqref{eq:learning_prob} alternatively as
\begin{equation}
\widehat{\boldsymbol w} = \argmin_{{\boldsymbol w} \in \mathcal{W}} \quad  r({\boldsymbol w});  \quad
\widehat{\boldsymbol X} = {\boldsymbol X}_{\rm min}(\widehat{\boldsymbol w}) \label{eq:learning_prob2}
\end{equation}
with
\[
r({\boldsymbol w}) = \|{\boldsymbol Y} - {\boldsymbol X}_{\rm min}({\boldsymbol w})\|_F^2
+ \gamma \, {\rm tr}\{{\boldsymbol X}_{\rm min}^T({\boldsymbol w}){\boldsymbol L}_s({\boldsymbol w}){\boldsymbol X}_{\rm min}({\boldsymbol w})\}
\]
and
\begin{equation}
\label{eq:learning_xmin}
 [{\boldsymbol I} + \gamma {\boldsymbol L}_s({\boldsymbol w})] {\boldsymbol X}_{\rm min}({\boldsymbol w}) ={\boldsymbol Y}.
\end{equation}
The computational complexity of solving the linear system of equations \eqref{eq:learning_xmin} decreases as the sparsity in ${\boldsymbol w}$ increases. Furthermore, the estimates $\{\widehat{\boldsymbol x}_k\}_{k=1}^L$ and $\widehat{\boldsymbol w}$ in \eqref{eq:learning_prob2} are still the same as in \eqref{eq:learning_prob}.
\begin{figure*}[!t]
\centering
       		\psfrag{celcius}{\footnotesize $^\circ$ C}
		\psfrag{5}{\tiny 5} \psfrag{6}{\tiny 6} \psfrag{7}{\tiny 7} \psfrag{8}{\tiny 8} \psfrag{9}{\tiny 9} \psfrag{10}{\tiny 10} \psfrag{11}{\tiny 11}
\begin{subfigure}[t]{0.4\textwidth}
                \includegraphics[width=\columnwidth, height=1.5in]{noiseless_graph.eps}
                \caption{Noiseless setting}
                \label{fig:sparsegraph_noisless}
\end{subfigure}%
~
\begin{subfigure}[t]{0.4\textwidth}
                \includegraphics[width=\columnwidth,height=1.5in]{noisy_graph.eps}
                \caption{Noisy setting}
                \label{fig:sparsegraph_noisy_am}
\end{subfigure
        \caption{\footnotesize{\emph{Sparse graph learning}: The colored dots indicate the temperature values. (a) Noiseless case. Graph with $K=110$ edges recovered by solving \eqref{eq:graph_learning_1}. (b) Noisy case: Convex relaxation is used to recover a graph with $K=110$ edges using \eqref{eq:convexrelax}.}}
        \label{fig:sparsegraph}
\vskip-4mm
\end{figure*}

Plugging the solution to \eqref{eq:learning_xmin} in $r({\boldsymbol w})$ and after some straightforward matrix algebra, we can express the regularized residual squared, $r({\boldsymbol w})$,
as
\begin{equation}
\begin{aligned}
r({\boldsymbol w})
&={\rm tr}\left\{{\boldsymbol Y}^T [{\boldsymbol I} + \gamma {\boldsymbol L}_s({\boldsymbol w})]^{-1}{\boldsymbol Y}\right\} \\
&\hskip12mm+ \gamma {\rm tr}\left\{{\boldsymbol Y}^T {\boldsymbol L}_s({\boldsymbol w}){\boldsymbol Y}\right\} - \|{\boldsymbol Y}\|_F^2.
\end{aligned}
\end{equation}

Relaxing the cardinality constraint $\|{\boldsymbol w}\|_0 = K$ with ${\boldsymbol 1}^T{\boldsymbol w} = K$ and the Boolean constraints $\{0,1\}^M$ with linear inequality constraints related to the box constraint $[0,1]^M$, the optimization problem  \eqref{eq:learning_prob2} will be convex on ${\boldsymbol w} \in [0,1]^M$. To see this, we introduce a variable
\[
{\boldsymbol Z} = {\boldsymbol Y}^T [{\boldsymbol I} + \gamma {\boldsymbol L}_s({\boldsymbol w})]^{-1}{\boldsymbol Y} + \gamma {\boldsymbol Y}^T {\boldsymbol L}_s({\boldsymbol w}){\boldsymbol Y} \in \mathbb{R}^{L \times L}
\]
and obtain a semidefinite program:
\begin{equation}
\begin{aligned}
\label{eq:convexrelax}
&\argmin_{{\boldsymbol Z},{\boldsymbol w}} \quad  {\rm tr}\{{\boldsymbol Z}\}\\
& {\rm s.to} \quad\left[\begin{array}{cc}{\boldsymbol Z} - \gamma {\boldsymbol Y}^T{\boldsymbol L}_s({\boldsymbol w}){\boldsymbol Y} & {\boldsymbol Y}^T \\{\boldsymbol Y} & {\boldsymbol I} + \gamma {\boldsymbol L}_s({\boldsymbol w}) \end{array}\right] \succeq {\boldsymbol 0}_{L+N},
\\
&\hskip10mm{\boldsymbol 1}^T{\boldsymbol w}= K,  \,\, 0 \leq w_m \leq 1, m=1,2,\ldots,M,
\end{aligned}
\end{equation}
with variables ${\boldsymbol w}$ and ${\boldsymbol Z}$, and recall that ${\boldsymbol L}_s({\boldsymbol w}) = \sum_{m=1}^Mw_m {\boldsymbol a}_m{\boldsymbol a}_m^T$. A standard off-the-shelf solver can be used for solving the semidefinite program in \eqref{eq:convexrelax}. For large-scale
problems, computationally cheaper first-order (and online) methods for solving \eqref{eq:convexrelax} can be derived as the size of the linear matrix inequality in \eqref{eq:convexrelax} depends on the size of the training data and the number of nodes.
\begin{figure}[tb]
\centering
       		\psfrag{smoothness}{\hskip-3mm\footnotesize${\rm tr}\{{\boldsymbol X}^T{\boldsymbol L}_s({\boldsymbol w}){\boldsymbol X}$\}}
		\psfrag{Noisless prop12345678}{\tiny proposed (noiseless, optimal)}
		\psfrag{Noisless exist12345678}{\tiny primal-dual~\cite{kalofolias2016learn}}
		\psfrag{no of edges}{\footnotesize No. of edges}
\begin{subfigure}[t]{0.8\columnwidth}
                \includegraphics[width=\columnwidth,height=1.5in]{smoothness.eps}
                \caption{Noiseless setting: smoothness for different graphs}
                \label{fig:smooth_pereval}
\end{subfigure}%
\\[1em]
\psfrag{noise level}{\footnotesize Noise level, $\sigma$}
\psfrag{error}{\hskip-8mm\footnotesize mean squared error}
\psfrag{Noisless prop12345678}{\tiny proposed (noiseless, sorting)}
\psfrag{Noisy am12345678}{\tiny alternating min. in Sec.~\ref{sec:am}}
\psfrag{Noisy exist12345678}{\tiny alternating min. from~\cite{dong2014learning}}
\psfrag{Noisy cvx12345678}{\tiny one-step cvx opt in Sec.~\ref{sec:cvx}}
\begin{subfigure}[t]{0.8\columnwidth}
                \includegraphics[width=\columnwidth, height=1.5in]{mse.eps}
                \caption{Denoising error for different noise levels}
                \label{fig:mse_pereval}
\end{subfigure
        \caption{\footnotesize{Performance evaluation.}}
        \label{fig:per_eval}
\vskip-4mm
\end{figure}
\vskip-2mm
\section{Numerical results}
\vskip-1mm
We use temperature measurements collected across $32$ weather stations in the French region of Brittany and the aim is to learn the graph that explains the observed data; see Fig.~\ref{fig:sparsegraph}. There are $744$ observations per weather station available, out of which we use $L=50$ snapshots as the training set and the remaining ones as the evaluation set. One such observation (i.e., a graph signal) on a graph with $N=32$ nodes is shown in Fig.~\ref{fig:sparsegraph}, where the colored dots indicate different temperature readings. The convex optimization problems are solved using the {\tt CVX} toolbox, which internally calls {\tt SDPT3}~\cite{cvx}\let\thefootnote\relax\footnote{\noindent Software and datasets to produce results of this paper can be downloaded from  {\url{http://cas.et.tudelft.nl/~sundeep/sw/icassp17Graphlearning.zip}}}. The candidate graph with $N=32$ will have $M=496$ edges, from which we aim to learn a subgraph with $K=110$ edges.

To begin with, we consider the noiseless case, where the true graph signal is assumed to be known, and graph learning in this case amounts to solving a sorting problem. As shown in Fig.~\ref{fig:sparsegraph_noisless}, we can see that in the learnt graph with $K=110$ edges, edges are present  between nodes that share similar values. Although the proposed approach doesn't always (e.g., for low values of $K$) ensure a well-connected graph, it clusters entities (or correlated nodes) with similar values. Fig.~\ref{fig:smooth_pereval} shows that the cost (i.e., smoothness) of the proposed closed-form sorting solution, which is optimal, is lower than the existing iterative solution~\cite{kalofolias2016learn}.

Next, we consider the noisy setting with the same training data as before, where we perform joint graph learning and denoising. In Fig.~\ref{fig:sparsegraph_noisy_am}, we show the learnt graph with $K=110$ edges based on the convex relaxation approach explained in Sec.~\ref{sec:cvx}. In Fig.~\ref{fig:mse_pereval}, we evaluate the denoising performance based on the learnt graph using the evaluation set. In particular, we show the mean squared error for different values of the noise level, where the mean squared error is computed from $1000$ independent Monte Carlo experiments.
The one-step solution based on convex optimization (cf. Sec.~\ref{sec:cvx}) leads to a lower error as compared to the alternating minimization approaches, which in general converge only to a stationary point. This also holds for our method developed in Sec.~\ref{sec:am}, however, we stress the fact that the proposed alternating minimization (cf. Sec.~\ref{sec:am}) is computationally much less expensive (involving two simple known solutions per iteration) as compared to the iterative solution in~\cite{dong2014learning}. The graph learnt under the noiseless setting does not perform well for denoising. Nevertheless, due its simple solution, it can be used to generate a base graph, which can be further refined for specific graph inference problems.

\vskip-4mm
\section{Conclusions}
\vskip-2mm
We have studied the problem of learning a sparse graph that adequately explains the data under a smoothness prior. We model the graph learning problem as the design of a sparse edge sampling function. In other words, we express the graph Laplacian in terms of an edge selection vector. We have considered both the noiseless and noisy setting. In the noiseless setting, designing the edge selection vector is elegant, and it boils down to a simple low-complexity sorting problem. However, in the presence of noise, we propose a computationally cheap alternating minimization algorithm as well as a one-step convex relaxation based solution.

\pagebreak
\bibliographystyle{IEEEtran}


\section{Background}

\subsection{Speeding up Python computation}

Several methods exist to overcome the inherent interpreter overhead in
Python, most of which focus on using pre-compiled functions to perform
computationally intense loops. The most prominent example of this is
the NumPy{}~\cite{vanDerWalt2011} array interface, which manages
array data in contiguous blocks to avoid memory copies and provides
vectorized calculations, where operations are applied element-wise
over arrays using pre-compiled loops. The same technique is used by
the SciPy package that adds a wealth of pre-compiled mathematical
operators that can be applied in this vectorized way.

The effectiveness of NumPy{} vectorization, however, is limited by
the need to create arrays with matching data layout, which often
results in the ``outer loop'' being a pure Python loop. A more
aggressive pre-compilation can be achieved with
Cython~\cite{Behnel2011}, which allows users to annotate Python code
with additional type information and invoke an explicit compilation
step at install time at the expense of flexibility. Even more dynamic
compilation can be achieved with Numba~\cite{Lam2015}, a Just-In-Time
(JIT) compiler that focuses on array-oriented computation and utilises
LLVM to compile a subset of the Python language at runtime.

The Symbolic Python package SymPy{} used in Devito also provides
JIT compilation via the \lstinline!lambdify! utility function that
compiles a symbolic expression into either a Python function that
supports NumPy{} vectorization or pure C or Fortran kernels.

\subsection{DSLs embedded in Python}

Interest in building generic DSLs for solving PDEs is not new with
early attempts dating back as far as
1970~\cite{Cardenas1970,Grant1988,Yukio1985,VanEngelen1996}. More
recently, two prominent finite element software packages,
FEniCS~\cite{Logg2012} and Firedrake~\cite{Rathgeber2015}, have
demonstrated the power of symbolic computation using the DSL
paradigm. Both packages implement the Universal Form Language
(UFL)~\cite{Alnaes2014} that allows scientists to express complex
finite element problems symbolically in the weak form. Reducing the
complexity of a problem to its symbolic definition then enables
further symbolic abstraction, such as the automated generation of
adjoint models, as demonstrated by Dolfin-Adjoint~\cite{Farrell2013}.
It is noteworthy that FEniCS, Firedrake and Dolfin-Adjoint all provide
a native Python interface to the user, with the DSL embedded in the
Python language itself.

The optimization of regular grid and stencil computations has also
produced a vast range of libraries and DSLs that aim to ease the
efficient automated creation of high-performance
codes~\cite{Brandvik2010,Henretty2013,Zhang2012,Datta2008,PATUS,STARGATES,HiPACC}.
Most stencil DSLs, however, define their own custom input language,
which limits their scope and applicability for practical
applications. Since solving realistic problems often requires more
than just a fast and efficient PDE solver, scientific applications can
benefit from DSLs embedded in Python directly due to the native
compatibility with the Python software ecosystem.

\subsection{Abstraction layers in DSL design}
\label{sec:dsl-abstraction-layers}

While DSL-based approaches promise great flexibility and productivity,
care needs to be taken when designing such domain-specific software
stacks. In particular the choice of abstraction layers and interfaces
is crucial, since future development and practical use can often be
hampered by too narrow abstractions, while all-encompassing packages
often become too generic, complicated and hard to maintain. The
separation of the high-level symbolic problem definition layer from
the low-level implementation details of the generated code is of
crucial importance here, since it provides a separation of concerns
between the domain scientists and HPC specialists. This not only
increases productivity on both sides, but also ensures maintainability
of the generated software stack.

Active libraries, such as OP2~\cite{Mudalidge2012} for unstructured
mesh applications or OPS~\cite{Reguly2014} for structured grids, aim
to decouple the problem specification from the underlying execution by
providing an abstraction of data movement following the access-execute
paradigm~\cite{Howes2008}. A particular implementation of such an
abstraction is provided by PyOP2~\cite{Rathgeber2012}, which acts as
the underlying unstructured mesh traversal engine for the previously
mentioned Firedrake finite element system. Similar mechanisms
are provided by the Halide project~\cite{Halide}, a domain specific
language, compiler, and auto-tuning system that was originally designed
for image processing pipelines

\section{Acknowledgements}

This work was financially supported in part by the Natural Sciences
and Engineering Research Council of Canada Collaborative Research and
Development Grant DNOISE II (CDRP J 375142-08), EPSRC grant
EP/L000407/1, the Imperial College London Intel Parallel Computing
Centre, SENAI CIMATEC and the MCTI (Federal Govt of Brazil)
scholarship MCTI/FINEP/CNPq 384214/2015-0. This research was carried
out as part of the SINBAD II project with the support of the member
organizations of the SINBAD
Consortium.

\bibliographystyle{IEEEtran}

\section{Devito design and API}
\label{sec:devito}

In contrast to traditional Python-embedded DSLs, Devito does not
define its own high-level language {\em per se}. Instead, it relies
on SymPy{} to provide the building blocks for users to express a
finite difference operator as a symbolic equation. Devito, however,
annotates these objects with additional meta-data and provides new
abstractions specific to defining finite difference problems. This
allows Devito to automate the transformation of high-level symbolic
expressions to optimized executable C code by maintaining and
manipulating SymPy{} expression objects throughout the code generation
process, while leveraging symbolic manipulation utilities provided by
SymPy{}.

The automated code generation process is performed in multiple steps
and follows the {\em Principle of Graceful Degradation}, which
suggests that users should be able to circumvent limits to the top
layer abstraction by utilising a lower-level API. A common example for
this is sparse point interpolation, which is required by many
scientific applications to extract solution values at specific points,
but can not easily be expressed as a high-level PDE. For such cases
Devito provides users with a lower level API to allow custom
expressions to be added to the generated code in a C-like manner.

\figref{fig:devito_layers} highlights the overall design hierarchy of
Devito and shows the individual layers of transformations. At the top,
user can define data object that associate \lstinline!sympy.Function!
symbols with data buffers and use them to derive stencil equations in
a concise symbolic format. The resulting stencil expressions are
automatically expanded and transformed into explicit array accesses
when the user builds a \lstinline!devito.Operator! object from the
SymPy{} expression. A secondary low-level kernel API is then utilised
to generate optimized C code from the intermediate representation, and
a set of compiler presets is used to perform Just-in-Time (JIT)
compilation before executing the kernel operator over the given data.

\input{figures/devito_layers}

\subsection{High-level API: Derivatives and equations}
\label{sec:devito_high}

One of the core ideas behind Devito is to associate simulation data
with symbolic objects to leverage the symbolic manipulation
capabilities provided by SymPy{} during code generation. For this
purpose Devito provides \lstinline!DenseData! and \lstinline!TimeData!
objects that behave like \lstinline!sympy.Function! objects, but also
encapsulate data buffers and present them to the user as wrapped
NumPy{} arrays. Importantly, Devito is thus able to manage data
allocation internally to provide data caching on disk and first touch
allocation to ensure data locality on NUMA systems.

To ensure associated meta-data is not lost during symbolic
manipulation, Devito maintains its own symbol cache that
``re-attaches'' Devito-specific information whenever a new instance of
a symbol is created. Since SymPy{} uses frequent re-instantiation of
\lstinline!Function! objects to facilitate symbolic substitutions in
expressions, this enables Devito symbols to behave opaquely during
symbolic manipulation, for example when forming finite difference
stencils. The following example demonstrates how a two-dimensional
symbolic data object \lstinline!f(x, y)!  maintains the shape
information of its associated data buffer (\lstinline!f.data!), while
expanding the first derivative in dimension \lstinline!x! using
SymPy{}'s \lstinline!as_finite_diff! utility.

\begin{minipage}{\linewidth}
\begin{lstlisting}
In [1]: from Devito import DenseData
In [2]: from sympy import as_finite_diff
In [3]: from sympy.abc import x, h

In [4]: f = DenseData(name='f', shape=(10, 12))
Out[4]: f(x, y)

In [5]: eqn = as_finite_diff(f.diff(x), [x, x+h])
Out[5]: -f(x, y)/h + f(h + x, y)/h

In [6]: f_h = eqn.args[0].args[1]
Out[6]: f(h + x, y)

In [7]: f_h.data.shape
Out[7]: (10, 12)
\end{lstlisting}
\end{minipage}

Moreover, Devito objects also provide shorthand notation for common
finite difference formulations, such as first, second and cross
derivatives in the space and time dimensions. The order of the stencil
discretization is hereby defined on the symbolic object itself, which
allows us to express the second derivative of \lstinline!f! in
\lstinline!x! as:

\begin{minipage}{\linewidth}
\begin{lstlisting}
In [1]: from devito import DenseData

In [2]: f = DenseData(name='f', shape=(10, 10),
                      space_order=2)
Out[2]: f(x, y)

In [3]: f.dx2
Out[3]: -2*f(x, y)/h**2 + f(-h + x, y)/h**2
        + f(h + x, y)/h**2
\end{lstlisting}
\end{minipage}

In the resulting stencil expression Devito inserts the symbol
\lstinline!h! to denote grid spacing in the \lstinline!x! dimension.
It is important to note here that this notation allows users to
quickly increase the order of the stencil discretization by simply
changing a single parameter during symbol construction.

Another utility property of Devito's symbolic data objects is the
\lstinline!f.laplace! operator that expands to \lstinline!(f.dx2 + f.dy2)!
for two-dimensional problems and \lstinline!(f.dx2 + f.dy2 + f.dz)! for
three-dimensional ones. This allows the dimensionless expression of
complex equations, such as a simple wave equation, to be defined very
concisely:

\begin{minipage}{\linewidth}
\begin{lstlisting}
In [1]: from devito import DenseData
In [2]: from sympy import Eq

In [3]: u = TimeData(name='u', space_order=2,
                     time_order=2, shape=(10, 12))
In [4]: m = DenseData(name='m', space_order=2
                      shape=(10, 12)))

In [5]: eqn = Eq(m * u.dt2, u.laplace)
Out[5]: Eq((-2*u(t, x, y)/s**2
            + u(-s + t, x, y)/s**2
            + u(s + t, x, y)/s**2)*m(x, y),
            -4*u(t, x, y)/h**2
            + u(t, x, -h + y)/h**2
            + u(t, x, h + y)/h**2
            + u(t, -h + x, y)/h**2
            + u(t, h + x, y)/h**2)
\end{lstlisting}
\end{minipage}

The most important abstraction provided by Devito however is the
\lstinline!Operator! class that allows users to trigger the automated
code generation and compilation cascade. Users simply need to provide an
expanded stencil expression and optional compiler settings via
\lstinline!op = Operator(eqn, compiler=IntelCompiler)!. The resulting
operator object can then be used to execute the generated kernel code
via \lstinline!op.apply(u, m)!.

\subsection{Low-level API: Expressions and loops}
\label{sec:devito_low}

One of the key transformations performed during the code generation
process is to resolve the grid spacing variables, for example
\lstinline!h!, into explicit grid indices. During this transformation
Devito converts symbolic functions like \lstinline!f(x, y)! into
\lstinline!sympy.Indexed! expressions of the form \lstinline!f[x, y]!,
before converting the dimension symbols to loop counters to yield
\lstinline!f[i1, i2]!. It is important to note here that the
intermediate indexed representation with explicit dimensions is also
accepted input when creating \lstinline!Operator! instances, allowing
users to input custom expressions with irregular grid indexing. The
low-level API is also a useful tool for developers to rapidly test
and develop additional higher-level abstractions.

Moreover, a new low-level API is currently under development that
will allow the automatic inference of iteration spaces and loop
constructs from the indexed stencil expression itself. The key
to this are \lstinline!Dimension! objects that act like index symbols
(\lstinline!x!, \lstinline!y! or \lstinline!t!) and define the
iteration space over which to loop in the generated kernel code, while
defining additional meta-information like boundary padding layers or
buffer switching for time loops. This new API will allow the generic
and automated construction of arbitrary loop nests via provided
\lstinline!Iteration! and \lstinline!Expression! classes, enabling the
rapid development of further mathematical abstractions in the
high-level API. This will also provide the possibility of closer
integration with dedicated access-execute frameworks (see
\secref{sec:dsl-abstraction-layers}) in the future to further leverage
automated loop-level performance optimizations.

\subsection{Performance optimization}
\label{sec:devito_performance}

Devito provides a set of automated performance optimizations during
code generation that allow user applications to fully utilise the
target hardware without changing the model specification. The
optimizations provided by Devito range from symbolic manipulations
applied to the provided stencil expressions to hardware-specific
implementation choices like loop blocking.

\subsubsection{Shared-memory Parallelism}

Devito provides thread-parallel execution through automated code
annotation with OpenMP pragmas inserted at the appropriate locations
in the code, with a single thread pool allocated outside of the
primary timestepping loop. On multi-socket systems Devito performs so
called first touch memory allocation on its symbolic data objects
ensuring high data locality on NUMA architectures.

\subsubsection{Vectorization}

SIMD parallelism and vectorization are a crucial performance
optimization for compute intensive stencil codes to fully exploit
modern HPC architectures. Devito ensures efficient vectorization by
inserting toolchain-specific compiler hints into the generated code
and enforcing data alignment on page boundaries during allocation.

\subsubsection{Loop Blocking}

Loop or cache blocking is a well-known performance optimization for
stencil codes that aims to increase data locality and thus cache
utilisation by splitting iteration spaces (loops) into spatially
aligned blocks. The effectiveness of the technique is highly dependent
on the cache size of the target architecture. Devito therefore
provides an interfaces for users to define block sizes explicitly, as
well as an ``auto-tune'' mode where a defined \lstinline!Operator!
object is used to identify the optimal block size through a
brute-force search.

\subsubsection{Common sub-expression elimination}

On top of architecture-specific optimizations Devito also aims to
leverage symbolic optimizations to speed up computation. An important
example of this is common sub-expression elimination (CSE) that aims
to reduce the amount of redundant computation within complex stencil
expressions by factoring out common terms. While it may be argued that
most modern optimizing compilers are able to do this, it is important
to note here that Devito utilises SymPy{}'s internal CSE functionality
to apply the technique at a much higher abstraction level.

As a result, processing times of symbolic expressions are reduced
during the code generation phase, and the compilation time of the
generated code itself is improved. This is important since complex
stencils can increase compilation time to several hours if CSE is
performed by the compiler, in particular for high-order seismic
inversion kernels that often contain very large numbers of terms.

\section{Discussion}

In this paper we present the Devito finite difference DSL and showcase
its multi-level Python API to efficiently generate complex
mathematical operators from SymPy{} expressions. We highlight Devito's
practical applicability in creating forward and adjoint operators for
seismic inversion operators from concise symbolic definitions and
verify their correctness by showing that the generated code does in
fact behave mathematically like the adjoint of the equation. Moreover,
we demonstrate that the optimized code generated by Devito is capable
of achieving a significant percentage of peak performance on a
traditional CPU, as well as a \Intel Xeon Phi\texttrademark  architecture.

When showcasing Devito's symbolic Python API we also compare it to
traditional approaches to accelerating Python computation and show
significant speed-ups on even a single core CPU. On top of the evident
performance benefits, a more concise formulation of simple
mathematical operators can be achieved using Devito, which increases
the productivity of domain scientists. Conversely, our total reliance
on the SymPy{} package demonstrates the importance of the Python
programming language and its vast software ecosystem in creating
performance-oriented DSLs and utilising them to build actual
scientific applications.

The use of Devito's second-layer API to create custom point-grid
interpolation routines further illustrates the complexities that often
prevent stencil DSLs from being used in practice. In the context of
seismic inversion, sparse point interpolation is the most important
feature that other finite difference DSLs lack. In addition, Devito 's
full compatibility with other scientific Python packages through the
NumPy{} array API entails that it is easy to implement a full-scale
inversion workflow in Python - including additional features such as
MPI parallelism between multiple inversion operators.

\section{Future Work}

Even though the primary target application for Devito is seismic
inversion, we aim to develop Devito towards becoming a generic finite
difference DSL. The diffusion example presented in
\secref{sec:diffusion_devito} shows that the current abstraction
generalises well to problems beyond the wave equation and its
additional low-level API provides an important degree of flexibility
when adding new features. Additional finite difference mechanisms,
such as staggered grids and new boundary condition types are of great
importance here to generalise Devito's applicability to other
application areas.

Additional improvements to the high and low-level symbolic APIs are
also planned that include customizable dimension symbols and the
automated derivation of iteration spaces. This additional automation
will in turn enable further integration with optimising compilers and
other external stencil frameworks and potentially even allow the
automated use of polyhedral compilation tools.

The seismic inversion operators demonstrated in this paper, and in
particular the sparse point interpolation features discussed in
\secref{sec:interpolation}, are still in an early development
stage. The seismic inversion example presented here uses the acoustic
isotropic wave equation, although, for real world applications, more
complex physical models involving anisotropy and elastic wave
propagation are commonly seen. In the future we plan to implement
these complex physical models using Devito, expanding Devito's
repertoire of features wherever necessary.

\section{Seismic Inversion Example}

The key motivating application behind the creation of the Devito
finite fifference DSL are seismic exploration problems including
advanced modelling and adjoint-based inversion. These not only require
a fast and accurate wave propagation model, but also efficient,
easy-to-define and rigorous adjoint operators. In the following
example we demonstrate the symbolic definition of modelling and
adjoint modelling operations for the acoustic wave equation denoted
as:
\begin{equation} \label{eqn:acoustic}
  m \frac{\partial^2 u}{\partial t^2} + \eta \frac{\partial u}{\partial t} - \nabla^2 u = q
\end{equation}
where $u$ denotes the pressure wave field, $m$ is the square slowness,
$q$ is the source term and $\eta$ denotes the spatially varying
dampening factor used to implement an absorbing boundary condition.

\subsection{Acoustic wave equation}
\label{sec:acoustic}

The implementation of the forward operator with an explicit time marching
scheme is shown in \lstref{lst:adjoint}. First, the symbolic data
objects are created, where \lstinline!u! is defined as a time-varying
function, whereas \lstinline!m! and \lstinline[mathescape]!$\eta$! are
spatial grids only. Notably, the source term \lstinline!q! is missing
from the high-level PDE description, since the injection of the source
requires sparse point-to-grid interpolation. Since interpolation does
not fall within a finite difference abstraction, it is dealt with via
separate expressions, as further explained in
\secref{sec:interpolation}

It is also important to note here that the spatial order of the
derivative stencils in the operator definition is provided as an
argument and can thus be changed easily by the user, and that the
entire formulation is independent of problem dimension. From the
individual component symbols we create the wave equation in a purely
symbolic form, leaving Devito's \lstinline!dt!, \lstinline!dt2!  and
\lstinline!laplace! operators to perform the stencil expansion
according to the spatial and temporal discretization order. The
resulting expression is then rearranged using the SymPy{}
\lstinline!solve!  routine to denote forward propagation, from which
the corresponding \lstinline!Operator! can be created and executed.

\begin{lstlisting}[
    caption=Example code to solve the wave equation,
    label={lst:wave}, mathescape=true, float, floatplacement=HT]
def forward(model, nt, dt, h, order=2):
    shape = model.shape
    m = DenseData(name="m", shape=shape,
                  space_order=order)
    m.data[:] = model  # Set model data
    u = TimeData(name='u', shape=model.shape,
                 time_dim=nt, time_order=2,
                 space_order=order)
    $\eta$ = DenseData(name="$\eta$", shape=shape,
                  space_order=order)

    # Derive stencil from symbolic equation
    eqn = m * u.dt2 - u.laplace + $\eta$ * u.dt
    stencil = solve(eqn, u.forward)[0]

    op = Operator(stencils=Eq(u.forward, stencil),
                  nt=nt, subs={s: dt, h: h},
                  shape=shape, forward=True)
    # Source injection code omitted for brevity

    op.apply()
\end{lstlisting}

A very similar symbolic definition can be used to define the adjoint
operator for the wave equation, as demonstrated in Listing
\ref{lst:adjoint}. The key difference to the forward implementation
is that the dampening term \lstinline[mathescape]! $\eta$ * u.dt! is subtracted rather
than added, and that we rearrange the stencil expression for execution
backward in time. The shorthand notation \lstinline!u.forward! and
\lstinline!u.backward! hereby denote the highest and lowest stencil
point in the second-order time discretization stencil,
\lstinline!t + s! and \lstinline!t - s! respectively. These can be
used to ensure the correct time marching direction once the symbolic
function offsets \lstinline!s! and \lstinline!h! have been resolved to
explicit grid indices in the final stencil expression.

\begin{lstlisting}[
    caption=Example code to solve the adjoint wave equation,
    label={lst:adjoint}, mathescape=true, float, floatplacement=HT]
def adjoint(model, nt, dt, h, spc_order=2):
    shape = model.shape
    m = DenseData(name="m", shape=shape,
                  space_order=order)
    m.data[:] = model
    v = TimeData(name='v', shape=shape,
                 time_dim=nt, time_order=2,
                 space_order=spc_order)
    $\eta$ = DenseData(name="$\eta$", shape=shape,
                  space_order=order)

    # Derive stencil from symbolic equation
    eqn = m * v.dt2 - v.laplace - $\eta$ * v.dt
    stencil = solve(eqn, v.backward)[0]

    op = Operator(stencils=Eq(u.backward, stencil),
                  nt=nt, subs={s: dt, h: h},
                  shape=shape, forward=False)
    # Receiver interpolation omitted for brevity

    op.apply()
\end{lstlisting}

\subsection{Sparse point interpolation}
\label{sec:interpolation}

One key restriction of the high-level PDE notation used in the
previous examples is the inability to express interpolation between
irregularly spaced sparse points and the grid. This feature is
frequently required to inject source terms by
point-to-grid interpolation and sample the resulting wave field at
non-aligned sparse points via grid-to-point interpolation.

\begin{lstlisting}[
    caption=Example code for point-to-grid interpolation.,
    label={lst:point2grid}, float, floatplacement=HT]
eqns = [Eq(u[t, x, y], u[t, x, y] + expr),
        Eq(u[t, x+1, y], u[t, x+1, y] + expr),
        Eq(u[t, x, y+1], u[t, x, y+1] + expr),
        Eq(u[t, x+1, y+1], u[t, x+1, y+1] + expr)]
it = Iteration(eqns, index=p, limits=source.shape)]
\end{lstlisting}

\lstref{lst:point2grid} demonstrates how to construct a set of
expressions to perform point-to-grid interpolation, where each grid
value enclosing the sparse point is updated according to a source
term. Most notably, the update to the grid value in the wave field
\lstinline!u! is created using the low-level indexed variable format,
allowing us to define the explicit offsets needed to update every
enclosing grid point with the relevant source expression
\lstinline!expr!. A similar set of expressions can be generated to
perform the inverse interpolation from grid values to sparse points,
but an example of this is omitted here for brevity. Grid-to-point
interpolation, however, is of particular importance for seismic
inversion problems, since the behaviour of sparse points changes
between the forward and the adjoint runs.

Since there are potentially many sources in a wave field model, we can
also construct an additional loop around the expression set to iterate
over all sparse points contained in \lstinline!source!. For this,
Devito provides an additional low-level \lstinline!Iteration! class
that adds a user-defined loop construct into the auto-generated code.
This enables custom iterations over source source points and their
associated coordinate data, allowing users to create
coordinate-dependent symbolic expressions that define the interpolation
for each reference cell point.

\subsection{Validation}

We pointed out the requirement to have rigorous adjoint in order to
use our framework in seismic inversion. To validate our implementation
we use the simple definition of an adjoint as a test. For a given
matrix $A$, its adjoint is defined as the matrix $A^T$ such that:
\begin{eqnarray} \label{eqn:adjtest}
	&\forall (x, y),\nonumber\\
	&< Ax, y> - <x, A^Ty> = 0,\\
	&\frac{< Ax, y>}{<x, A^Ty>} = 1,\nonumber\\
\end{eqnarray}

where $<.,.>$ denotes the inner product. In the case of the wave
equation, $A$ represent the modelling operator while $A^T$ represents
the adjoint modelling operator, and from the definition in Equation
\ref{eqn:adjtest} we can compare the output of our kernels to verify
the accuracy of the adjoint. We validate our implementation with this
test in a two-dimensional and three-dimensional setting for multiple
spacial discretization order in Table \ref{adjres} .

\begin{table*}
\centering
\caption{Adjoint test for different discretization orders in 2D and
  3D. The test is computed on a two layer model.}
\label{adjres}
\begin{tabular}{|c|c|c||c|c|c|}
Order & Dimension & $<Ax,y>$ & $<x,A^Ty>$ &
           Difference & ratio \\
\hline
2nd order      &  2D:  & 373323.7976& 373323.7975434& 6.07169350e-05  &  1.0\\
4th order      &  2D:  & 340158.1486& 340158.1485252&  0.00012756     &  1.0\\
6th order      &  2D:  & 341557.3948& 341557.3947399&  0.00014287     &  1.0\\
8th order      &  2D:  & 358240.8513& 358240.8511931&  0.00016741     &  1.0\\
10th order     &  2D:  & 393488.5561& 393488.5559269&  0.00023841     &  1.0\\
12th order     &  2D:  & 439561.4005& 439561.4002034&  0.00035794     &  1.0\\
2nd order      &  3D:  & 2.17496552 & 2.174965534979&  -1.23030883e-08&  0.99999999\\
4th order      &  3D:  & 3.64447937 & 3.644479393901&  -2.13132316e-08&  0.99999999\\
6th order      &  3D:  & 3.78730372 & 3.787303745086&  -2.22477072e-08&  0.99999999\\
8th order      &  3D:  & 3.80286229 & 3.802862312852&  -2.23545817e-08&  0.99999999\\
10th order     &  3D:  & 3.80557957 & 3.805579596334&  -2.23736993e-08&  0.99999999\\
12th order     &  3D:  & 3.80318675 & 3.803186773144&  -2.23587757e-08&  0.99999999\\
\end{tabular}
\end{table*}

\subsection{Performance}
\label{sec:acoustic_performance}

The performance of the finite difference code generated by Devito is
demonstrated on two target architectures, a \Intel Xeon\texttrademark E5-2690v2
architecture with 10 physical cores running at 3GHz and a \Intel Xeon Phi\texttrademark
accelerator card, in \figref{fig:acoustic_performance}. The model used
for the benchmark is a gradient test with forward and adjoint
operators on a three-dimensional grid with $201 \times 201 \times 70$
grid points and $40$ PML grid points on each side, resulting in a
computational grid of size $281 \times 281 \times 150$. The grid size
is $15 m$ and the source term is a Ricker wavelet at $10 Hz$. The wave
field is modelled for 1 second with spatial discretizations of varying
order from 2 to 16. The results indicate that Devito is able to utilise
both architectures to a high degree of efficiency, while maintaining
the ability to increase accuracy by switching to higher order stencil
discretizations dynamically.

\begin{figure}[ht] \centering
\begin{subfigure}{.5\textwidth} \centering
\includegraphics[width=\linewidth]{./figures/xeonorder.pdf}
\caption{\Intel Xeon\texttrademark E5-2690v2 10C @ 3GHz}
\end{subfigure}

\begin{subfigure}{.5\textwidth} \centering
\includegraphics[width=\linewidth]{./figures/xeonphiorder.pdf}
\caption{\Intel Xeon Phi\texttrademark }
\end{subfigure}

\caption{Roofline plots for two target systems}
\label{fig:acoustic_performance}
\end{figure}

\section{2D Diffusion Example}

The use of symbolic notations to define finite difference operators
strongly enhances the efficiency of model developers and allows well
known problems to be expressed in a concise manner. In this section we
will first demonstrate the use of symbolic Python notation in comparison
to pure Python code on a simple two-dimensional diffusion example. We
will also analyse the performance benefits of Devito operators over
vectorized NumPy{} array notation, before demonstrating the
construction of seismic inversion operators using Devito's symbolic
interfaces in the next section.

The general diffusion or heat equation in two dimensions is commonly
notated as
\begin{equation} \label{eqn:diffusion}
  \frac{\partial u}{\partial t} = \alpha\nabla^2 u
\end{equation}
where $u$ denotes the unknown function to solve for, $\alpha$ is the
diffusion coefficient and $\nabla^2u$ denotes the Laplace operator, as
defined in \secref{sec:devito_high}. To keep the following examples
short we will be using a forward Euler approach with constant grid
spacing in $x$ and $y$, while the the timestep size is computed as
\begin{equation} \label{eqn:diffusion_dt}
\Delta t = \frac{\Delta x^2\Delta y^2}{2\alpha (\Delta x^2 + \Delta
  y^2)}
\end{equation}
where $\Delta x$, $\Delta y$ denotes the grid spacing and $\Delta t$
is the timestep size.

\subsection{Pure Python}

A first direct implementation of the described equation using an
alternating buffer scheme for \lstinline!u! in pure Python is shown in
\lstref{lst:diff_python}. The implementation consists of three loops
that apply a 5-point star stencil - a 3-point stencil in \lstinline!x!
(\lstinline!uxx!) and a 3-point stencil in \lstinline!y!
(\lstinline!uyy!) to compute the update. It is important to note here
that the stencil is only applied to elements in the interior of the
grid in order to avoid out-of-bound array accesses.

\begin{lstlisting}[
    caption= Diffusion code using pure Python loops.,
    label={lst:diff_python}, float, floatplacement=HT]
_, nx, ny = u.shape
for t in range(timesteps):
    t0 = t
    t1 = (t + 1)
    for i in range(1, nx-1):
        for j in range(1, ny-1):
            uxx = (u[t0, i+1, j] - 2*u[t0, i, j]
                   + u[t0, i-1, j]) / dx2
            uyy = (u[t0, i, j+1] - 2*u[t0, i, j]
                   + u[t0, i, j-1]) / dy2
            u[t1, i, j] = u[t0, i, j]
                          + dt * a * (uxx + uyy)
\end{lstlisting}

\subsection{Vectorized NumPy{} arrays}

Since the loop performance overhead of the Python interpreter is well
known, the most common optimization of simple arithmetic computations
on array and grid data is to utilise NumPy{} vectorization. This
applies arithmetic operations element-wise over all elements in an
array using pre-compiled kernels. For stencil expressions it is
therefore important to adjust the bounds of the array view that
represents each stencil point, as demonstrated in
\lstref{lst:diff_numpy}.

\begin{lstlisting}[
    caption=Diffusion code with vectorized NumPy{} arrays.,
    label={lst:diff_numpy}, float,floatplacement=HT]
for t in range(timesteps):
    t0 = t
    t1 = (t + 1)
    uxx = (u[t0, 2:, 1:-1] - 2*u[t0, 1:-1, 1:-1]
           + u[t0, :-2, 1:-1]) / dx2
    uyy = (u[t0, 1:-1, 2:] - 2*u[t0, 1:-1, 1:-1]
           + u[t0, 1:-1, :-2]) / dy2
    u[t1, 1:-1, 1:-1] = u[t0, 1:-1, 1:-1]
                        + a * dt * (uxx + uyy)
\end{lstlisting}

\subsection{Symbolic Python - lambdified}

An alternative approach to specifying the loop kernel that implements
a particular equation directly is to define it symbolically via
SymPy{}, as demonstrated in \lstref{lst:diff_sympy}. The diffusion
equation can thus be expressed in a manner that is mathematically
concise, while SymPy{} utilities are used to expand the derivative
stencils and create an executable function via \lstinline!sympy.lambdify!.
It is noteworthy that the \lstinline!'numpy'! mode argument allows us to utilise
vectorized operations on sub-arrays, similar to the pure NumPy{}
variant.

\begin{lstlisting}[
    caption=Diffusion code using symbolic SymPy objects.,
    label={lst:diff_sympy}, float, floatplacement=HT]
from sympy import Eq, Function, lambdify
from sympy.abc import x, y, t, s, h
p = Function('p')
dx2 = as_finite_diff(p(x, y, t).diff(x, x),
                     [x - h, x, x + h])
dy2 = as_finite_diff(p(x, y, t).diff(y, y),
                     [y - h, y, y + h])
dt = as_finite_diff(p(x, y, t).diff(t), [t, t+s])
eqn = Eq(dt, a * (dx2 + dy2))
stencil = solve(eqn, p(x, y, t + s))[0]
subs = (p(x, y, t), p(x+h, y, t), p(x-h, y, t),
        p(x, y+h, t), p(x, y-h, t), s, h, a)
kernel = lambdify(subs, stencil, 'numpy')

for ti in range(timesteps):
    t0 = ti
    t1 = (ti + 1)
    u[t1, 1:-1, 1:-1] = kernel(
        u[t0, 1:-1, 1:-1], u[t0, 2:, 1:-1],
        u[t0, :-2, 1:-1], u[t0, 1:-1, 2:]
        u[t0, 1:-1, :-2], dt, dx)
\end{lstlisting}

\subsection{Symbolic Python - Devito}
\label{sec:diffusion_devito}
A similar symbolic approach is provided by Devito, where the expansion
of the derivative stencils is automated by the provided symbolic data
objects. As shown in \lstref{lst:diff_devito}, this allows the stencil
discretization to be adjusted via a single parameter change, while the
equation definition and reorganisation of the stencil expression via
\lstinline!sympy.solve! is identical to the pure SymPy{} approach.
The \lstinline!devito.Operator! created from the expression is then
used to generate low-level code and execute it for a given number of
timesteps.

\begin{lstlisting}[
    caption=Diffusion code using the symbolic Devito API.,
    label={lst:diff_devito}, float, floatplacement=HT]
from devito import TimeData, Operator
from sympy.abc import s, h
u = TimeData(name='u', shape=(nx, ny),
             time_order=1, space_order=2)
u.data[0, :] = ui[:]

eqn = Eq(u.dt, a * (u.dx2 + u.dy2))
stencil = solve(eqn, u.forward)[0]
op = Operator(stencils=Eq(u.forward, stencil),
              subs={h: dx, s: dt}, nt=timesteps)
op.apply()
\end{lstlisting}

An example of the auto-generated C code, is shown in
\figref{lst:diff_devito_c}. The example demonstrates Devito's use of
parallelisation and vectorization pragmas, as well as its internal
resolution of alternating buffer access. Moreover it can be seen that
Devito has replaced all numerical constants, such as the diffusion
coefficient \lstinline!a! and spacing parameters \lstinline!dx! and
\lstinline!dt!, and inserted all known loop dimensions to aid compiler
optimizations.

\begin{lstlisting}[
    basicstyle={\scriptsize\ttfamily},
    caption=Auto-generated C code to solve diffusion example.,
    label={lst:diff_devito_c}, float, floatplacement=HT]
extern ``C'' int Operator(float *u_vec)
{
  float (*u)[1000][1000] = (float (*)[1000][1000]) u_vec;
  {
    int t0;
    int t1;
    #pragma omp parallel
    for (int i3 = 0; i3<500; i3+=1)
    {
      #pragma omp single
      {
        t0 = (i3
        t1 = (t0 + 1
      }
      {
        #pragma omp for schedule(static)
        for (int i1 = 1; i1<999; i1++)
        {
          #pragma omp simd aligned(u:64)
          for (int i2 = 1; i2<999; i2++)
          {
            u[t1][i1][i2] = 2.5e-1F*u[t0][i1][i2-1]
                          + 2.5e-1F*u[t0][i1][i2+1]
                          + 2.5e-1F*u[t0][i1-1][i2]
                          + 2.5e-1F*u[t0][i1+1][i2];
          }
        }
      }
    }
  }
\end{lstlisting}

\subsection{Performance comparison}

To demonstrate the performance advantages of custom code generation
and JIT compilation as performed in Devito over traditional approaches
to Python acceleration, we have compared the single core performance
of the vectorized NumPy{} and SymPy{} implementations of the diffusion
example. For this purpose the diffusion benchmark was run for 500
timesteps with a grid size of 0.001 on a single core of a Intel\textsuperscript{\tiny{\textregistered}} Core
i3-4030U CPU. The results presented in
\figref{fig:diffusion_performance} show a clear performance advantage
of the DSL approach, as well as highlighting a significant overhead of
the vectorized SymPy{} implementation over the pure NumPy{}
approach. It is important to note here that the Devito DSL is the only
benchmarked approach that also encapsulates the timestepping loop
in the compilation step and that further performance advantages can be
expected when shared-memory parallelism is utilised, as we will show in
\secref{sec:acoustic_performance}.

\begin{figure}[ht] \centering
\includegraphics[width=.4\textwidth]{./figures/DiffusionModes.pdf}
\caption{Serial runtime comparison between NumPy{}, SymPy{} and
  Devito for the diffusion example.}
\label{fig:diffusion_performance}
\end{figure}

\section{Introduction}

Python is one of the most popular high-level programming languages for
academics and scientists due to its clean syntax, great
expressiveness and a vast ecosystem of open source
packages~\cite{Oliphant2007,Millman2011,Perez2011}. Its direct use in
High Performance Computing (HPC), though, is complicated by the slow
execution within the Python interpreter itself. Several mechanisms
exist to overcome the interpreter overhead, ranging from static
compilation via Cython~\cite{Behnel2011} to generalised Just-In-Time
(JIT) compilation via Numba~\cite{Lam2015}.

However, for specific, often compute-intensive application areas such as
solving partial differential equations (PDE) with the finite element
method, domain specific languages (DSL) have also proven
successful~\cite{Logg2012}. DSLs aim to decouple the problem
specification from its low-level implementation to create a
separation of concerns between domain scientists and HPC specialists,
resulting in a direct payoff in productivity. Python is particularly
useful for creating DSLs, since its runtime type inference and
operator overloading via ``magic functions'' allows domain-specific
constructs to be easily embedded in the language itself.

Moreover, due to their restriction to a single problem domain,
code-generating DSLs can also augment the generated code to a
particular computational architecture, allowing hardware-specific
performance optimizations. In the age of specialized HPC
architectures, such as \Intel Xeon Phi\texttrademark  or GPU accelerators, this not only
promises performance portability but also provides a valuable
alternative to the continuing search for new all-encompassing parallel
programming paradigms.

In this paper we present
Devito\footnote{\url{https://github.com/opesci/devito}}, a new finite
fifference DSL that utilises the SymPy{}~\cite{Meuer2016} package to
generate optimized parallel C code from high-level symbolic operator
definitions. Devito is primarily targeted at generating fast wave
propagation kernels for seismic inversion problems and therefore
offers an abstraction hierarchy that enables not only the generation
of fast stencil kernels, but also a range of domain-specific features,
such as sparse point interpolation, required in a real scientific
context. For this purpose Devito provides a two-level API that allows
users to define operator kernels as a mixture of high-level symbolic
equations and explicit C-like SymPy{} expressions to solve complex
inversion problems.