\section{Introduction}
\label{introduction}

The VC dimension is a measure of what is called the ``capacity'' of a classification algorithm: that is, roughly speaking, its flexibility or expressive power.  The practical interest of this quantity arises from its role in the \emph{Probably Approximately Correct} learning model and related models (see e.g. \citet[ch. 3]{Shais}, \citet{Holden95}), where it appears in results relating the size of the training set to the classifier's accuracy.  The VC dimension is defined in terms of the concept of \emph{shattering}~\citep{Vapnik98}. For a two-class problem in $d$-dimensional space $\mathbb{R}^d$, a set $A$ of functions $\mathbb{R}^d \to \{+1,-1\}$ is said to shatter a set of points (in $\mathbb{R}^d$) if any labelling of that set can be given by an element of $A$.  The VC dimension of $A$ is the size of the largest set which can be shattered by $A$.

The nearest-neighbour (1NN) classification rule is a classic technique of supervised learning, first introduced by \cite{Fix52}.  The 1NN rule determines a label for (``classifies'') a given point in a metric space by assigning it the label of the nearest point in a previously determined set of labelled reference points. In the original and simplest algorithm, the reference set is the set of all the training data.  If the amount of training data is \emph{a priori} unbounded, then the VC dimension of the set of associated 1NN-rule classifying functions is infinite: trivially, a set of points of any size is labelled correctly by a classifier whose reference set is the same set of labelled points.

However, it is often not practical to store all the training data, especially in an era of ``big data''.  Therefore  many algorithms have been proposed which learn a smaller reference set from the training data: \cite{Garcia12} and \cite{Triguero12} survey more than 75 such algorithms between them.  What is the VC dimension associated with these ``editing'' NN algorithms?  If the reference set may grow without bound, then the VC dimension is infinite, as for the na{\"i}ve classifier described above.  But if the reference set is constrained not to exceed a given size, then the VC dimension is finite. The main results of the present work are lower and upper bounds for the VC dimension of the set of all 1NN-rule classifiers which use a reference set of given fixed size, in Euclidean space.

It should be noted that, in general, the size of the reference set to be formed by an editing algorithm is not fixed \textit{a priori}.  There are exceptions, an important one being the case where a classifier for a data stream is kept current by using a fixed window of the most recent $N$ points as its reference set \citep{GunnNeurocomputing}.  But it is very common to impose a \emph{maximum} size limit on the reference set, if only by hand in an \textit{ad hoc} fashion.  Our upper bounds will apply to any algorithm for which the size of the reference set is bounded, though our lower bounds will apply only where the size is prescribed, and the algorithm is capable of returning any reference set of that size.

We believe that there is significant interest in the VC dimension of prototype classifiers with  fixed-size reference sets, and that this is shown by the fact that an incorrect purported result \cite[Proposition 2]{Karacali} has been cited more than 50 times as giving the VC dimension of such classifier sets.  Our study corrects the record.

We feel there is significant value in bringing together the existing theoretical results, some of which seem otherwise likely to remain obscure to practitioners.  Beyond our deductions from existing theory, our novel contributions are 1) A new lower bound for the two-dimensional case, higher than that implied by previous results for polytopes (Proposition \ref{OneBetter}),  and 2) an upper bound for all dimensions, slightly less tight than the best that can be deduced from existing theory, but avoiding the use of exotic functions and thereby facilitating a discussion of the asymptotic behaviour of the limits (Corollary \ref{LooseCor}).

We will discuss the theoretical framework in section \ref{NNinVC}, and briefly review the relevant literature in section \ref{lit}. In sections \ref{lower2d} and \ref{lowerdd} we derive lower bounds for the VC dimension, and in section \ref{upper} we determine upper bounds.  Our results are summarised in section \ref{conclusions}.

\section{1NN classifiers in the VC theory}
\label{NNinVC}

We use ``classifier'' to mean a function from the feature space to the set of classes, which is used to classify unlabelled examples. An algorithm which learns such a function from training data is a ``classification algorithm''. The set of all classifiers which the algorithm might produce in response to all sets of training data is called the ``hypothesis class'' of that algorithm.  (For example, the original Rosenblatt Perceptron selects among all functions which map one half-space to one class, and the other half-space to the other class: this set of functions is the hypothesis class of the Perceptron.)  When we speak of the VC dimension of a classification algorithm, we mean the VC dimension of the hypothesis class from which that algorithm learns a classifier.

NN classification algorithms
are not usually thought of as selecting a classifier from a fixed hypothesis class in this way~\citep[see e.g.][ch.\ 19]{Shais}.
This is because the great practical advantage of 1NN classification is that a classifier function does not need to be explicitly evaluated when classifying an unlabelled example; the new point is simply assigned the class of the nearest prototype in the reference set, which can be identified efficiently using a $k$-d tree.  However, this process is equivalent to classifying the new point according to a classifier function which maps the Voronoi cell of each prototype to the label of that prototype.  Thus, the hypothesis class of the 1NN rule with a reference set of $m$ prototypes in $d$-dimensional space is the set of all labellings of all $m$-cell Voronoi diagrams in the space; it is parameterised by the $md$ co-ordinates of the $m$ prototypes, and the $m$ choices of label.

We consider only features which take real values, thus classifiers whose domain is $\mathbb{R}^d$, and consider only the 1NN rule with the Euclidean metric. $\mathrm{1NN}(d,m)$ will denote the set of all classifiers $g:\mathbb{R}^d \to \{+1,-1\}$ which use the (Euclidean) nearest-neighbour rule with a reference set of size $m$.  $\mathrm{VCdim}(A)$ will denote the VC dimension of a set $A$ of classifiers.  The purpose of this paper is to give lower and upper bounds for $\mathrm{VCdim}(\mathrm{1NN}(d,m))$.

\section{Related Work}
\label{lit}

A number of existing theoretical results have implications for the VC dimension of the NN classifier with arbitrary reference set of fixed size.  However, these results are scattered in the literature and their implications have not been made explicit. In particular, the results we use to establish lower bounds were developed for a class of polytope classifiers, without mention of NN classifiers. In the present section we give a brief overview of relevant theoretical work, starting with some brief historical context and going on to include the work whose implications we will directly use in subsequent sections.

Questions of ``separating capacities'' of families of decision surfaces were already considered before the advent of the Vapnik-Chervonenkis theory (see e.g.\ \cite{Cover65} and references therein).  During and shortly after the years of the initial development of the VC theory, several authors used approaches from classical combinatorial geometry to derive results about the VC dimension (or related separability properties) of several simple sets.  Readers interested in this literature will need to be aware of the distinction between the VC dimension of a set of classifiers, as we have defined it above, and the VC dimension of a set of subsets of the Euclidean space in question (called by some authors a ``concept class'').  See our discussion at the start of section \ref{upper}, and~\citet[pp. 196, 199, and 215]{Devroye96}.

To give an example of the results obtained, \citet{Dudley79} reports that the set of balls in $\mathbb{R}^d$ has VC dimension $d+1$.  (N.B. Dudley's quantity $V$ is one greater than the quantity defined as the VC dimension in more recent literature.)  Similarly, it may be shown that the set of all half-spaces in $\mathbb{R}^d$ has VC dimension $d+1$~\citep[see e.g.][ch. 13]{Devroye96}.

More recent authors have considered the intersection or union of half-spaces, which is to say, sets which are the interior or exterior of (possibly unbounded) polytopes. \citet{Blumer89} show that the set of interiors of $N$-gons has VC dimension $2N + 1$.  This result is quoted by \citet{TakacsPataki} as the starting point for their work in which they find upper and lower bounds for the VC dimension of convex polytope classifiers, from which we will derive a lower bound in section~\ref{lowerdd}.

NN classifiers, like convex polytope classifiers, have decision boundaries which are the union of subsets of hyperplanes.  However, the decision region for a NN classifier is not in general a simple intersection or union of half-spaces; it may be a complicated union of the interiors of polytopes formed by such intersections.  This may explain why the question of the VC dimension of the NN classifier with an arbitrary reference set of fixed size $m$ has not previously, to our knowledge, been addressed, other than in~\cite{Karacali}.  (\citet{Devroye96} consider the closely related property of the \emph{shatter coefficient} of such classifiers; we will make explicit the implications of this work in section \ref{upper}.)

\citet[Proposition 2]{Karacali} claim that the VC dimension of the NN classifier with reference set of size $m$, $\mathrm{VCdim}(\mathrm{1NN}(d,m))$ in our notation, is exactly $m$. Their argument consists of exhibiting a set of $m+1$ points which cannot be correctly classified by $m$ prototypes. They do not argue that a general set of $m+1$ points cannot be shattered.  This apparently reflects a misunderstanding of the definition of VC dimension.  The VC dimension is the largest number for which some set of that size can be shattered, not the largest number such that all sets of that size can be shattered.  The latter quantity is called the Popper dimension, a quantity
which thus far has not found a r\^{o}le in statistical learning theory~\citep{Corfield09}.

\section{Lower bounds - two dimensions}
\label{lower2d}
Tak{\'a}cs and Pataki \citep{Takacs,TakacsPataki} prove bounds for the VC dimension of sets of classifiers whose decision boundaries are convex polytopes.  These sets can easily be related to sets of nearest-neighbour classifiers, giving lower bounds for the latter.  This is because a decision boundary which is an $N$-faceted convex polytope can be obtained as the decision boundary of a 1NN classifier with $N+1$ prototypes, by placing a prototype of one label inside the polytope, and the remaining $N$ prototypes, with the opposite label, as the reflections of the first prototype in the $N$ facets of the desired decision boundary.  See Figure \ref{fig:PolygonVoronoi} for an example construction with $N=5$, $d=2$. We formalise this observation as Proposition \ref{subset} below, and give a formal proof in the appendices.

\begin{figure}[htb]
 \centering
 \includegraphics[width=0.6\linewidth]{PolygonVoronoiRegion}
 \caption{Illustration of the construction of a data set whose classification region is determined by a given convex polygon $P$.  The required prototype set (filled and empty circles) is constructed by reflecting an arbitrary interior point in (the lines containing) the edges of $P$. Voronoi boundaries are shown; the Voronoi cell of the interior point is the decision region of the classifier, coincident with $P$.}
 \label{fig:PolygonVoronoi}
 \end{figure}

We will use $G(d,N)$ to denote the set of classifiers $g:\mathbb{R}^d \to \{+1,-1\}$ whose decision boundary is a convex $N$-faceted polytope.

\begin{lemma}
 \label{subset}
 The set of convex $N$-faceted polytope classifiers is a subset of the set of NN classifiers with reference set of size $N+1$, in the same Euclidean space.  That is,
  \begin{equation}
 G(d,N) \subseteq \mathrm{1NN}(d,N+1).
 \end{equation}
\end{lemma}

\begin{corollary}
\label{impliesLowerBound}
A lower bound for the VC dimension of $G(d,N)$ is also a lower bound for the VC dimension of $\mathrm{1NN}(d,N+1)$.
\end{corollary}

The corollary follows immediately: any set of points which can be shattered by an element of $G(d,N)$ can be shattered by an element (the same element) of $\mathrm{1NN}(d,N+1)$.

Tak{\'a}cs gives the following result for two-dimensional Euclidean space:

\begin{lemma}
 (Tak{\'a}cs)
 \label{Takacs1}

 \begin{equation}
 h(G(2,N)) \geq 2N + 2,
 \end{equation}
 for $N \geq 2$.
\end{lemma}

We give a sketch of the proof given in \cite{Takacs} in the appendices.  By a more complicated argument, Tak{\'a}cs establishes that this lower bound is also an upper bound for the VC dimension of polytope classifiers.  In general, upper bounds for polytope classifiers are of no relevance to the more general set of NN-rule classifiers.  But the $N=2$ case is an exception, and we make the following brief side remark about this case:

\begin{remark}
With three prototypes, the only decision surfaces an NN-rule classifier can form are an open 2-gon, or a pair of parallel lines.  But any finite set of points which can be correctly dichotomised by a pair of parallel lines can also be correctly split by a digon formed by a suitably small adjustment of the lines such that they are non-parallel.  So for this case, the set of NN-rule classifiers has no greater separating power than the related set of polygon classifiers.  Thus a precise value for the VC dimension of this set of NN-rule classifiers is established:

 \begin{equation}
 \mathrm{VCdim}(\mathrm{1NN}(2,3)) = 6.
 \end{equation}

\end{remark}

Now, for general $m$, Lemma \ref{Takacs1}, with Corollary \ref{impliesLowerBound}, implies a lower bound for $\mathrm{VCdim}(\mathrm{1NN}(2,m))$:

\begin{equation}
 \mathrm{VCdim}(\mathrm{1NN}(2,m)) \geq 2m
\end{equation}

However, we can do better than this.  For $m \geq 4$, the NN-rule classifier can create a larger class of decision surfaces than a single convex polytope. In Appendix~\ref{appendix1} we present a new argument, inspired by the elementary geometrical approach of Tak{\'a}cs but considerably more involved, demonstrating a stronger lower bound for $\mathrm{VCdim}(\mathrm{1NN}(2,m))$, $m \geq 4$:

\begin{proposition}
\label{OneBetter}
\begin{equation}
 \mathrm{VCdim}(\mathrm{1NN}(2,m)) \geq 2m + 1.
\end{equation}
\end{proposition}

Though this improvement is the smallest possible, it establishes the principle that the VC dimension of 1NN classifiers with $m$ prototypes in $\mathbb{R}^2$ is larger than that of the relevant comparable class of polygon decision boundaries (recall that Tak{\'a}cs' lower bound, Lemma \ref{Takacs1}, is also an upper bound for that class).  That is, Proposition \ref{OneBetter} establishes that 1NN classifiers in $\mathbb{R}^2$ gain in expressivity from their ability (for $m \geq 4$) to form boundaries other than convex polygons.

\section{Lower bound -- higher dimensions}
\label{lowerdd}
Tak{\'a}cs' result was extended by Tak{\'a}cs and Pataki to Euclidean spaces of dimension higher than two:

\begin{proposition}
 (\cite{TakacsPataki})
 \label{Takacs2}
 \begin{equation}
 \mathrm{VCdim}(G(d,N)) \geq dN + 2,
 \end{equation}
 for $d \geq 2$, $N \geq 2$.
\end{proposition}
\begin{proof}
 The geometrical arguments are less simple than for the two-dimensional case; we refer readers to \cite{TakacsPataki} for the proof.
\end{proof}
\begin{remark}
 Tak{\'a}cs and Pataki also offer slightly stronger lower bounds for the special cases $d=3$ and $d=4$: $h(G(3,N)) \geq 3N + 3$ and $h(G(4,N)) \geq 4N + 5$; but these require a higher minimum value for $N$ than Lemma \ref{Takacs2} does.
\end{remark}

As in the two-dimensional case, we can deduce a lower bound for the VC dimension of the NN classifier:
\begin{corollary}
  \label{LowerBoundProp}
  \begin{equation}
  \mathrm{VCdim}(\mathrm{1NN}(d,m)) \geq dm + 2 - d,
  \end{equation}
  for $d \geq 2$, $m \geq 3$.
\end{corollary}
\begin{remark}
The relation does not hold for $d=2$, $m=2$: an NN-rule classifier with two prototypes in the plane has for its decision boundary a line, so cannot shatter four points: linear classifiers are, famously, unable to solve the XOR problem.
\end{remark}

\section{Upper bounds}
\label{upper}

The \emph{shatter coefficient} of a family of sets $B$ (for our purposes, $B$ is a set of subsets of $\mathbb{R}^d$) is a number closely related to VC dimension.  It is called the ``growth function'' by some authors, but other authors give that term a different definition.  The $n$th shatter coefficient of $B$, denoted $S(B,n)$, is the maximum number of different subsets of $n$ points that can be formed by intersection of the $n$ points with elements of $B$: that is, the number of subsets that can be ``picked out'' using elements of $B$.  (The maximum is taken over all sets of $n$ points.)

The VC dimension of $B$ is then the largest $n$ such that $S(B,n) = 2^n$.  This is an expression of the concept of shattering for families of sets $B$ rather than classifiers: if the set of all subsets of the $n$ points which can be formed by intersection of the $n$ points with elements of $B$ is all $2^n$ possible subsets of the $n$ points, then the $n$ points are shattered by $B$.

The VC dimension of a family of classifiers, as we defined it in section \ref{introduction}, is equal to the VC dimension of the associated family of decision regions, as just defined, and the shatter coefficient of a family of classifiers is defined equal to the shatter coefficient of the associated family of decision regions see \citep[see][ch.12]{Devroye96}.

Devroye et al. give upper bounds for the shatter coefficients of the class $C(d,m)$:
\begin{lemma}
 (Devroye, Gy{\"o}rfi, and Lugosi) For $m \geq 3$,
 \label{DevroyeLemma}
 \begin{align}
  S(C,n) &\leq 2^m n^{9(m-2)} &\text{for } d=2\\
  S(C,n) &\leq 2^m n^{(d+1)m(m-1)/2} &\text{for } d \geq 3.
 \end{align}
\end{lemma}
\begin{proof}
 See \cite[p. 312]{Devroye96}
\end{proof}
As before, $d$ is the dimension of the space and $m$ is the number of prototypes used by the classifier.

These upper bounds are based simply on the observation that for a reference set with $m$ prototypes there are at most $m(m-1)/2$ Voronoi cell boundaries, so the number of points which can be shattered by the set of Voronoi diagrams with $m$ centres is bounded above by the number of points which can be shattered by $m(m-1)/2$ hyperplanes.  The stronger result for $d=2$ comes from a restriction on the number of edges of a planar graph, applied to the Delaunay triangulation which is the dual of the Voronoi diagram.

These bounds imply the following result for the VC dimension of the NN classifier:

\begin{proposition}
\label{Wbound}
For $m \geq 3$,
\begin{equation}
 \mathrm{VCdim}(\mathrm{1NN}(d,m)) \leq -\frac{q}{\log{2}} W_{-1} \left(-\frac{\log 2}{q} 2^{-\frac{m}{q}} \right),
 \end{equation}
where
\begin{align}
  q &= 9(m-2) &\text{for } d = 2\\
  q &= (d+1)m(m-1)/2 &\text{for } d \geq 3.
\end{align}
$W$ is the Lambert W function; the $W_{-1}$ branch is the relevant branch. The logarithms are natural logarithms.
\end{proposition}

\begin{proof}
 See appendix \ref{LambertResult}.
\end{proof}

We can obtain from this a looser but more easily interpretable upper bound by using a recent result which gives a lower bound for the $W_{-1}$ branch of the Lambert function:
\begin{lemma}
(Chatzigeorgiou)
  \begin{equation}
    W_{-1}(-e^{-u-1}) > -1 - \sqrt{2u} - u,
  \end{equation}
  for $u>0$.
\end{lemma}
\begin{proof}
  See \cite{Chatz13}.
\end{proof}

 Using this result with proposition \ref{Wbound} gives the following looser bound on the VC dimension:
\begin{corollary}
\label{LooseCor}
 Let $q' = q/\log{2}$, where $q$ is defined as in proposition \ref{Wbound}. Then for $m \geq 3$,
\begin{equation}
 \label{looseupper}
 \mathrm{VCdim}(\mathrm{1NN}(d,m)) < q' \left( \sqrt{2\left(\frac{m}{q'} + \log{q'} -1\right)} + \frac{m}{q'} + \log{q'} \right).
\end{equation}

\end{corollary}
This upper bound enables us to bound the the asymptotic rate of growth of $C$.  Now, $q'$ grows monotonically with $m$ and with $d$, and grows at least as fast as $O(m)$ for increasing $m$.  So whether considering growth with increasing $m$ or growth with increasing $d$, the fastest-growing term within the brackets of (\ref{looseupper}) is $\log{q'}$.  So we have
\begin{equation}
\label{simq}
 \mathrm{VCdim}(\mathrm{1NN}(d,m)) \lesssim q' \log{q'},
\end{equation}
for large $m$ or $d$.

Table~\ref{tab:results} summarises the VC dimension bounds derived in this study.

\section{Discussion of asymptotic behaviour}
\label{discussion}

Considering first the rate of growth of the VC dimension with $m$ for fixed $d$, equation \ref{simq} implies
\begin{equation}
\label{Om2}
\mathrm{VCdim}(\mathrm{1NN}(D,m)) \in O(m^2 \log(m^2)) = O(m^2 \log(m)),
\end{equation}
for fixed $d=D>2$, with the better result
\begin{equation}
\mathrm{VCdim}(\mathrm{1NN}(2,m)) \in O(m \log(m)),
\end{equation}
for $d = 2$. Figure~\ref{fig:VCBoundsFigure} illustrates the upper bounds as a function of $m$, for the cases $d = 2$ and $d = 3$.

\begin{figure}[htb]
\centering
\includegraphics[width=0.65\linewidth]{VCBoundsFigure}
\caption{Upper bounds for the VC dimension of NN-rule classifiers.  ``Accurate'' upper bounds are the best we have obtained, given in Proposition~\ref{Wbound}. ``Approximate'' upper bounds are the looser bounds given in Corollary~\ref{looseupper}.}
\label{fig:VCBoundsFigure}
\end{figure}

It is interesting to compare the log-linear growth in $m^2$ given by equation (\ref{Om2}) with recent results for neural networks given by \cite{Harvey}.  Consider a neural network with $d$ input neurons (the real coordinates of the feature space), and $N$ neurons in a single hidden layer with binary threshold activation functions. Each neuron in the hidden layer of such a network encodes a hyperplane decision boundary: the neuron will be in one or the other of its binary states depending on which side the input vector lies of a plane normal to the vector of weights of that neuron.  So a network with $N=\frac{1}{2}m(m-1)$ hidden neurons is comparable to 1NN classifiers with $m$ prototypes, in the sense that both build their decision boundaries from parts of $\frac{1}{2}m(m-1)$ hyperplanes.  \cite{Harvey} prove that a piecewise-linear neural network with $W$ parameters has VC dimension with asymptotic growth $O(W \log W)$ in the number of parameters.  For the network just described, the number of parameters is $ \sim Nd = \frac{1}{2}m(m-1)d$, meaning the asymptotic growth of its VC dimension is $O(m^2d \log m^2d)$, just as relation (\ref{simq}) gives as an upper bound for our 1NN classifiers.  That is, the neural network achieves the (asymptotically) highest value it can for its VC dimension, given the number of hyperplane decision surfaces it has to work with.  It is an interesting open question whether the same is true for 1NN classifiers.

Turning now to consider the rate of growth with $d$ for fixed $m=M$, equation (\ref{simq}) implies
\begin{equation}
 \mathrm{VCdim}(\mathrm{1NN}(d,M)) \in O(d \log(d)).
\end{equation}
Thus, the VC dimension grows polynomially in $d$ (asymptotically slower than $d^2$).  This has implications for learnability: for example, polynomial growth of the VC dimension of a class with the dimension of the space is a necessary condition for the class to be properly polynomially learnable~\cite[Theorem 3.1.1.]{Blumer89}.

\section{Conclusions}
\label{conclusions}
\begin{table}
\caption{Summary of the VC dimension bounds. $m$ denotes the number of prototypes in the reference set, and $d$ denotes dimensionality. $W_{-1}$ denotes the $-1$ branch of the Lambert W function.}
\label{tab:results}
\centering

\begin{tabular}{ccc}
Type of bound&Dimensionality&Expression\\
\hline
Lower&$d = 2$&$2m + 1$\\
Lower&$d \geq 2$&$dm + 2 - d$\\
Upper&$d = 2$&$-\frac{9(m-2)}{\log{2}} W_{-1} \left(-\frac{\log 2}{9(m-2)} 2^{-\frac{m}{9(m-2)}} \right)$\\
Upper&$d \geq 2$&$-\frac{(d+1)m(m-1)}{2\log{2}} W_{-1} \left(-\frac{2\log 2}{(d+1)m(m-1)} 2^{-\frac{2m}{(d+1)m(m-1)}} \right)$
\end{tabular}

\end{table}

The VC dimension for the set of all NN-rule classifiers in $d$-dimensional Euclidean space with a reference set of size $m$ grows at least as fast as $dm$ and not faster than $O(m^2 \log{m}$) as $m$ increases.  For the case of two-dimensional Euclidean space, the VC dimension grows not faster than $O(m \log{m}$).

Considering instead growth with $d$, the VC dimension for the set of all NN-rule classifiers in $d$-dimensional Euclidean space with a reference set of size $m$ grows at least as fast as $d(m-1)$ and not faster than $O(d \log{d}$) as $d$ increases.

Precise lower and upper bounds for this VC dimension are given in our Corollary \ref{LowerBoundProp} and Proposition \ref{Wbound} respectively, and summarised in Table~\ref{tab:results}.

The consequence of these bounds that is of interest to practitioners is the implication for the size of sample needed to learn an accurate classifier.  In the Probably Approximately Correct learning model, the \emph{sample complexity} is the number of training examples needed to learn a classifier of given accuracy with given probability.  The sample complexity of a family of classifiers is known to depend linearly on the VC dimension~\cite[Theorem 6.8]{Shais}.  Therefore, the bounds we give above for the asymptotic growth of the VC dimension are also bounds on the asymptotic growth of the size of the training data set needed to learn (with given probability) an accurate NN classifier with reference set of given size.  The lower bound applies only to classification algorithms which produce reference sets of given fixed size (and can produce \emph{any} reference set of that size). The upper bound is significantly more broadly applicable. The upper bound applies to any NN classification algorithm which may not have more than $m$ points in its reference set.

Thus we may conclude: the size of the training set needed to learn an accurate NN-rule classifier with reference set of size $m$ in $d$-dimensional Euclidean space grows not faster than $O(m^2 \log{m}$) as $m$ increases.  For the case of two-dimensional Euclidean space, the size of the training set required grows not faster than $O(m \log{m}$).  Considering instead growth with $d$ for fixed $m$, the size of the training set required grows not faster than $O(d \log{d}$).

The fact that the growth rate of the upper bound is asymptotically faster (with $m$) for $d \geq 3$ than for $d=2$ raises the interesting possibility that there may be something fundamentally different about the behaviour of NN-rule classifiers in 3 dimensions and higher from their behaviour in two-dimensional space.  If future work were to establish an $O(m^2)$ lower bound for $d \geq 3$, this would be confirmed. If, instead, an $O(m \log(m))$ upper bound were established, implying the same behaviour for the 1NN classifier in higher dimensions as in 2 dimensions, this would imply instead an interesting discrepancy between the behaviour of 1NN classifiers and the behaviour of neural networks with access to an equal number of hyperplanes from which to construct their decision boundaries, as discussed in section \ref{discussion}.

\section*{Acknowledgment}
This work was done under project RPG-2015-188 funded by The Leverhulme Trust, UK.  While preparing the paper for publication, IG received support from the European Union's Horizon 2020 research and Innovation programme under grant agreement No 731593.

\bibliographystyle{natbib}


\section{Introduction}
\label{sec:introduction}
Coupled systems of free flow and flow through a porous medium can be found ubiquitously in various kinds of natural and industrial contexts, including
soil water evaporation \cite{vanderborght2017a}, fuel cell water management \cite{gurau2009a}, food processing \cite{verboven2006a}, and evaporative cooling for turbomachinery \cite{dahmen2014a}.
Despite every-growing computational capacities, a discrete numerical modeling these kind of systems, including the porous geometry, is only feasible for small-scale problems. For larger domain
sizes, averaging techniques involving the concept of a representative elementary volume (REV) \cite{bear1972a}, are used in order to yield upscaled models for the description of the porous
domain \cite{whitaker1999a}. These models can then be coupled to the free-flow region either using a single-domain or a two-domain approach. For the former, both the porous medium and the
free flow are described by a single set of equations, as first introduced by Brinkman \cite{brinkman1949a}. The two domains are then discerned by a spatial variation of material parameters.
On the other hand, the two-domain approach decomposes the problem into two disjoint subdomains. The free-flow region is then governed by the Navier-Stokes equations while Darcy's or Forchheimer's law is used in the porous medium subdomain
\cite{ochoa1995a, layton2002, jamet2009a, mosthaf2011a}. In order to maintain thermodynamic consistency, appropriate coupling conditions have to be formulated which enforce the conservation
of mass, momentum and energy across the interface between the two domains \cite{hassanizadeh1989a}.
The aim of this work is to couple two different discretization schemes, and we will focus on the two-domain approach since it provides this flexibility more readily.
We will consider laminar single-phase flow in the following but the extension to compositional multi-phase flow in the porous domain \cite{mosthaf2011a} and the use of turbulence models in the free flow part \cite{fetzer2017a} is possible.

Being an active area of research, various mathematical and numerical models for the
coupling of (Navier-) Stokes and Darcy/Forchheimer have been developed during recent years.
Examples using the same spatial discretization scheme for both domains range from the staggered-grid finite volume method \cite{iliev2004a, harlow1965a},
finite element method \cite{discacciati2009a} or the box-scheme, a vertex-centered finite volume method \cite{mosthaf2011a, huber2000a}. Furthermore,
combinations of colocated and staggered-grid schemes have been developed \cite{rybak2015a, masson2016a, fetzer2017a}. For more discretization approaches
 to the modeling of coupled free flow and porous medium flow, we refer the reader to \cite{discacciati2009a} and references therein.

In this work, the free flow equations are discretized using the staggered-grid method, which forms a stable numerical method for such problems. In turn, no additional stabilization techniques are necessary, in contrast to colocated schemes \cite{versteeg2007a}, for example. In the porous medium, we employ a multi-point-flux-approximation method (MPFA) \cite{aavatsmark2002a}. This method was developed to overcome the shortcomings of the classical two-point-flux approximation.
In particular, the MPFA scheme does not require the grid to be $\mathbf{K}$-orthogonal, meaning that the grid cells need not be in line with the principal directions of the permeability
tensor $\mathbf{K}$. This is especially important for skewed and unstructured grids or in case the principal directions of the permeability tensor are inclined, in the case of
layered porous structures or faults, for example. The MPFA method has been applied previously for solving Brinkman's equation \cite{iliev2014a, brinkman1949a} for a coupled system of free flow and porous medium. The referenced works thus use the same discretization scheme for both subdomains.
The novelty of this work is that we employ a staggered-grid method for the free-flow model (Navier-Stokes), and couple it with a MPFA method for the porous medium model (Darcy).

The paper is organized as follows. Section~\ref{sec:equations} introduces the coupled model by presenting the governing equations in the free flow and porous medium flow subdomains, respectively, as well as the coupling conditions at the interface.
The discretization schemes are introduced in Section~\ref{sec:discretization}, including the newly proposed coupling at the interface between the two subdomains.
Numerical results are presented in Section~\ref{sec:numresults} with the use of two test cases.
Finally, Section~\ref{sec:conclusions} focuses on the conclusions.

\section{Governing Equations}
\label{sec:equations}

In order to present the governing equations for the coupled model, we first introduce the assumptions on the geometry and the notational conventions.
After these preliminary definitions, we continue with the equations governing free flow and porous medium flow, followed by the coupling conditions.

Let $\Omega\subset\mathbb{R}^d$, $d\in \{2, 3\}$, be an open, connected
Lipschitz domain with boundary $\partial\Omega$ and $d$-dimensional measure $\meas{\Omega}$.
Furthermore, let $\Omega^\mathrm{pm}$ and $\Omega^\mathrm{ff}$ be a disjoint partition of $\Omega$ representing the porous medium and free flow subdomains, respectively.
The subdomain boundaries are given by the interface $\Gamma^\mathrm{if} := \partial \Omega^\mathrm{ff} \cap \partial \Omega^\mathrm{pm}$ as well as the remainders $\Gamma^\mathrm{pm} := \partial \Omega^\mathrm{pm} \setminus \Gamma^\mathrm{if}$ and $\Gamma^\mathrm{ff} := \partial \Omega^\mathrm{ff} \setminus \Gamma^\mathrm{if}$. For brevity, the superscripts are often omitted when no ambiguity arises.

The external boundary $\partial \Omega = \Gamma^\mathrm{ff} \cup \Gamma^\mathrm{pm}$ is further decomposed such that $\Gamma^\mathrm{ff} = \Gamma_\mathrm{v}^\mathrm{ff} \cup \Gamma_\mathrm{p}^\mathrm{ff}$ and $\Gamma^\mathrm{pm} = \Gamma_\mathrm{v}^\mathrm{pm} \cup \Gamma_\mathrm{p}^\mathrm{pm}$ disjointly. Here, the subscript denotes whether the velocity or the pressure is prescribed as a boundary condition.
To ensure unique solvability of the resulting system, we assume that $\meas{\Gamma_\mathrm{p}^\mathrm{pm} \cup \Gamma_\mathrm{p}^\mathrm{ff}} > 0$, i.e. that a pressure boundary condition is imposed on a subset of the boundary $\partial \Omega$ with positive measure.

Let $\mathbf{n}$ denote the unit normal vector on $\partial \Omega$ oriented outward with respect to $\Omega$. We abuse notation and let $\mathbf{n}$ moreover denote the unit normal vector on $\Gamma^\mathrm{if}$ oriented outward with respect to $\Omega^\mathrm{ff}$.

\subsection{Free Flow}
\label{sub:free_flow}

In our model, the Navier-Stokes equations govern the free flow in subdomain $\Omega^\mathrm{ff}$. These equations are given by:
\begin{subequations}
  \label{eq:navierstokes}
  \begin{align}
      \frac{\partial \varrho}{\partial t} + \nabla \cdot (\varrho \mathbf{v}) &= q, \\
  \frac{\partial (\varrho \mathbf{v})}{\partial t} + \nabla \cdot \left(\varrho \mathbf{v} \mathbf{v}^{\mathrm{T}}
      - \mu (\nabla \mathbf{v} + (\nabla \mathbf{v})^{\mathrm{T}})
      + p \mathbf{I} \right)
      &= \varrho \textbf{g}, &&\mathrm{in} \, \Omega^\mathrm{ff}, \\
       \mathbf{v} &= \mathbf{v}_\Gamma, &&\mathrm{on} \, \Gamma^\mathrm{ff}_\mathrm{v}, \\
       p &= p_\Gamma, &&\mathrm{on} \, \Gamma^\mathrm{ff}_\mathrm{p}.
   \end{align}
\end{subequations}

The unknown variables are the velocity $\mathbf{v}$ and the pressure $p$. Here, $\varrho$ and $\mu$ denote the potentially pressure-dependent density and viscosity, respectively, while $q$ is a source (or sink) term, and $\mathbf{g}$ describes the influence of gravity. $\mathbf{I}$ is the identity tensor in $\mathbb{R}^{d \times d}$. The boundary conditions are given by known quantities $\mathbf{v}_\Gamma$ and $p_\Gamma$, representing the velocity or pressure at the corresponding boundary.

\subsection{Porous Medium Flow}
\label{sub:porous_medium_flow}

The equations governing single-phase flow in the porous-medium $\Omega^\mathrm{pm}$ are given by
\begin{subequations}
    \label{eq:darcy}
  \begin{align}
      \frac{\partial \varrho}{\partial t} + \nabla \cdot \left( \varrho \mathbf{v} \right) &= q, \\
      \mathbf{v} + \frac{1}{\mu} \mathbf{K} \left( \nabla p - \varrho \mathbf{g}\right) &= 0, &&\mathrm{in} \, \Omega^\mathrm{pm}, \label{eq: darcys law}\\
   \mathbf{v} \cdot \mathbf{n} &= v_\Gamma, &&\mathrm{on} \, \Gamma^\mathrm{pm}_\mathrm{v}, \\
       p &= p_\Gamma, &&\mathrm{on} \, \Gamma^\mathrm{pm}_\mathrm{p}.
  \end{align}
\end{subequations}
Equation \eqref{eq: darcys law} states that the momentum balance in the porous medium is given by Darcy's law, i.e. the Darcy velocity is calculated as $\mathbf{v} = -\frac{1}{\mu} \mathbf{K} \left( \nabla p - \varrho \mathbf{g}\right)$, with $\mathbf{K}$ being the permeability tensor. Similar to  the free-flow equations \eqref{eq:navierstokes}, $q$ denotes a source or sink term. Finally $v_\Gamma$ and $p_\Gamma$ are known quantities representing the normal flux and pressure on the corresponding boundaries, respectively.

\subsection{Coupling Conditions}
\label{sub:coupling}

In order to derive a thermodynamically consistent formulation of the coupled problem,
 conservation of mass and momentum has to be guaranteed at the interface between the porous
 medium and the free-flow domain. We therefore impose the following interface conditions:

\begin{subequations}\label{eq:interfaceConditions}
 \begin{align}
        \label{eq:interfaceConditionsFlux}
    \mathbf{v}^\mathrm{ff} \cdot \mathbf{n} &=  - \mathbf{v}^\mathrm{pm} \cdot \mathbf{n} , \\
    \label{eq:interfaceConditionsMomentum}
    \left(\varrho \mathbf{v} \mathbf{v}^\mathrm{T} - \mu (\nabla \mathbf{v} + (\nabla \mathbf{v})^\mathrm{T}) + p \mathbf{I} \right)^\mathrm{ff} \mathbf{n} &=
      p^\mathrm{pm} \mathbf{n}, \\
     \label{eq:interfaceConditionsBJS}
     \left( - \frac{\sqrt{\mathbf{t} \cdot \mathbf{K} \mathbf{t}}}{\alpha_{\mathrm{BF}}} (\nabla \mathbf{v}) \mathbf{n} - \mathbf{v} \right)^\mathrm{ff} \cdot \mathbf{t} &=  0, && \text{on} \, \Gamma^\mathrm{if}.
 \end{align}
\end{subequations}

The momentum transfer normal to the interface is given by \eqref{eq:interfaceConditionsMomentum} \cite{layton2002}.
Condition \eqref{eq:interfaceConditionsBJS} is the commonly used Beavers-Joseph-Saffman slip condition \cite{beavers1967a, saffman1971a}. Here, $\mathbf{t}$ denotes any unit vector from the tangent bundle of $\Gamma^\mathrm{if}$ and $\alpha_{\mathrm{BF}}$ is a parameter to be determined experimentally. We remark that this condition is technically a boundary condition for the free flow, not a coupling condition between the two flow regimes. Furthermore, it has been developed for free flow strictly parallel to the
interface and might lose its validity for other flow configurations.

\section{Discretization}
\label{sec:discretization}

This section is devoted to giving an outline of the numerical schemes used in the individual subdomains
 and the incorporation of the interface conditions \eqref{eq:interfaceConditions}. However, we first introduce some notational conventions concerning the partition of $\Omega$ in the following definition.

\begin{definition}[Grid discretization]
\label{def:griddisc}
The tuple \mbox{$\mathcal{D} := (\mathcal{T},\mathcal{E},\mathcal{P},\mathcal{V})$} denotes the grid discretization, in which
  \begin{enumerate}[label=(\roman*)]
  \item $\mathcal{T}$ is the set of control volumes (cells) such that
    \mbox{$\overline{\Omega}= \cup_{K \in \mathcal{T}} \overline{K}$}.
    For each cell $K\in\mathcal{T}$, $\meas{K}>0$ denotes the cell volume.
 \item
    $\mathcal{E}$ is the set of faces such that each face $\sigma$ is a $(d-1)$-dimensional hyperplane
    with measure
    \mbox{$\meas{\sigma}>0$}.
    For each cell $K \in \mathcal{T}$, $\mathcal{E}_K$ is the subset of
    $\mathcal{E}$ such that \mbox{$\partial K  = \cup_{\sigma \in\mathcal{E}_K}{\sigma}$}. Furthermore, $\mathbf{x}_\sigma$ denotes
    the face evaluation points and $\mathbf{n}_{K,\sigma}$ the unit vector that is normal to
    $\sigma$ and outward to $K$.

  \item
    $\mathcal{P} := \lbrace \mathbf{x}_K\rbrace_{K \in \mathcal{T}}$ is the set of \emph{cell centers} (not required to be the barycenters) such that
    $\mathbf{x}_K\in K$ and $K$ is star-shaped with respect to $\mathbf{x}_K$.
    For all $K\in\mathcal{T}$ and $\sigma\in\mathcal{E}_K$, let $d_{K,\sigma}$ denote the Euclidean distance between $\mathbf{x}_K$ and $\sigma$.

  \item
    $\mathcal{V}$ is the set of vertices of the grid, corresponding to the corners of the cells.
  \end{enumerate}
\end{definition}

For ease of exposition, we assume $d = 2$, however, this is not a limitation and the model can be readily extended for three dimensions.

\subsection{Staggered grid scheme}
\label{sub:staggered}

A staggered-grid finite volume scheme, also known as MAC scheme \cite{harlow1965a}, is used in the free-flow subdomain.
Here, scalar quantities including pressure and density are stored on the cell centers $\mathcal{P}$ while the velocity degrees of freedom are located on the primary control volumes' faces $\mathcal{E}$.
The resulting scheme is stable, hence oscillation-free solutions are guaranteed without the need for additional stabilization techniques. This is in contrast with colocated schemes, in which all unknowns are defined at the same location \cite{versteeg2007a}. $\mathcal{D}$ is a uniform Cartesian grid with mesh size $h$ in both directions.

We start with the momentum balance equations. For each face in $\mathcal{E}$, we construct a secondary control volume $K^*$ with boundary $\partial K^*$, as depicted in Figure~\ref{pc:staggeredgrid}. On each $K^*$, we integrate the momentum balance equation and apply Gauss divergence theorem to obtain
\begin{equation}
  \begin{aligned}
    \int_{K^*} \frac{\partial (\varrho \mathbf{v})}{\partial t} \, \mathrm{d} x + \int_{\partial K^*} (\varrho \mathbf{v} \mathbf{v}^{\mathrm{T}}) \cdot \mathbf{n} \, \mathrm{d} \Gamma
    &- \int_{\partial K^*} (\mu (\nabla \mathbf{v} + (\nabla \mathbf{v})^{\mathrm{T}})) \cdot \mathbf{n} \, \mathrm{d} \Gamma
    + \int_{\partial K^*} p \mathbf{n} \, \mathrm{d} \Gamma
    = \int_{K^*} \varrho \mathbf{g} \, \mathrm{d} x.
  \end{aligned}
\end{equation}

The first and second components of this vector equation are considered separately. Due to the different locations at which the degrees of freedom are defined, we require interpolation operators in order to continue.
Let us introduce the average ($\{ \cdot \}$) and jump quantities ($\llbracket \cdot \rrbracket$) on cell centers ($\mathcal{P}$) and vertices ($\mathcal{V}$) of the primal grid for $\mathbf{v} = [v_x, v_y]^{\text{T}}$ (see Figure \ref{pc:staggeredgrid}):
\begin{equation}
\begin{aligned}
  \{ \mathbf{v} \}|_{\mathcal{P}} &= \frac{1}{2}
  \begin{bmatrix}
    v_x^E + v_x^W \\ v_y^N + v_y^S
  \end{bmatrix}
  & \quad
  \{ \mathbf{v} \}|_{\mathcal{V}} &= \frac{1}{2}
  \begin{bmatrix}
    v_x^N + v_x^S \\ v_y^E + v_y^W
  \end{bmatrix} \\
  \llbracket \mathbf{v} \rrbracket|_{\mathcal{P}} &= \frac{2}{h}
  \begin{bmatrix}
    v_x^E - v_x^W \\ v_y^N - v_y^S
  \end{bmatrix}
  & \quad
  \llbracket \mathbf{v} \rrbracket|_{\mathcal{V}} &= \frac{1}{h}
  \begin{bmatrix}
    v_x^N - v_x^S + v_y^E - v_y^W \\
    v_x^N - v_x^S + v_y^E - v_y^W
  \end{bmatrix}
\end{aligned}
\end{equation}
Here, the superscript $\{E, N, W, S\}$ refers to the closest degree of freedom East, North, West, or South of the evaluation point, see Figure \ref{pc:staggeredgrid}.
\begin{figure} [ht!]
  \centering
  \includegraphics[width=1.0\linewidth,keepaspectratio]{pics/staggeredgrid.pdf}
  \caption{Grid and notations used for the staggered-grid discretization. $K_x^*, K_y^*$ denote the dual cells. The picture on the left illustrates the situation where $\mathbf{x}_{\sigma^*}$ coincides with a cell center, whereas the picture on the right shows the case where the center $\mathbf{x}_{\sigma^*}$ of a dual face $\sigma^*$ coincides with a vertex of the primary grid.}
  \label{pc:staggeredgrid}
\end{figure}
In the following, the superscript ``up'' denotes the upwind quantity relative to the velocity $\mathbf{v}$.
Moreover, we introduce $\mu^\mathrm{avg}$ such that in all cell centers $\mathcal{P}$, it denotes the corresponding viscosity whereas in each vertex $\mathcal{V}$, it is the viscosity averaged over the adjacent cells.

With the operators defined, we discretize the momentum balance equation for each secondary control volume $K_i^*$, in which the subscript $i \in \{ x, y \}$ denotes whether the control volume surrounds a vertical or horizontal face, respectively. The discretized equation for component $i$ is then given by
\begin{equation}
  \begin{aligned}
    \int_{K_i^*} \frac{\partial (\varrho v_i)}{\partial t} \, \mathrm{d} x
    + \int_{\partial K_i^*} (\{ \mathbf{v} \} (\varrho v_i)^\mathrm{up}) \cdot \mathbf{n} \, \mathrm{d} \Gamma
    - \int_{\partial K_i^*} \mu^\mathrm{avg} \llbracket \mathbf{v} \rrbracket \cdot \mathbf{n} \, \mathrm{d} \Gamma
    &+ \int_{\partial K_i^*} p n_i \, \mathrm{d} \Gamma
    = \int_{K^*} \varrho g_i \, \mathrm{d} x.
  \end{aligned}
\end{equation}

We emphasize that the boundary integrals are computed numerically using the following rules for a scalar-valued function $f$ and a vector-valued function $\mathbf{f}$:
\begin{subequations}
  \begin{align}
    \int_{\partial K^*} f \, \mathrm{d} \Gamma
    &= \sum_{\sigma^* \in \mathcal{E}_{K^*}} \meas{\sigma^*} f(\mathbf{x}_{\sigma^*})
    = h \sum_{\sigma^* \in \mathcal{E}_{K^*}} f(\mathbf{x}_{\sigma^*}), \\
    \int_{\partial K^*} \mathbf{f} \cdot \mathbf{n} \, \mathrm{d} \Gamma
    &= \sum_{\sigma^* \in \mathcal{E}_{K^*}}   \meas{\sigma^*}\mathbf{f}(\mathbf{x}_{\sigma^*}) \cdot \mathbf{n}_{K^*, \sigma^*}
    = h \sum_{\sigma^* \in \mathcal{E}_{K^*}} \mathbf{f}(\mathbf{x}_{\sigma^*}) \cdot \mathbf{n}_{K^*, \sigma^*}.
  \end{align}
\end{subequations}
By definition, each $\mathbf{x}_{\sigma^*}$ will either be a cell center ($\mathcal{P}$) or a vertex ($\mathcal{V}$) of the grid.

Finally, the mass balance is evaluated on each cell of the grid, i.e. we compute for each $K \in \mathcal{T}$:
\begin{equation}
  \int_{K} \frac{\partial \varrho}{\partial t} \, \mathrm{d} x
  + \int_{\partial K} \varrho^\mathrm{up} \mathbf{v} \cdot \mathbf{n} \, \mathrm{d} \Gamma
  =
  \int_{K} q \, \mathrm{d} x.
\end{equation}

The above equations fully define the staggered-grid discretization scheme for all internal faces of the grid $\mathcal{D}$. For incorporation of boundary conditions, we refer the reader to \cite{versteeg2007a}.

\subsection{Cell-centered finite volume scheme}
\label{sub:ccfvm}

A cell-centered finite volume scheme is employed in the Darcy subdomain, i.e. the
 grid elements are used as control-volumes and the degrees of freedom are associated
 with the cell centers. Typically,
 the finite volume formulation is obtained by integrating the first equation of
 \eqref{eq:darcy} over a control volume $K \in \mathcal{T}$ and by applying the Gauss divergence theorem:
 \begin{equation}
 \int_K \frac{\partial \varrho}{\partial t} \, \mathrm{d}x +   \sum_{\sigma \in \mathcal{E}_K }\int_{\sigma} \frac{\varrho^\mathrm{up}}{\mu^\mathrm{up}}\left( - \mathbf{K} \left( \nabla p - \varrho \mathbf{g}\right) \right) \cdot \mathbf{n} \, \mathrm{d} \Gamma = \int_K q \, \mathrm{d}x.
     \label{eq:darcyIntegrated}
 \end{equation}

Replacing, the exact fluxes by an  approximation, i.e. \mbox{$F_{K, \sigma} \approx \int_{\sigma} \left( - \mathbf{K}_K \left( \nabla p - \varrho \mathbf{g}\right) \right) \cdot \mathbf{n} \, \mathrm{d} \Gamma$}
 (here $\mathbf{K}_K$ is the value of $\mathbf{K}$ associated with control volume $K$), yields
\begin{equation}
 \int_K \frac{\partial \varrho}{\partial t} \, \mathrm{d}x + \sum_{\sigma \in \mathcal{E}_K} \frac{\varrho^\mathrm{up}}{\mu^\mathrm{up}} F_{K, \sigma} = Q_K, \quad \forall \, {K \in \mathcal{T}},
  \label{eq:darcyCCdiscrete}
\end{equation}
where $F_{K, \sigma}$ is the discrete flux through face $\sigma$ flowing out of cell
 $K$, $Q_K := \int_K q \, \mathrm{d}x$ is the integrated source/sink term, and $(\cdot)^\mathrm{up}$ denotes upwinding with respect to the sign of the flux $F_{K,\sigma}$.

Finite volume schemes primarily differ in the approximation of the term
$(\mathbf{\mathbf{K}}_K \nabla p) \cdot \mathbf{n}$
(i.e. the choice of the fluxes $F_{K, \sigma}$). The widely used linear two-point flux approximation (TPFA), for example, is a simple but robust scheme.
However, it is well-known that it is inconsistent on grids that are not $\mathbf{K}$-orthogonal
(see e.g.\ \cite{edwards1998finite}). In this work we consider anisotropic permeability tensors in
the porous medium and $\mathbf{K}$-orthogonality of the grid can thus not be guaranteed.
Therefore, we employ a multi-point flux approximation (MPFA) scheme for the formulation of the discrete
fluxes, which has been presented in \cite{aavatsmark2002a}. This particular scheme is termed MPFA-O and
 is only one among many methods that fall into the family of MPFA schemes (\cite{aavatsmark2008compact, edwards2010mpfafps}).
 Please note that we will omit the suffix ``-O'' throughout this document wherever it would affect the readability.

 \begin{figure} [ht!]
  \centering
  \includegraphics[width=0.8\linewidth,keepaspectratio]{pics/interactionvolume_interior.pdf}
  \caption{Interaction region for the MPFA-O method. The parameter $\xi$, $0 \le \xi < 1$ is used to define the location of the intermediate
            face pressure unknowns $p_{\sigma_i^v}$. Here, the situation for $\xi = 0.5$ is illustrated.}
  \label{pc:interactionRegion_interior}
\end{figure}

For the computation of the fluxes, a dual grid is created
 by connecting the barycenters of the cells with the barycenters of the faces ($d=2$)
 or the barycenters of the faces and edges ($d=3$). This divides each cell into sub-control
 volumes $K^v$ with $K \in \mathcal{T}$ and $v \in \mathcal{V}$. Analogously, each face is sub-divided into sub-control volume faces $\sigma^v$,
 see Figure \ref{pc:interactionRegion_interior}. Expressions for
the face fluxes $F_{K, \sigma^v}$ are obtained by introducing the face
pressures $p_{\sigma^v}$. The location of these face pressures along the sub-control volume faces $\sigma^v$
is parameterized by $\xi$, $0 \le \xi < 1$, and is the center of the original face $\sigma$ for
$\xi = 0$ and would be the position of the vertex $v$ for $\xi = 1$. These face pressures are then eliminated by enforcing
the continuity of fluxes across each sub-control volume face.
I.e., for each face $\sigma^v$ between $K^v$ and $L^v$, we  impose:
\begin{subequations}
 \begin{align}
     &F_{K, \sigma^v} + F_{L, \sigma^v} = 0.

     \label{eq:sigmaConditions}
 \end{align}
\end{subequations}

We allow for piecewise constant $\mathbf{K}$
 (denoted as $\mathbf{K}_K$ for each cell $K$) and construct discrete gradients
 $\nabla_\mathcal{D}^{K^v} p$, per sub-control volume $K^v$, depending on its two embedded
 sub-control volume faces.
 Let us consider $K^v$ in Figure~\ref{pc:interactionRegion_interior} with faces $\sigma_1^v$ and $\sigma_3^v$. Here, the discrete
 gradients are constructed to be consistent such that the following holds for $i \in \{1, 3\}$:
\begin{equation}
  \nabla_\mathcal{D}^{K^v} p \cdot (\mathbf{x}_{\sigma^v_i}- \mathbf{x}_{K}) =
  p_{\sigma^v_i} - p_K.
\label{eq:piecewiselin}
\end{equation}

Thus, a discrete gradient (for sub-control volume $K^v$) that fulfills the
two conditions \eqref{eq:piecewiselin} is defined by
\begin{equation}
  \nabla_\mathcal{D}^{K^v} p  = \mathbb{D}^{-T}_{K^v}
    \begin{bmatrix}
      p_{\sigma^v_1} - p_K \\
      p_{\sigma^v_3} - p_K
    \end{bmatrix}, \qquad \text{ with }\; \mathbb{D}_{K^v} :=
    \begin{bmatrix}
      \mathbf{x}_{\sigma^v_1}- \mathbf{x}_K & \mathbf{x}_{\sigma^v_3} - \mathbf{x}_K
    \end{bmatrix}.
    \label{eq:MPFAGradientRecons}
\end{equation}

This enables us to write the discrete flux across $\sigma^v$ between $K^v$ and $L^v$ as follows:
\begin{equation}
    F_{K, \sigma^v} := - |\sigma^v| \mathbf{n}_{K, \sigma^v}^T \mathbf{\mathbf{K}}_K \nabla_\mathcal{D}^{K^v} p + \gamma_{K, \sigma^v},
    \label{eq:discreteFlux}
\end{equation}
where we introduced $\gamma_{K, \sigma^v} = \rho_{\sigma^v} \meas{\sigma^v} \mathbf{n}_{K, \sigma^v}^T \mathbf{\mathbf{K}}_K \mathbf{g}$,
with $\rho_{\sigma^v} = \frac{\rho_K+\rho_L}{2}$, to incorporate the effect of gravity.

To deduce a cell-centered scheme, the face pressures $p_{\sigma^v_i}$
 are eliminated. This is done by enforcing flux continuity \eqref{eq:sigmaConditions}
 within each interaction
 volume and by solving a local system of equations.

We rewrite these conditions in matrix form, and introduce the sans serif font to denote the corresponding matrices and vectors. All local
 face pressures of an interaction region are collected in the vector $\mathsf{p_{\sigma}}$, cell pressures in the vector
 $\mathsf{p_K}$, and all terms related to gravity in the vector $\mathsf{g}$. Flux continuity then allows us to rewrite the face pressures in terms of the cell pressures:
\begin{equation}
  \mathsf{A} \, \mathsf{p_{\sigma}} = \mathsf{B} \, \mathsf{p_K} + \Delta \mathsf{g}.
  \label{eq:ivLocalSystem}
\end{equation}

Here, the $\Delta$ represents the difference in contributions due to gravity over each face.
Let $\mathsf{f}$ denote the vector of all fluxes across the sub-control volume faces of the interaction
 region. These can be expressed in matrix form using equation \eqref{eq:discreteFlux}:
\begin{equation}
  \mathsf{f} = \mathsf{C} \, \mathsf{p_{\sigma}} + \mathsf{D} \, \mathsf{p_K} + \mathsf{g}.
  \label{eq:ivLocalFluxes}
\end{equation}

With these introduced matrices the final expressions for the local
 sub-control volume face fluxes, related to the interaction region, read:
\begin{equation}
  \mathsf{f} = \underbrace{\left( \mathsf{C} \mathsf{A}^{-1} \mathsf{B}
         + \mathsf{D} \right)}_{=: \mathsf{T}} \, \mathsf{p_K}
         + \mathsf{C} \mathsf{A}^{-1} \left( \Delta \mathsf{g} \right) + \mathsf{g}.
  \label{eq:ivLocalFluxesTrans}
\end{equation}
The entries of the matrix $\mathsf{T}$ are often referred to as the transmissibilities.

\subsection{Coupling}
\label{sub:discrete_coupling}

In this section, we consider the realization of the coupling conditions
\eqref{eq:interfaceConditions}. As depicted in Figure~\ref{pc:interactionRegion_interface}, the
grids are chosen to be non-matching at the interface such that each sub-control volume face
coincides with a face $\sigma \in \mathcal{E}^\mathrm{ff}$ of the free-flow domain. In turn, a natural coupling
arises between the staggered grid discretization and the MPFA method, due to the coinciding degrees
of freedom (for $\xi = 0.5$, see Figure \ref{pc:interactionRegion_interface}).
 We emphasize that each cell in the porous-medium subdomain has two neighboring cells in
the free-flow subdomain (two-dimensional setup).

We start with the flux continuity condition \eqref{eq:interfaceConditionsFlux}. As depicted in
Figure~\ref{pc:interactionRegion_interface}, let us consider a sub-control volume $K^v$ located at the
interface such that $\sigma_3^v$ is located on $\Gamma^\mathrm{if}$. We then let the velocity from the free-flow
 domain determine the flux over $\sigma_3^v$:
\begin{equation}
\begin{aligned}
    F_{K, \sigma_3^v} &= - \mu^{up} \meas{\sigma_3^v} \veli{y, \sigma_3^v}^\mathrm{ff}.

  \label{eq:interactionRegionConditions}
\end{aligned}
\end{equation}

Collecting the right-hand side of \eqref{eq:interactionRegionConditions} in the matrix-vector product
 $\mathsf{Wv_\mathrm{ff}}$, we can rewrite \eqref{eq:ivLocalSystem} for interaction regions
 that are located at the interface to the free-flow domain:
\begin{equation}
  \mathsf{A} \, \mathsf{p_{\sigma}} = \mathsf{B} \, \mathsf{p_K} + \mathsf{Wv_\mathrm{ff}} + \Delta \mathsf{g}.
  \label{eq:ivLocalSystemInterface}
\end{equation}

In the situation shown in Figure~\ref{pc:interactionRegion_interface}, the face pressures
$p^{\mathrm{pm}}_{\sigma^v_1}, p^{\mathrm{pm}}_{\sigma^v_2}, p^{\mathrm{pm}}_{\sigma^v_3}$ are thus dependent on the
primary unknowns $p^\mathrm{pm}_K, p^\mathrm{pm}_L$ of the porous-medium domain and the face velocities
$v^\mathrm{ff}_{\sigma^v_2},v^\mathrm{ff}_{\sigma^v_3}$ of the free-flow domain (i.e. $p^{\mathrm{pm}}_{\sigma^v_i} =
p^{\mathrm{pm}}_{\sigma^v_i}(p^\mathrm{pm}_K, p^\mathrm{pm}_L,v^\mathrm{ff}_{\sigma^v_2},v^\mathrm{ff}_{\sigma^v_3})$). Insertion of
\eqref{eq:ivLocalSystemInterface} in \eqref{eq:ivLocalFluxes} leads to the expression
\begin{equation}
  \mathsf{f} = \mathsf{T} \, \mathsf{p_K}
         + \mathsf{C} \mathsf{A}^{-1} \left( \mathsf{Wv_\mathrm{ff}} + \Delta \mathsf{g} \right) + \mathsf{g}
  \label{eq:ivLocalFluxesInterface}
\end{equation}
for the fluxes across sub-control volume faces within interaction regions located at the interface to the
 free-flow domain.

\begin{figure} [ht!]
  \centering
  \includegraphics[width=1.0\linewidth,keepaspectratio]{pics/interactionvolume.pdf}
  \caption{Interaction region for the MPFA-O method at the interface to the free-flow subdomain. The graphic illustrates the
            non-matching grids at the interface and the choice of $\xi=0.5$ for the MPFA scheme such that the degrees of freedom
            for the face velocities in the free-flow domain coincide with the intermediate face pressure unknowns introduced on the interface.}
  \label{pc:interactionRegion_interface}
\end{figure}

The remaining coupling conditions are imposed as followed. The momentum balance
\eqref{eq:interfaceConditionsMomentum} is enforced using the reconstructed face pressures from
\eqref{eq:ivLocalSystemInterface}. On the other hand, the Beavers-Joseph-Saffman condition
\eqref{eq:interfaceConditionsBJS} is technically a boundary condition for the free-flow problem, as
previously noted in Section~\ref{sec:equations}, and is implemented accordingly

Finally, we remark that in the case of compressible fluids, the density and viscosity are
pressure-dependent. The upwind terms $\mu^{\mathrm{pm},\mathrm{up}}, \varrho^{\mathrm{pm},\mathrm{up}}, \varrho^{\mathrm{ff},\mathrm{up}}$ are
then evaluated using the cell-pressure unknowns, i.e. for a face $\sigma$ between porous-medium
cell $K$ and free-flow cell $L$, these terms therefore depend on the pressures $p^\mathrm{pm}_K$ and
$p^\mathrm{ff}_L$.

\section{Numerical results}
\label{sec:numresults}
All simulations are performed using the open-source simulator {Du\-Mu$^\text{x}$ } \cite{dumux}, which comes in the form of an additional DUNE module \citep{blatt.ea:2016}.
We employ a monolithic approach, where both sub-problems are assembled into one system of equations and use an implicit Euler method for the time discretization.
Newton's method is used to solve the non-linearities involved in the systems of equations.
For all test cases, the compressible fluid ``air'' (see the {Du\-Mu$^\text{x}$ } documentation \cite{dumux}) is used. We consider two-dimensional setups, however, the implementation is also able
to handle three-dimensional domains.

\subsection{Test case 1}

The first test case is similar to the one that has been presented in \cite{iliev2004a}.
However, the authors of \cite{iliev2004a} used a Navier-Stokes-Brinkman-type system for
both domains and this system is discretized by using a staggered-grid scheme in both domains.
For anisotropic permeability tensors, this requires the interpolation of all velocity components at grid faces.

In this work, we use a different approach, where a staggered-grid scheme is used to discretize
the free-flow system and a cell-centered finite volume scheme to discretize the porous-medium system.
Thus, no additional velocity degrees of freedom are needed for the porous-medium domain.
However, for anisotropic permeability tensors or unstructured grids this requires more sophisticated cell-centered finite volume schemes,
as for example the MPFA scheme that has been presented in Section \ref{sub:ccfvm}.

Air is flowing through a two-dimensional channel which is partially blocked by a rectangular porous medium as shown in Figure \ref{pc:settingCaseOne}. The first test case involves small Reynolds
numbers ($Re \ll 1$) with respect to the average velocity in the narrow section in the channel above the porous medium. In this case, only the stationary solution (which is reached after a few time
steps, starting from a resting fluid) is investigated. The top and bottom of the domain are considered as rigid, impermeable walls with $\mathbf{v} = 0$ (including the wall part below the porous box). Flow is driven by a pressure difference between the left and the right
boundary which is set to $\Delta p = \unit[10^{-6}]{Pa}$.

\begin{figure} [ht!]
	\centering
	\includegraphics[width=0.7\linewidth,keepaspectratio]{pics/setting_testOne.pdf}
	\caption{Setting for the first test case.}
	\label{pc:settingCaseOne}
\end{figure}

The permeability tensor in the porous medium is given as
\begin{equation}
\mathbf{K} =  \mathbf{R}(\alpha)
\begin{pmatrix}
\frac{1}{\beta}k & 0 \\
0 & k \\
\end{pmatrix} \mathbf{R}^{-1}(\alpha), \quad
\text{with} \quad \mathbf{R}(\alpha) = \begin{pmatrix}
\cos \alpha & -\sin \alpha \\
\sin \alpha & \cos \alpha \\
\end{pmatrix},
\label{eq:RotatedPermTensor}
\end{equation}
where $\alpha$ is the rotation angle, $\beta \geq 1$ the anisotropy ratio, and $\unit[k=10^{-6}]{m^2}$. In the following, the influence of the anisotropy on the total mass fluxes crossing the boundaries $\Gamma^\mathrm{pm}_\mathrm{in}, \Gamma^\mathrm{pm}_\mathrm{out}, \Gamma^\mathrm{pm}_\mathrm{top}$ is investigated for $\alpha \in \lbrace  -45^\circ, -30^\circ, 0^\circ, 30^\circ, 45^\circ  \rbrace$ and $\beta \in \lbrace  10 ,100  \rbrace$.
In the free-flow domain the same grid is used for the MPFA and the TPFA schemes, whereas in the porous medium the grid used for the MPFA scheme is coarsened such that we obtain a non-matching interface, as shown in Figure \ref{pc:interactionRegion_interface}.

Figure \ref{pc:testOneVelMpfa} and \ref{pc:testOneVelTpfa} show the velocity magnitude and the pressure profile
of the MPFA and TPFA schemes for $\beta = 100$ and $\alpha \in \lbrace -45^\circ, 45^\circ\rbrace$.
For the TPFA scheme, the results are only shown for $\alpha = 45^\circ$ because the results for $\alpha = - 45^\circ$ are identical.
Due to the flow resistance imposed by the porous medium, most of the air passes this obstacle through the constricted section above the block,
thus leading to the highest flow velocities there. While the gas passes the porous block virtually parallel for both $\alpha = 45^\circ$ and $\alpha = 0^\circ$
when applying the TPFA method, the effect of anisotropy is clearly visible in the MPFA results. Here, the flow follows the inclined principal direction
of the permeability tensor, exiting or entering the porous domain at the top. Small regions of local recirculation can be found within the obstacle which are caused by the medium's
anisotropy and the closed wall at the bottom. For $\alpha = - 45^\circ$, the upward flow in the box creates a recirculation at the right part of the bottom, where a small amount of
gas is actually pulled from the free-flow channel into the porous medium. Analogously for $\alpha = 45^\circ$, the downward flow causes a recirculation at the bottom left, where the
gas cannot exit through the solid wall and thus leaves the domain towards the left, in opposition to the general flow field.

Due to the small pressure differences, we subtract a reference value
of $\unit[p_\mathrm{ref}=10^{5}]{Pa}$ for improved visibility as shown on the right side of Figures \ref{pc:testOneVelMpfa} and \ref{pc:testOneVelTpfa}.
Almost the entire pressure drop along the channel is observed at the porous domain.
Again, the influence of anisotropy is clearly visibly for the MPFA results in terms of an inclined pressure profile which also reflects the closed bottom of the domain. None of these effects can be
captured by the TPFA method.

\begin{figure}[ht!]
	\centering
	\includegraphics[width=0.45\linewidth,keepaspectratio]{pics/mpfa_lowre_100_-45_vel.png}
	\hspace{0.04\linewidth}
	\includegraphics[width=0.45\linewidth,keepaspectratio]{pics/mpfa_lowre_100_-45_press.png} \\
	\vspace{0.02\linewidth}
	\includegraphics[width=0.45\linewidth,keepaspectratio]{pics/mpfa_lowre_100_45_vel.png}
	\hspace{0.04\linewidth}
	\includegraphics[width=0.45\linewidth,keepaspectratio]{pics/mpfa_lowre_100_45_press.png}
	\caption{Velocity (left column) and pressure (right column) profiles of MPFA scheme for an ansisotropy ratio of $\beta = 100$ and for angles $\alpha = -45^\circ$ (upper row) and $\alpha = 45^\circ$ (lower row). The reference pressure is set to $p_{\mathrm{ref}} = \unit[10^5]{Pa}$.}
	\label{pc:testOneVelMpfa}
\end{figure}

\begin{figure}[ht!]
	\centering
	\includegraphics[width=0.45\linewidth,keepaspectratio]{pics/tpfa_lowre_100_45_vel.png}
	\hspace{0.04\linewidth}
	\includegraphics[width=0.45\linewidth,keepaspectratio]{pics/tpfa_lowre_100_45_press.png}
	\caption{Velocity (left column) and pressure (right column) profiles of TPFA scheme for an ansisotropy ratio of $\beta = 100$ and for angle $\alpha = 45^\circ$. The reference pressure is set to $p_{\mathrm{ref}} = \unit[10^5]{Pa}$.}
	\label{pc:testOneVelTpfa}
\end{figure}

Figure \ref{pc:massFluxesTestOne} shows the total mass fluxes over the porous-medium free-flow interfaces. For $\alpha = 0$, the TPFA and MPFA scheme result in almost
the same solutions (small differences occur because a finer mesh is used for the TPFA scheme in the porous medium). For this angle, the TPFA scheme is also consistent
and produces the correct results. However, it is observed that for all other angles the total mass fluxes significantly differ. The off-diagonal terms of $\mathbf{K}$
are not considered in the TPFA transmissibilities, which explains why the TPFA results are independent from the direction of rotation. Furthermore, the total fluxes
over $\Gamma^\mathrm{pm}_\mathrm{top}$ are small for the TPFA scheme, in contrast to the MPFA scheme, where the total mass fluxes increase with increasing rotation angle
due to the contribution from the non-parallel porous-medium flow as mentioned above.

\begin{figure}[ht!]
	\centering
	\includegraphics[width=0.45\linewidth,keepaspectratio]{pics/testOne_fluxInOutBeta10.pdf}
	\hspace{0.04\linewidth}
	\includegraphics[width=0.45\linewidth,keepaspectratio]{pics/testOne_fluxTopBeta10.pdf}
    \includegraphics[width=0.45\linewidth,keepaspectratio]{pics/testOne_fluxInOutBeta100.pdf}
	\hspace{0.04\linewidth}
	\includegraphics[width=0.45\linewidth,keepaspectratio]{pics/testOne_fluxTopBeta100.pdf}
	\caption{Total mass fluxes of MPFA and TPFA schemes over the porous-medium boundaries $\Gamma^\mathrm{pm}_\mathrm{in}, \Gamma^\mathrm{pm}_\mathrm{out}$ (left column) and  $\Gamma^\mathrm{pm}_\mathrm{top}$ (right column) for an ansisotropy ratio of $\beta = 10$ (upper row) and $\beta = 100$ (lower row). In the left pictures, the dashed lines correspond to the mass fluxes at the inlet boundary $\Gamma^\mathrm{pm}_\mathrm{in}$, whereas the solid lines to the fluxes at the outlet boundary $\Gamma^\mathrm{pm}_\mathrm{out}$. Positive fluxes mean that fluid flows into the porous medium.}
	\label{pc:massFluxesTestOne}
\end{figure}

\subsection{Test case 2}
The next test case investigates the solution behavior for higher Reynolds numbers.
It uses a similar setting as the previous one with the difference that the channel is elongated in $x$-direction. The computational domains are given by $\Omega = [0,2.5] \times [0,0.25]$ m and $\Omega^\mathrm{pm} = [0.4,0.6] \times [0,0.2]$ m such that $\Omega^\mathrm{ff} = \Omega \setminus \Omega^\mathrm{pm}$. A pressure difference of $\Delta p = \unit[2 \cdot 10^{-3}]{Pa}$ between the left and the right boundary results in $Re \approx 130$ in the channel right atop the porous block.

\begin{figure}[ht!]
	\centering
	\includegraphics[width=0.7\linewidth,keepaspectratio]{pics/mpfa_highre_100_45_vel_time20.png}
	\includegraphics[width=0.7\linewidth,keepaspectratio]{pics/mpfa_highre_100_45_vel_time40_noColorbar.png}
	\includegraphics[width=0.7\linewidth,keepaspectratio]{pics/mpfa_highre_100_45_vel_time80_noColorbar.png}
	\includegraphics[width=0.7\linewidth,keepaspectratio]{pics/mpfa_highre_100_45_vel_time200_noColorbar.png}
	\includegraphics[width=0.7\linewidth,keepaspectratio]{pics/mpfa_highre_100_45_vel_time1000_noColorbar.png}
	\caption{Velocity profiles of MPFA scheme for an ansisotropy ratio of $\beta = 100$ and for angle $\alpha = 45^\circ$ for times $t \in \lbrace \unit[20]{s}, \unit[40]{s}, \unit[80]{s}, \unit[200]{s}, \unit[1000]{s} \rbrace$. The domain is scaled by a factor of 2 in $y$-direction for a better visualization. The porous-medium boundary is represented by the black lines. }
	\label{pc:testHighReVel}
\end{figure}

\begin{figure}[ht!]
	\centering
	\includegraphics[width=0.7\linewidth,keepaspectratio]{pics/tpfa_highre_100_45_vel_time1000.png}
	\caption{Velocity profile of TPFA scheme for an ansisotropy ratio of $\beta = 100$ and for angle $\alpha = 45^\circ$ for times $t = \unit[1000]{s}$. The domain is scaled by a factor of 2 in $y$-direction for a better visualization. The porous-medium boundary is represented by the black lines. }
	\label{pc:testHighReVelTPFA}
\end{figure}

For this test case, it takes much longer until a stationary solution is reached, which is why the solutions are investigated at different time steps.
In the following, we focus on the discussion of the case where $\alpha = 0^\circ$ or $\alpha = 45^\circ$ and $\beta = 100$. Figure \ref{pc:testHighReVel} shows the resulting velocity fields of the MPFA method at
different times for $\alpha = 45^\circ$. As before, the gas hits the front of the porous block and is forced mainly through the narrow channel section above it which acts as
some sort of duct through which a jet of high velocity fluid streams in the open channel on the right side. As the velocity gradually increases over time, the formation of vortex structures within
the free-flow channel can be observed after around $\unit[20]{s}$. Having reached equilibrium at around $\unit[1000]{s}$, the system features two stable countercurrent, larger vortices downstream
the porous obstacle and one small recirculation zone in front of the block. Again, the flow within the obstacle is clearly influenced by the anisotropy, a feature not reproduced by the TPFA method
as seen in Figure \ref{pc:testHighReVelTPFA}. Here, the flow within the block does not follow $\mathbf{K}$'s orientation. Instead, the fluid immediately strives for the top of the obstacle once it
has entered it which is due to the strongly increased vertical permeability ($\beta = 100$). Again, the TPFA scheme is not able to capture this relevant anisotropy effect as it only considers the
main diagonals of $\mathbf{K}$. Figure \ref{pc:massFluxesTestTwo} depicts the cumulative mass fluxes per unit depth across the porous medium's boundaries (top row) and across the plane at the center
of the narrow channel section ($x = \unit[0.5]{m}$, bottom row) for a case with $\alpha = 0^\circ$ (left column) and one with $\alpha = 45^\circ$ (right column). For the former, the results of
the TPFA and MPFA are very similar. With increasing time, the total mass flux entering the obstacle's front (red line) increases in accordance with the global velocity field (see Figure
\ref{pc:testHighReVel}) until it approaches a constant value after approximately $\unit[100]{s}$. The same applies for the fluxes leaving the obstacle's top (blue line), which approaches the same
value as the incoming fluxes but with a different sign. This indicates that no mass flux occurs over the obstacle's back which can also be seen from the black line. This line (and thus, the fluxes
over the obstacle's top) only deviates from zero at the very beginning of the simulation where it has to balance out the initial disequilibrium between the red and blue curve. Due to the imposed
anisotropy ratio of $\beta = 100$, the  horizontal permeability is 100 times lower than the vertical one and fluid is immediately pushed towards the top, where it exits the domain again.
Even though the fluid's density is actually pressure-dependent, significant compressibility effects could not be
observed for the given setup. The differences between TPFA and MPFA  with respect to the blue and red curves (front and back of the obstacle) are most likely due to the different
discretization width in the porous domain as explained before. In total, less gas seems to enter the porous block for the MPFA method. The plot on the lower left of Figure
\ref{pc:massFluxesTestTwo} shows the temporal evolution of the total mass flux through the channel above the obstacle. Both the TPFA and the MPFA scheme converge to the same result,
the grid resolution of the free-flow domain is identical for both cases. Note that even with a constant flux through the constriction after $\unit[100]{s}$, the global velocity field downstream the
obstacle still changes, including the formation of vortices as described before.

\begin{figure}[ht!]
	\centering
	\includegraphics[width=0.45\linewidth,keepaspectratio]{pics/testHighReFluxOverPMTimeAngle0.pdf}
	\hspace{0.04\linewidth}
    \includegraphics[width=0.45\linewidth,keepaspectratio]{pics/testHighReFluxOverPMTimeAngle45.pdf}
	\includegraphics[width=0.44\linewidth,keepaspectratio]{pics/testHighReFluxOverFFTimeAngle0.pdf}
	\hspace{0.06\linewidth}
	\includegraphics[width=0.44\linewidth,keepaspectratio]{pics/testHighReFluxOverFFTimeAngle45.pdf}
	\caption{Total mass fluxes of MPFA and TPFA schemes over the porous-medium boundaries $\Gamma^\mathrm{pm}_\mathrm{in}, \Gamma^\mathrm{pm}_\mathrm{top}, \Gamma^\mathrm{pm}_\mathrm{out}$ (upper row) and  at the center of the constricted section of the free-flow channel (over the line segment connecting the points $(\unit[0.5]{m},\unit[0.2]{m})^T$ and $(\unit[0.5]{m},\unit[0.25]{m})^T$) (lower row) for an ansisotropy ratio of $\beta = 0$ (left column) and $\beta = 100$ (right column). Positive fluxes mean that fluid flows into the porous medium.}
	\label{pc:massFluxesTestTwo}
\end{figure}

Considering the right column of Figure \ref{pc:massFluxesTestTwo} ($\alpha = 45^\circ$), the differences between the two schemes are significantly higher. For the MPFA method, we observe a
considerable inflow across the obstacle's top boundary, coming from the constricted free-flow channel and drawn by the downwards inclined flow field within the porous medium. At the same time,
there is only a very limited inflow through the obstacle's front which is due the obstacle's anistropy.

This is in strong contrast to the TPFA method's result, where air is still leaving the obstacle's top as also described before and seen in Figure \ref{pc:testHighReVelTPFA}.
The MPFA method results in a higher flux through the constricted channel (see bottom right of Figure \ref{pc:massFluxesTestTwo}) as there is less fluid entering the porous obstacle as explained
above.

\subsection{Test case 3}

Another advantage of using MPFA in the porous medium domain is the ability to use unstructured grids while maintaining
 consistency of the scheme. The test case presented in this section considers a porous medium domain in which geometrical
 constraints favor the use of triangles (unstructured) over quadrilaterals for its discretization. This situation can arise,
 for example, in environmental applications considering the exchange processes between the atmosphere and the subsurface, where
 the latter can be composed of complex shaped geological layers. An illustration of the setup of this test case is given in
 Figure \ref{pc:settingCaseTwo}, while a detailed view on the discretization of the two compartments is depicted in Figure \ref{pc:testTwoGrid}.

\begin{figure}[ht!]
	\centering
	\includegraphics[width=0.95\linewidth,keepaspectratio]{pics/setting_testtwo.pdf}
	\caption{Setting for the third test case.}
	\label{pc:settingCaseTwo}
\end{figure}

\begin{figure}[ht!]
	\centering
	\includegraphics[width=0.55\linewidth,keepaspectratio]{pics/grid_testtwo.png}
	\caption{Detailed view on the grid used in test case three.}
	\label{pc:testTwoGrid}
\end{figure}

Two sets of simulations were performed using the rotation angles $\alpha = 45^\circ$ and $\alpha = -45^\circ$, together with $\beta = 10$,
 $\unit[k=10^{-6}]{m^2}$ and $\Delta p = \unit[1 \cdot 10^{-6}]{Pa}$. The results using both MPFA and TPFA in the porous medium domain are shown
 in the Figures \ref{pc:testTwoPMVelocity_-45} and \ref{pc:testTwoPMVelocity_45}, respectively. Please note that unlike before, the velocity vectors are now scaled by magnitude.
 As expected, the solutions using TPFA exhibit
 a more distorted velocity field originating from the scheme being inconsistent on both unstructured grids and for anisotropic tensors.
 On the other hand, the solutions using MPFA provide velocity fields that follow the geometrical features of the porous medium.
 For $\alpha = -45^\circ$, two main regions in which the flux from the porous medium into the free-flow domain is concentrated,
 and where the maximum velocities occur, can be observed (see Figure \ref{pc:testTwoPMVelocity_-45}). The direction of these
 maximum velocities coincides with the direction of highest permeability. In contrast to that, the inflow from the free-flow domain into the porous
 medium occurs in the direction of the lowest permeability and thus at smaller velocities.
 In the case of $\alpha = 45^\circ$, the highest velocities are observed in the regions where an inflow from the free-flow domain into the porous
 medium occurs, again following the direction of the highest permeability (see Figure \ref{pc:testTwoPMVelocity_45}). With TPFA being used in the porous
 medium domain, the differences between the velocity fields obtained from the two angles turn out to be smaller in comparison to the solutions obtained with MPFA.
 Furthermore, the low velocity regions seem to be generally overestimated. This effect shows itself also in the integrated transfer flux across the interface, which
 results in $\unit[1.3 \cdot 10^{-8}]{kg/s}$ for MPFA and $\unit[2.2 \cdot 10^{-8}]{kg/s}$ for TPFA, thus, around $\unit[70]{\%}$ higher.

\begin{figure}[ht!]
	\centering
	\includegraphics[width=0.85\linewidth,keepaspectratio]{pics/test2_mpfa_45.png}
	\includegraphics[width=0.85\linewidth,keepaspectratio]{pics/test2_tpfa_45.png}
	\caption{Velocity distribution in the porous medium for the third test case with $\beta = 10$ and $\alpha = -45^\circ$.
             The upper and lower pictures depict the result using MPFA and TPFA in the porous medium domain, respectively.}
	\label{pc:testTwoPMVelocity_-45}

\end{figure}

\begin{figure}[ht!]
	\centering
    \includegraphics[width=0.85\linewidth,keepaspectratio]{pics/test2_mpfa_-45.png}
	\includegraphics[width=0.85\linewidth,keepaspectratio]{pics/test2_tpfa_-45.png}
	\caption{Velocity distribution in the porous medium for the third test case with $\beta = 10$ and $\alpha = 45^\circ$.
             The upper and lower pictures depict the result using MPFA and TPFA in the porous medium domain, respectively.}
	\label{pc:testTwoPMVelocity_45}

\end{figure}

\section{Conclusions}
\label{sec:conclusions}

In this work, a discretization method is proposed for problems concerning free flow coupled to
porous medium flow. The method combines the stability of the staggered grid finite volume method
for the free-flow equations with the consistency of MPFA finite volume methods for flows in
anisotropic porous media. We have shown how appropriate alignment of the grids results in a natural
coupling between the two discretization schemes due to coinciding degrees of freedom.

The stability and consistency of the method are emphasized with the use of numerical experiments.
Especially in the presence of anisotropy in the porous medium, a significant difference is
observed with respect to the widely used, but inconsistent, TPFA finite volume method. We moreover
emphasized the use of unstructured grids in the porous medium, allowing for computations on more
general geometries. The use of unstructured grids in the free-flow domain \cite{eymard2014a} is a topic for future
investigation.

Future work will moreover address the extension of the presented model to multi-phase flow simulations including compositional and non-isothermal effects,
 with a special focus on evaporative processes at the interface between the porous medium and the free-flow domain. This extension
 would make the model applicable to a large variety of applications, as e.g.\ the drying of soil due to evaporation \cite{mosthaf2011a} and subsequent soil
 salinization \cite{jambhekar2015a}. This drying process can lead to the creation of fractures within the porous medium. Therefore, we want to employ
 a discrete fracture model for the description of the porous medium domain, for which there is a preexisting implementation available
 in {Du\-Mu$^\text{x}$ } (see \cite{glaser2017discrete}). To increase efficiency, we will also use other finite volume schemes for the porous domain, as for example those that have been recently presented in \cite{Schneider.ea:2018}.

\section*{Acknowledgements}
We thank the German Research Foundation (DFG) for supporting this work by funding SFB 1313, Research Project A02.


\subsubsection*{\bibname}}

\usepackage[english]{babel}
\addto{\captionsenglish}{%
  \renewcommand{\bibname}{References}
}
\usepackage[utf8x]{inputenc}
\usepackage[T1]{fontenc}
\usepackage[noend]{algpseudocode}
\usepackage{fancyhdr,amssymb}
\usepackage[ruled,vlined]{algorithm2e}
\usepackage{array}
\usepackage{stmaryrd}
\usepackage{amsfonts}

\usepackage{amsfonts,amsmath,amssymb,amsthm,epsfig,psfrag}
\usepackage[T1]{fontenc}
\usepackage{babel}

\usepackage{amsmath}
\usepackage{mathtools}
\usepackage{dsfont}
\usepackage{systeme}
\usepackage{graphicx}
\usepackage[colorinlistoftodos]{todonotes}
\usepackage{blindtext}
\usepackage{enumitem}
\usepackage{xcolor}

\newcommand{\cmmnt}[1]{}

\newcolumntype{M}[1]{>{\centering\arraybackslash}m{#1}}
\newcolumntype{N}{@{}m{0pt}@{}}

\setlist[description]{leftmargin=\parindent,labelindent=\parindent}

\begin{document}

\maketitle

\begin{abstract}
    We address the problem of graph classification based only on structural information. Inspired by natural language processing techniques (NLP), our model sequentially embeds information to estimate class membership probabilities. Besides, we experiment with NLP-like variational regularization techniques, making the model predict the next node in the sequence as it reads it. We experimentally show that our model achieves state-of-the-art classification results on several standard molecular datasets. Finally, we perform a qualitative analysis and give some insights on whether the node prediction helps the model better classify graphs.

\end{abstract}

\section{Introduction}

Many natural or synthetic systems have a natural graph representation where entities are described through their mutual connections: chemical compounds, social or biological networks, for example. Therefore, automatic mining of such structures is useful in a variety of applications.

Graph classification raises several difficulties to leverage standard machine learning algorithms. Indeed, most of these algorithms take vectors of fixed size as inputs. In the case of graphs, usual representations such as edge list or adjacency matrix do not match this constraint. The size of the representations is graph dependent (number of edges in the first case, number of nodes squared in the second) and these representations are index dependent: up to indexing of its nodes, a same graph admits several equivalent representations. In a classification task, the label of a graph is independent from the indices of its nodes, so the model used for prediction should be invariant to node ordering as well. Handling discrete inputs with variable size and ordering is a well known problem in natural language processing (NLP). This is why we adapt NLP techniques to tackle graph classification.

In this paper, we propose a method to sequentially embed graph information in order to perform classification. By construction, this recurrent graph classifier overcomes the common difficulties listed above. Besides, we propose to use an additional node prediction block to help the model to capture the intrinsic structure of the graphs. The complete model is denoted \textit{variational recurrent graph classifier} (VRGC). Experiments show that this leads to better classification results for larger datasets. For clarity of the contribution of our work, we use neither node attributes nor edge attributes.

Related works use either graph kernels \citep{nikolentzos2017kernel, nikolentzos2017matching, nikolentzos2018degeneracy, neumann2016propagation, shervashidze2011weisfeiler, yanardag2015deep}, sequential methods \citep{callut2008classification, xu2012protein, jin2018learning, you2018graphrnn} or graph features \citep{barnett2016feature, DBLP:journals/corr/NarayananCVCLJ17, gomez2017dynamics, dutta2017high}.

\section{Model}
\label{sec:model}

We propose to use a sequential approach to embed graphs with a variable number of nodes and edges into a vector space of a chosen dimension. This latent representation is then used for classification. Node index invariance is approximated through specific pre-processing and aggregation.

Let $G=(V, E)$ be an undirected and unweighted graph with $V$ a set of nodes and $E$ a set of edges. The graph $G$ can be represented, modulo any permutation $\pi$ over its nodes $\{n_i\}_{i=1}^{|V|}$, by its boolean adjacency matrix $A^{\pi}$ such that $A^{\pi}_{ij} = 1$ if nodes indexed by $i$ and $j$ are connected in the graph and $A^{\pi}_{ij} = 0$ otherwise. We use this adjacency matrix as a raw representation of the graph.

Our VRGC is composed of three main parts: node ordering and embedding, classification and regularization with variational auto-regression (VAR). See Figure \ref{fig:macro} for an illustration.

\begin{figure*}
    \centering
    \includegraphics[width=1\textwidth]{archi_macro.png}
    \caption{Macroscopic representation of VRGC.}
    \label{fig:macro}
\end{figure*}

\paragraph{Node ordering and embedding}
\label{subsec:node ordering}

Before being processed by the neural network, the adjacency matrix of a graph is transformed on-the-fly \citep{you2018graphrnn}. First, a node is selected at random and used as root for a breadth first search (BFS) over the graph. The rows and columns of the adjacency matrix are then reordered according to the sequence of nodes returned by the BFS. Next, each row $i$ (corresponding to the $i^{th}$ node in the BFS ordering) is truncated to keep only the connections of node $n_{i}$ with the $\min(i, d)$ nodes that preceded in the BFS. This way, each node is $d$-dimensional, and each truncated matrix is zero-padded in order to have dimensions $(|V|_{max}, d_n)$. Throughout the rest of the paper, we use the notation $n_G$ for $|V|_{max}$.

After node ordering and pre-embedding, each graph is processed as a sequence of $d$-dimensional nodes by a gated recurrent unit (GRU) neural network \citep{cho2014learning}. The GRU is a special RNN able to learn long term dependencies by solving vanishing gradient effect\footnote{The choice of GRU over Long Short Term Memory networks is arbitrary as they have equivalent long-term modeling power \citep{chung2014empirical}.}. In order to help the recurrent network training, we propose to add a simple fully connected network between pre-embedding and recurrent embedding. Therefore, the node will be presented to the GRU in the shape of continuous vectors instead of binary adjacency vectors.

Finally the GRU sequentially embeds each node $n_i$ by using $n_{i-1}$ and information contained in a memory cell $h_{i-1}$ that theoretically embeds all previously seen information. The embedded node sequence $\{h_i\}_{i=1}^{n_G}$ then feeds both the VAR and the classifier as discussed in subsequent sections. See top line of Figure \ref{fig:archi}.

\paragraph{Classification}
\label{subsec:classification}

After the embedding step, we use an additional GRU dedicated to classification that takes $\{h_i\}_{i=1}^{n_G}$ as input.
Its last memory cell, denoted $\tilde{h}_{n_G}$, feeds a softmax multilayer perceptron (MLP) which performs class prediction.
Formally, let $c$ be the class index, the classifier is trained by minimizing the cross-entropy loss between ground-truth and $\hat{p}(G, r)$ the softmax class membership probability vector for a given graph $G$ that has been sorted by a BFS rooted with node $r$. We call this objective term $\mathcal{L}_{classif}$.
As discriminating patterns might be spread across the whole graph, the network is required to model long-term dependencies. By construction, GRUs have such ability. See middle line of Figure \ref{fig:archi}.

\paragraph{Regularization with variational auto-regression}
\label{subsec:vap}

As the structure of a graph is the concatenation of the interactions between all nodes and their respective neighbors, learning a good representation without using node attributes requires for the model to capture the structure of the graph while classifying. Accordingly, we add an auto-regression block to our model: at each node, the network makes a prediction for the next node adjacency. Multi-task learning is a powerful leverage to learn rich representation \citep{sanh2018hierarchical} in NLP. In particular, such representation for sequence classification has already been used for sentiment analysis \citep{latif2017variational, xu2017variational}.

We use a variational auto-encoder (VAE) \citep{kingma2013auto} to learn a representation of each node $n_i$ given $h_{i-1}$. The first layers of the encoder are shared with the classifier and corresponds to the graphs preprocessing (blue part in Figures \ref{fig:macro} and \ref{fig:archi}). The subsequent encoder layers, the latent sampling and the decoder constitute the VAR. For each graph $G$ with embedded nodes $\{ n_{i} \}_{i \in \llbracket 1, n_G \rrbracket}$, the fully connected variational auto-encoder takes $h\coloneqq\{ h_{i} \}_{i \in \llbracket 1, n_G-1 \rrbracket }$ as input. Let $\{ z_{i} \}_{i \in \llbracket 1, n_G \rrbracket}$ be the latent random variables for the following model

\begin{align*}
    p(G,z|h) = p(n_1, z_1)\prod_{i=2}^{n_G}{p_\theta(n_i|z_{i})q_\phi(z_{i}|h_{i-1})}.
\end{align*}

In practice, $p_\theta$ and $q_\phi$ are modelled by neural networks parametrized by $\theta$ and $\phi$, which require differentiable functions for training. However, $p_\theta(n_i|z_i)$ models a binary adjacency vector representing the connections between node $n_i$ and previously visited nodes $n_{j<i}$. Therefore, we use sigmoid continuous relaxation to train our model, and hard binary sampling at test time. We use a Gaussian variational posterior distribution. Training is done by maximizing the variational lower bound of the log-likelihood of the observation as in Kingma's VAE. The exact loss is displayed in Appendix \ref{app:elbo} and denoted $\mathcal{L}_{pred}$.

The regularization part is illustrated in the bottom line of Figure \ref{fig:archi}. In the end, the model is trained by minimizing the total loss $\mathcal{L} = \mathcal{L}_{classif} + \alpha \mathcal{L}_{pred}$, where $\alpha$ is a hyper-parameter.

\begin{figure*}
    \centering
    \includegraphics[width=1\textwidth]{VRGC.png}
    \caption{Architecture for VRGC. Top: node ordering and embedding. Middle: classification. Bottom: regularization with VAR plus final aggregation. FC stands for fully-connected.}
    \label{fig:archi}
\end{figure*}

\paragraph{Aggregation of the results at test time}
\label{sub:aggreg}

The node ordering step introduces randomness to our model. On the one hand, it helps learn more general graph representations during the training phase, but on the other hand, it might produce different outputs for the same graph during the testing phase, depending on the root of the BFS. In order to counter this side effect, we add the following aggregation step for the testing phase\cmmnt{ combining a soft and a hard vote}. Each graph is processed $N$ times by the model with $N$ different roots for BFS ordering. The $N$ class membership probability vectors are extracted and averaged. The average score vector is noted $\bar{p}$ and computed as follows with an element-wise sum:

$$
    \bar{p}(G) = \frac{1}{N} \sum_{\substack{i=1 \\ r \sim \mathcal{U}(\llbracket 1, n_G \rrbracket)}}^N {\hat{p}(G, r)}.
$$

This soft vote is repeated $K$ times resulting in $K$ probability vectors $\left\{\bar{p}_{.,k}(G)\right\}_{k=1}^K$ for each graph $G$. The final class attributed to a graph corresponds to the highest probability among the $K$ vectors.
This second hard vote enables to choose the batch of votes for which the model is the most confident.
\begin{align*}
    \hat{c}(G) & = \arg \max_{c\in \llbracket 1,C \rrbracket} \left\{ || \bar{p}_{c, .}(G) ||_\infty \right\}
\end{align*}

\section{Experiments}
\label{sec:experiments}

\paragraph{Datasets and results}

We evaluated our model against four standard datasets from biology: Mutag (MT), Enzymes (EZ), Proteins Full (PF) and National Cancer Institute (NCI1) \citep{KKMMN2016}. A detailed description of each dataset in provided in Appendix \ref{app:datasets}.

We compare our results to those obtained by Earth Mover's Distance \citep{nikolentzos2017matching} (EMD), Pyramid Match \citep{nikolentzos2017matching} (PM), Feature-Based \citep{barnett2016feature} (FB), Dynamic-Based Features \citep{gomez2017dynamics} (DyF), Stochastic Graphlet Embedding \citep{dutta2017high} (SGE), Truncated Laplacian Spectrum \citep{de2018simple} (TLS) and Family of Graph Spectral Distances \citep{verma2017hunt} (FGSD). All values are directly taken from the aforementioned papers as they use a setup similar to ours. For algorithms presenting results with and without node features, we reported the results without node features. For those presenting results with several sets of hyper-parameters, we reported the results for the parameters that performed best on the largest number of datasets. Results are reported in Table \ref{tab:results}. We obtain state-of-the-art results on three out of the four datasets used for this paper and the second best result on the fourth one.

\begin{table}
    \centering
    \renewcommand{\arraystretch}{1.2}
    \begin{tabular}{l|M{1.5cm} M{1.5cm} M{1.5cm} M{1.5cm} M{1.5cm} M{1.5cm}}
           & MT            &      EZ       & PF            & NCI1 \\
      \hline
      EMD  & 86.1          &     36.8      & -             & 72.7 \\
      PM   & 85.6          &     28.2      & -             & 69.7 \\
      FB   & 84.7          &     29.0      & 70.0          & 62.9 \\
      DyF  & 86.3          &     26.6      & 73.1          & 66.6 \\
      SGE  & 87.3          &     40.7      & 71.9          & - \\
      TLS  & 88.4          &     43.7      & 73.6          & 75.2 \\
      FGSD & \textbf{92.1} &     -         & 73.4          & 79.8 \\
      \hline
      RGC  & 89.5          & \textbf{48.7} & 72.5          & 78.1 \\
      VRGC & 86.3          & 48.4          & \textbf{74.8} & \textbf{80.7} \\
      \end{tabular}
  \caption{Experimental results of different models plus our own on four standard molecular datasets. RGC stands for recurrent graph classifier, VRGC for variational RGC. All other acronyms are defined in section \ref{sec:experiments}.}
  \label{tab:results}
\end{table}

\paragraph{Node indexing invariance}
Our model is designed to be independent from node ordering of the graph with respect to different BFS roots. Inputs representing the same graph (up to node ordering) should be close from one another in the latent embedding space. As the preprocessing is performed on each graph at each epoch, a same graph is processed many times by the model during training with different embeddings. This creates a natural regularization for the network. Indeed, as illustrated in Figure \ref{fig:projection}, the projections corresponding to the same graphs form a heap in the low dimensional representation of the latent space.

\paragraph{Contribution of the VAR to classification}
The variational regularization term seems to help the model finding a more meaningful latent representation for classification while graph dataset becomes larger. Note that the extra cost of training the VAR is marginal with respect to the training of the RNNs. We provide an illustration of the output of VAR block in Figure \ref{fig:reconstructions}.

\paragraph{Conclusion and room for improvement}
This paper proposed a recurrent graph classifier with variational regularization. The invariance to node indexing is greedy learned from numerous iterations on randomly rooted BFS-ordered graph. For future work, we should investigate the impact VAR capacity (number and size of hidden layers) on classification accuracy or generalization.

\paragraph{Acknowledgments} We would like to thank Thomas Bonald and Sebastien Razakarivony for their comments and help. We also would like to thank NVIDIA and its GPU Grant Program for providing the hardware we used in our experiments. This work is supported by the company Safran through the CIFRE convention 2017/1317.

\bibliographystyle{abbrvnat}