\section{Introduction}
Classical mixed finite element theory for the steady
Stokes problem
\begin{equation} \label{eq:cont:stokes}
\begin{split}
  -\nu \Delta {\bm v} + \nabla p & = \bm f, \\
   -{\rm div}\, {\bm v} & = g
\end{split}
\end{equation}
with inhomogeneous Dirichlet boundary data,
$\bm f \in \bm L^2(\Omega)$ and $g \in L^2_0(\Omega)$ emphasizes
that the divergence constraint
$-{\rm div}\,{\bm v}=g$
requires an appropriate discrete mimicking of the
{surjectivity} of the divergence operator
$\mathrm{div} \!\!: {\bm H}^1_0(\Omega) \to L^2_0(\Omega)$
in order
to guarantee optimal convergence properties, see e.g.\ \cite{MR3097958,jlmnr:sirev}.
Recently it has been stressed that
the divergence constraint in the Stokes
problem
 naturally induces a semi-norm and
corresponding equivalence classes of forces, which
require a \emph{second challenge} for an accurate
space discretization:
two forces
$\bm f_1 \in \bm L^2(\Omega)$ and
$\bm f_2 \in \bm L^2(\Omega)$ are \emph{velocity-equivalent}
\cite{2018arXiv180810711G}
\begin{equation}
  \bm f_1 \simeq \bm f_2,
\end{equation}
if they lead to the same velocity solution ${\bm v}$ in
the Stokes problem \eqref{eq:cont:stokes} --- and this
happens if and only if both forces differ by a
gradient field \cite{jlmnr:sirev,MR3824769}, i.e.,
\begin{equation}
  \bm f_1 \simeq \bm f_2  \qquad \Leftrightarrow
   \qquad \exists{\phi \in H^1(\Omega)/\mathbb{R}} : \bm f_2 = \bm f_1 + \nabla \phi.
\end{equation}
The argument is straightforward: denote by
$({\bm v}_1, p_1)$ and
$({\bm v}_2, p_2)$ the pairs of velocity and pressure
solutions corresponding to the forces $\bm f_1$ and
$\bm f_2=\bm f_1 + \nabla \phi$.
Then, the difference of the solutions
$(\delta {\bm v}, \delta p) := ({\bm v}_2 - {\bm v}_1, p_2-p_1)
\in \bm H^1_0(\Omega) \times L^2_0(\Omega)$
fulfills the incompressible Stokes equations
$-\nu \Delta (\delta {\bm v}) + \nabla (\delta p) = \nabla \phi$,
$\mathrm{div} (\delta {\bm v})=0$ with
homogeneous Dirichlet boundary data.
This problem has
the unique solution
$(\delta {\bm v}, \delta p) = (\bm 0, \phi)$, and
thus
$\bm f_1$ and $\bm f_2 = \bm f_1 + \nabla \phi$
are \emph{velocity-equivalent}
due to $\delta {\bm v} = \bm 0$.

In conclusion one observes that the
velocity solution ${\bm v}$ of \eqref{eq:cont:stokes}
is determined by the following data:
\begin{enumerate}
\item  Dirichlet boundary data,
\item the data $g$,
\item and the Helmholtz--Hodge projector of the data $\bm f$,
which is defined by
$$
  \mathbb{P}(\bm f)
    := \argmin_{\phi \in H^1(\Omega)} \| \bm f - \nabla \phi
      \|_{\bm L^2(\Omega)},
$$
\end{enumerate}
while the data term $\bm f - \mathbb{P}(\bm f)$
only influences the pressure.

The recently introduced notion \emph{pressure-robustness}
\cite{LM2016}
allows to discriminate between space discretizations
for \eqref{eq:cont:stokes}, whose discrete velocity
solutions ${\bm v}_h$ depend  on $\mathbb{P}(\bm f)$
and not on the entire data $\bm f$.  Such schemes lead to
a priori error estimates for the discrete velocity
that depend only on ${\bm v}$ and not on $({\bm v}, \frac{1}{\nu} p)$
--- as in nearly all classical mixed finite element methods
\cite{jlmnr:sirev}.

This contribution focuses now on applying the improved understanding of relevant data in the Stokes problem,
in order to derive a priori error estimates
for various discretely inf--sup stable mixed methods in
cases of \emph{minimal regularity}. A special focus is set
on a recent modified pressure-robust mixed method
\cite{LMT15, LLMS2017}, where the modification introduces a consistency
error that can be optimally estimated in a straightforward manner by $C h^k | \Delta {\bm v} |_{H^{k-1}(\Omega)}$
provided that ${\bm v} \in {\bm H}^{k-1}(\Omega)$. For the lowest--order methods $(k=1)$ this requires $\Delta {\bm v} \in {\bm L}^2(\Omega)$.
In situations of \emph{minimal regularity}, i.e.,
${\bm v} \in {\bm H}^{1+s}(\Omega)$ with $0 < s < 1$,
we provide an estimation of the consistency error
by a more sophisticated argument involving the
Helmholtz--Hodge projector of the data
$\nu^{-1}\mathbb{P}(\bm f)$. This term is obviously in ${\bm L}^2(\Omega)$, whenever it holds $\bm f \in {\bm L}^2(\Omega)$
and it is shown to be equal to $\mathbb{P}(- \Delta {\bm v})$.
Thus,
although it holds in general that
$\Delta {\bm v} \not\in {\bm L}^2(\Omega)$ one can exploit in the numerical analysis that at least the divergence--free part of
$\Delta {\bm v}$ is in ${\bm L}^2(\Omega)$.
{
This observation also leads to a seemingly new estimate
for classical mixed methods, which can be sharper than classical a priori estimates, see Theorem \ref{thm:ConsistencyErrorEstimate_classical}.}
Eventually, all classical conforming
finite element methods yield
an estimate of the form
\begin{align*}
\|\nabla ({\bm v}_h-{\bm S}_h({\bm v}))\|_{L^2(\Omega)} \le C_A  h^s \| \mathbb{P} (-\Delta {\bm v}) \|_{L^2(\Omega)}
+ \frac{C_B h}{\nu} \|\bm f - \mathbb{P} (\bm f) \|_{L^2(\Omega)},
\end{align*}
while their pressure-robust siblings allow for estimates of the form
\begin{align*}
\|\nabla ({\bm v}_h-{\bm S}_h({\bm v}))\|_{L^2(\Omega)} \le (C_1 h + C_A h^s) \|\mathbb{P} (-\Delta {\bm v})\|_{L^2(\Omega)},
\end{align*}
with $C_1 > 0$, $C_A > 0$ and $C_B > 0$ are constants that do not
depend on $h$. Note that for divergence-free conforming
methods, see e.g.\ \cite{scott:vogelius:conforming,MR3120580},
it holds $C_1=C_A=C_B=0$, but for them
the only nontrivial part of
the numerical analysis is the proof of the discrete inf-sup
stability.
Further, structurally identical results are obtained
for the classical and a modified pressure-robust
nonconforming Crouzeix--Raviart finite element method.

The rest of this paper is structured as follows. Section~\ref{sec:preliminaries} introduces the Stokes problem
as well as the framework for the modified finite element method
and the assumptions that are crucial for the theoretical results.
Section~\ref{sec:HLP} focusses on the Helmholtz--Hodge projector
and its application in stability estimates.
Section~\ref{sec:stoksprojectors} introduces the continuous and discrete Stokes projectors and their properties.
Section~\ref{sec:estimates_classical} applies the tools of the previous sections to obtain quasi-optimal estimates for classical finite element
methods that only depend on the data. Section \ref{sec:estimates_probust}
does the same for the modified pressure-robust finite element methods
where now the error is additionally independent of the pressure and the inverse of the viscosity $\nu$. Section \ref{sec:estimates_nonconforming} revisits
quasi-optimal and pressure-robust error estimates for the
nonconforming Crouzeix--Raviart finite element method.
Finally we perform some numerical experiments in Section \ref{sec:Numerics}
and compare these empirical results with the theory.

\section{Preliminary results}\label{sec:preliminaries}
This section introduces some notation,
recalls some preliminaries and formulates an assumption that
is fundamental for the presented theory.
We adopt standard space notation and denote
vector--valued functions and vector--valued function spaces
in boldface.  We use
$(\cdot,\cdot)$ to denote the $L^2$-inner product over $\Omega\subset \mathbb{R}^n$,
and by $\langle\cdot,\cdot\rangle$ the duality pairing between some Hilbert space and its dual.
We denote by $L^2_0(\Omega)$ the Hilbert-space of square-integrable scalar functions with zero average, and
\begin{align*}
{\bm H}({\rm div};\Omega) & = \{\bm w\in {\bm L}^2(\Omega):\ {\rm div}\, \bm w\in L^2(\Omega)\},\\
{\bm H}_0({\rm div};\Omega) & = \{\bm w\in {\bm H}({\rm div};\Omega):\ {\bm v}\cdot \bm n|_{\p \Omega} = 0\},
\end{align*}
where $\bm n$ denotes the outward unit normal of $\p\Omega$.

\subsection{Stokes problem and weak elliptic regularity assumption}
In the following, we study finite element methods
for the model problem:
for $\bm f \in \bm L^2(\Omega)$
seek $({\bm v}, p) \in {\bm H}^1_0(\Omega) \times L^2_0(\Omega)$ such that it holds
\begin{align}\label{eqn:Stokes_problem}
  - \nu \Delta {\bm v} + \nabla p = \bm f, \quad \text{and} \quad
  {\rm div}\, {\bm v} = 0\qquad \text{in }\Omega.
\end{align}
The extension to the more general divergence
constraint ${\rm div}\, {\bm v} = g$ with $g\in L^2_0(\Omega)$
is straightforward, and we refer to \cite{jlmnr:sirev}
for details.

A weak formulation of the problem is given by:
search for $({\bm v}, p) \in {\bm H}^1_0(\Omega) \times L^2_0(\Omega)$
such that it holds
\begin{equation} \label{eq:weak:stokes:problem}
\begin{split}
  \nu (\nabla {\bm v}, \nabla \bm w) - (p, {\rm div}\, \bm w) & =
    (\bm f, \bm w), \\
    ({\rm div}\,{\bm v}, q) & = 0
\end{split}
\end{equation}
for all $(\bm w, q) \in {\bm H}^1_0(\Omega) \times L^2_0(\Omega)$.

The space of divergence-free ${\bm H}^1_0(\Omega)$ vector fields is denoted as
\begin{equation}
{\bm V}^0 := \{\bm w\in {\bm H}^1_0(\Omega):\ \nabla \cdot \bm w = 0\}.
\end{equation}

\begin{assumption}\label{asmptn:regularity}
Throughout the paper, we assume that the Stokes problem inherits ${\bm H}^{1+s}(\Omega)\times H^s(\Omega)$ elliptic regularity for
some $s \in (0,1]$ and that
$\nu \| {\bm v} \|_{H^{1+s}(\Omega)} + \|p\|_{H^s(\Omega)}\le C_{\mathrm{ell},s} \| \bm f \|_{L^2(\Omega)}$.
\end{assumption}

\section{Helmholtz--Hodge projector}\label{sec:HLP}
According to the $\bm L^2$-orthogonal Helmholtz--Hodge decomposition
(see e.g.\ \cite{GR86})
any vector field $\bm f \in \bm L^2(\Omega)$ can be uniquely
decomposed
into
\begin{align} \label{eqn:HHLdecomposition}
  \bm f =

  \nabla \alpha + \mathbb{P}(\bm f),
\end{align}
where $\alpha \in H^1(\Omega)/\mathbb{R}$, and
\[
\mathbb{P}(\bm f)
\in \bm{L}^2_\sigma(\Omega) := \lbrace \bm w \in {\bm L}^2(\Omega) :
(\nabla q, \bm w) = 0 \text{ for all } q \in H^1(\Omega) \rbrace
\]
 is the Helmholtz--Hodge projector
of $\bm f$. Note, that the Helmholtz--Hodge projector of ${\bm f}$ is divergence--free and is the
orthogonal ${\bm L}^2$ projection
of $\bm f$ onto $\bm{L}^2_\sigma(\Omega)$,
i.e.,
\begin{align*}
  (\mathbb{P}(\bm f), \bm w) = (\bm f,\bm w) \qquad \text{for all } \bm w \in \bm{L}^2_\sigma(\Omega).
\end{align*}
Moreover, for the Stokes velocity solution ${\bm v}$
it holds
\begin{equation}\label{eqn:weak_solution}
  \nu (\nabla {\bm v},\nabla \bm w) = (\bm f, \bm w)  = (\mathbb{P}(\bm f),\bm w) \qquad \text{for all } \bm w \in {\bm V}^0.
\end{equation}
The
domain of the Helmholtz--Hodge projector
can be extended to ${\bm H}^{-1}(\Omega)$
with range in $({\bm V}^0)^*$, the space
of bounded linear functionals on ${\bm V}^0$.
Indeed, for every functional
$\bm f \in {\bm H}^{-1}(\Omega)$ the Helmholtz--Hodge projector
can be defined as the \emph{restriction} to ${\bm V}^0$, i.e.,
it holds
\begin{equation} \label{eq:HH:HminusOne}
  <\! \mathbb{P}(\bm f), \bm w \!> = \, <\! \bm f, \bm w \!>
   \qquad \text{for all } \bm w \in {\bm V}^0.
\end{equation}
Condition  \eqref{eq:HH:HminusOne} defines \emph{an extension} of
the Helmholtz--Hodge projector from $\bm L^2(\Omega)$
to ${\bm H}^{-1}(\Omega)$. Assume that the functional
$\hat{\bm f} \in {\bm H}^{-1}(\Omega)$ has a
representation
$\bm f\in {\bm L}^2(\Omega)$ with
$\bm f = \nabla \alpha + \mathbb{P}\bm f$.
Then it holds
for all $\bm w \in {\bm V}^0$
\begin{equation*}
  <\! \mathbb{P}(\hat{\bm f}), \bm w \!> = <\! \hat{\bm f}, \bm w \!> =
  (\bm f, \bm w) = (\mathbb{P}(\bm f), \bm w).
\end{equation*}

\begin{lemma} \label{lem:delta:v:reg}
Denote by $-\Delta : {\bm H}^1_0(\Omega) \to {\bm H}^{-1}(\Omega)$ via
\begin{equation} \label{eq:laplace:v}
  <\! -\Delta \bm w, \bm \psi \!> := (\nabla \bm w, \nabla \bm \psi)
  \qquad \text{for all } \bm \psi \in {\bm H}^1_0(\Omega).
\end{equation}
Then the weak velocity solution ${\bm v}$
of  \eqref{eq:weak:stokes:problem} satisfies
\begin{equation}
  \mathbb{P}(-\Delta {\bm v}) = \frac{1}{\nu} \mathbb{P}(\bm f).
\end{equation}
\end{lemma}
\begin{proof}
This follows directly from a combination of \eqref{eqn:weak_solution}
and \eqref{eq:HH:HminusOne}.
\end{proof}

Thus, although the regularity of the functional $-\Delta {\bm v}$
is not better in general than $-\Delta {\bm v} \in {\bm H}^{-1}(\Omega)$, its
divergence--free part $\mathbb{P}(-\Delta {\bm v})$ has the
better regularity ${\bm L}^2(\Omega)$.

\begin{remark}
We emphasize that Lemma \ref{lem:delta:v:reg} is of
central importance for the derivation
of \emph{pressure-robust} a priori error estimates
in case of minimal regularity.
We also stress that the quantity $\nu^{-1} \mathbb{P}(\bm f)$,
which appears naturally in the analysis of pressure-robust
methods,
does in fact \textit{not} scale with the inverse of $\nu$,
since it only depends on ${\bm v}$.
\end{remark}

An immediate consequence from Lemma
\ref{lem:delta:v:reg}
is the following result that bounds the norm of the
velocity field
by the norm of the Helmholtz--Hodge projector
of the data $\bm f$.
\begin{lemma}[Continuous stability estimate]
The exact solution of problem ~\eqref{eqn:Stokes_problem} satisfies
$$
  \|\nabla {\bm v} \|_{L^2(\Omega)} \leq \frac{C_{PF}}{\nu} \| \mathbb{P}(\bm f) \|_{L^2(\Omega)} = C_{PF}\| \mathbb{P} (\Delta {\bm v}) \|_{L^2(\Omega)}
$$
where $C_{PF}$ is the constant from the Poincar\'{e}--Friedrichs inequality.
\end{lemma}
\begin{proof}
The result follows directly from testing \eqref{eqn:weak_solution}
with $\bm w = {\bm v}$ and using the
Poincar\'{e}--Friedrichs inequality.\hfill
\end{proof}
\begin{remark}
Here, we emphasize that
the right hand side of the stability estimate is
given by a \emph{semi-norm of the data $\bm f$}.
This is a crucial point, which arguably has not been fully exploited
in  classical mixed theory \cite{MR3097958, GR86}.
\end{remark}

\section{Notation and setting for conforming finite element methods}
In the following, we introduce some notation for
the finite element methods used in this contribution.
We denote by ${\bm X}_h\times Q_h$,
a discretely inf-sup stable finite element pair
\cite{MR3097958} for the Stokes problem
with homogeneous Dirichlet boundary conditions
{with respect to a conforming, shape--regular and simplicial triangulation $\mathcal{T}_h$ with $h = \max_{T\in \mathcal{T}_h} {\rm diam}(T)$.}
The $L^2$ best approximation onto the discrete pressure space
$L^2_0(\Omega)$ is denoted by $\pi_h: L^2_0(\Omega) \to Q_h$, i.e., for all $r \in L^2_0(\Omega)$
it holds
\begin{equation}\label{eqn:L2ProjDef}
  (\pi_h r, q_h) = (r, q_h) \qquad \text{for all $q_h \in Q_h$.}
\end{equation}
{We assume that $Q_h$
has the approximation property
\begin{equation}\label{eqn:L2ProjApprox}
\|r-\pi_h r\|_{L^2(\Omega)} = \inf_{q_h\in Q_h} \|r-q_h\|_{L^2(\Omega)}\le C_{\pi_h, s} h^{s}\|r\|_{H^{s}(\Omega)
\end{equation}
for all $r\in H^{s}(\Omega)\cap L^2_0(\Omega)$ and $s\in (0,1]$.
}

Let $\mathrm{div}_h:{\bm X}_h \to Q_h$
with $\mathrm{div}_h = \pi_h {\rm div}\,$
denote
the discrete divergence operator.
Due to the assumed discrete inf--sup stability of the pair ${\bm X}_h\times Q_h$,
$\mathrm{div}_h$ is surjective with bounded right--inverse \cite{MR3097958}.
We define
the space of
discretely divergence--free functions as
\begin{equation*}
{\bm V}^0_h  :=\{{\bm v}_h\in {\bm X}_h : \mathrm{div}_h {\bm v}_h = 0 \}.
\end{equation*}

\subsection{Some modified finite element methods}
As shown in  \cite{Lin14, LMT15, LLMS2017},
a certain modification of the discrete right--hand side of
the incompressible Stokes problem  renders
 inf-sup stable mixed   methods
{pressure-robust}.
These {pressure-robust} finite element methods
employ a reconstruction operator with the properties stated in the following assumption.

\begin{assumption}\label{asmptn:reconstructionoperator}
We assume that there exists
an auxiliary finite element space $\bm Y_h\subset {\bm H}_0({\rm div};\Omega)$
and a  reconstruction operator ${\bm I}_h:{\bm H}^1_0(\Omega)\to \bm Y_h$ such that
\begin{align}\label{eqn:IhApprox_divfree}
& (i) &
{\rm div}\, ({\bm I}_h {\bm v}_h) & = \mathrm{div}_h {\bm v}_h & \text{for all } {\bm v}_h\in {\bm X}_h,\\
& (ii) & \|{\bm v}_h - {\bm I}_h {\bm v}_h\|_{L^2(\Omega)} &\le C_1 h \|\nabla {\bm v}_h\|_{L^2(\Omega)} & \text{for all } {\bm v}_h\in {\bm X}_h,\label{eqn:IhApprox}
\end{align}
where $C_1$ depends only on the shape regularity of the mesh.
\end{assumption}

The modified finite element method for the Stokes problem applies
the reconstruction operator in the right-hand side.
The resulting scheme  seeks $({\bm v}_h,p_h)\in {\bm X}_h\times Q_h$ such that
\begin{alignat}{2}
\nu (\nabla {\bm v}_h,\nabla \bm w_h) - (\mathrm{div}_h \bm w_h,p_h) & = ({\bm f},{\bm I}_h \bm w_h)\qquad &&\text{for all } \bm w_h\in {\bm X}_h,\label{eqn:discrete_problem}\\
(\mathrm{div}_h {\bm v}_h,q_h) & = 0 \qquad &&\text{for all } q_h\in Q_h.\nonumber
\end{alignat}

Testing \eqref{eqn:discrete_problem} with
discretely divergence-free velocity test functions yields
\begin{equation}
  \label{eqn:weak_solution_discrete}
  \nu  (\nabla {\bm v}_h,\nabla \bm w_h) = (\bm f, {\bm I}_h \bm w)  = (\mathbb{P}(\bm f), {\bm I}_h \bm w_h) \qquad  \text{for all } \bm w_h \in {\bm V}^0_h,
\end{equation}
since for $\bm w_h \in {\bm V}^0_h$ it holds ${\bm I}_h \bm w_h \in \bm{L}^2_\sigma(\Omega)$.
This last identity is characteristic for pressure-robustness
and in general \emph{not} true for non-divergence-free classical finite element methods. It tells us that the discrete velocity solution
${\bm v}_h$ of \eqref{eqn:discrete_problem} depends on the appropriate continuous data
$\nu^{-1} \mathbb{P}(\bm f)$ of the problem.

In the case of discontinuous pressure spaces $Q_h$, the
standard interpolation operators of the
Raviart-Thomas or Brezzi--Douglas--Marini finite element spaces can be employed as a reconstruction operator ${\bm I}_h$,
see \cite{jlmnr:sirev,LM2016, LMT15} for details.
For instance, in the case of the Bernardi--Raugel finite element method \cite{BernardiRaugel}, the standard interpolator into the BDM space of order one can be used.
For continuous pressure spaces, the design of the reconstruction operator is more involved; see \cite{LLMS2017}
for details in case of the Taylor--Hood or MINI finite element family.

\begin{remark}
Note, that for ${\bm I}_h = \bm 1$ (the identity operator) in \eqref{eqn:discrete_problem} the classical finite element method is obtained. However,
only divergence-free $H^1$-conforming classical finite element methods, see e.g.\ \cite{scott:vogelius:conforming,MR3120580}, satisfy Assumption~\ref{asmptn:reconstructionoperator}
with $C_1 = 0$. In the results below it will be
specified which results rely on this assumption.
\end{remark}

\begin{lemma}[Discrete stability estimates]
Let $({\bm v}_h,p_h)\in {\bm X}_h\times Q_h$ satisfy \eqref{eqn:discrete_problem}
and write $\bm f = \nabla \alpha+\mathbb{P}\bm f$.
Then if the discrete scheme satisfies Assumption \ref{asmptn:reconstructionoperator}, it holds the estimate
$$
  \|\nabla {\bm v}_h \|_{L^2(\Omega)} \leq  (C_{PF} + C_1h) \| \nu^{-1} \mathbb{P}(\bm f) \|_{L^2(\Omega)}
  = (C_{PF} + C_1h) \| \mathbb{P} (\Delta {\bm v}) \|_{L^2(\Omega)}.
$$
If the discrete scheme with ${\bm I}_h = \bm 1$ does not satisfy Assumption \ref{asmptn:reconstructionoperator}, it only holds
\begin{align*}
\|\nabla {\bm v}_h \|_{L^2(\Omega)}
&\leq C_{PF} \| \nu^{-1} \mathbb{P}(\bm f) \|_{L^2(\Omega)}
+ \frac{1}{\nu} \| \alpha - \pi_h \alpha \|_{L^2(\Omega)}\\
&\leq C_{PF} \| \mathbb{P} (\Delta {\bm v}) \|_{L^2(\Omega)}
+ \frac{C h}{\nu} \| \bm f - \mathbb{P}(\bm f) \|_{L^2(\Omega)}.
\end{align*}
\end{lemma}
\begin{proof}
Testing \eqref{eqn:weak_solution_discrete} with $\bm w_h = {\bm v}_h$, a discrete Poincar\'e--Friedrichs inequality and \eqref{eqn:IhApprox} yield
\begin{align*}
\nu \|\nabla {\bm v}_h \|_{L^2(\Omega)}^2
= (\mathbb{P}(\bm f), {\bm I}_h {\bm v}_h )
& \leq \| \mathbb{P}(\bm f) \|_{L^2(\Omega)}
\left( \| {\bm v}_h \|_{L^2(\Omega)} + \| {\bm v}_h - {\bm I}_h {\bm v}_h\|_{L^2(\Omega)} \right)\\
& \leq (C_{PF} + C_1h) \| \mathbb{P}(\bm f) \|_{L^2(\Omega)} \|\nabla {\bm v}_h \|_{L^2(\Omega)}.
\end{align*}
If ${\bm I}_h = \bm 1$ and Assumption \ref{asmptn:reconstructionoperator} is not satisfied then inserting the Helmholtz--Hodge decomposition of $\bm f$
and an integration by parts give
\begin{align*}
\nu \|\nabla {\bm v}_h \|_{L^2(\Omega)}^2
& = (\mathbb{P}(\bm f), {\bm v}_h ) + (\nabla \alpha, {\bm v}_h)\\
& = (\mathbb{P}(\bm f), {\bm v}_h )
 -(\alpha, {\rm div}\, {\bm v}) \\
& = (\mathbb{P}(\bm f), {\bm v}_h )
 -(\alpha - \pi_h \alpha, {\rm div}\, {\bm v}) \\
& \leq \left(C_{PF} \| \mathbb{P}(\bm f) \|_{L^2(\Omega)}
+ \| \alpha - \pi_h \alpha \|_{L^2(\Omega)}
 \right) \| \nabla {\bm v}_h \|_{L^2(\Omega)}.
\end{align*}
Property \eqref{eqn:L2ProjApprox}
shows $\| \alpha - \pi_h \alpha \|_{L^2(\Omega)} \leq C h \| \nabla \alpha \|_{L^2(\Omega)} = C h \| \bm f - \mathbb{P}(\bm f) \|_{L^2(\Omega)}$.
This concludes the proof.
\end{proof}

\section{Continuous and discrete Stokes projectors}\label{sec:stoksprojectors}
In preparation for the a priori error estimates,
this section studies the continuous and the discrete Stokes projectors. They are defined as the ${\bm H}^1$-seminorm
best-approximations into the (discretely) divergence-free functions, i.e.\ ${\bm S}_h:{\bm H}^1_0(\Omega)\to {\bm V}^0_h$
and ${\bm S}:{\bm X}_h\to {\bm V}^0$ are defined by
\begin{alignat}{2}
\label{eqn:DefSh}
(\nabla {\bm S}_h({\bm v}),\nabla \bm w_h) &= (\nabla {\bm v},\nabla \bm w_h)\qquad &&\forall \bm w_h\in {\bm V}^0_h,\\
\label{eqn:DefS}
(\nabla {\bm S}({\bm v}_h),\nabla \bm w) & = (\nabla {\bm v}_h,\nabla \bm w)\qquad &&\forall \bm w\in {\bm V}^0.
\end{alignat}
The rest of this section collects useful properties of these projectors.

\begin{lemma}[Stokes projector identity]\label{lem:StokesProjectoridentity}
For any ${\bm v} \in {\bm H}^1_0(\Omega)$ and ${\bm v}_h \in {\bm X}_h$, it holds the identity
\[
(\nabla {\bm S}_h({\bm v}),\nabla \bm w_h) = (\nabla {\bm v},\nabla {\bm S}(\bm w_h))\qquad \forall {\bm v}\in {\bm V}^0,\ \bm w_h\in {\bm V}^0_h.
\]
\end{lemma}
\begin{proof}
  This follows directly from the combination of the definitions of ${\bm S}_h$ and ${\bm S}$.
\end{proof}

\begin{lemma}\label{lem:Duality}
Suppose that the Stokes problem
satisfies Assumption~\ref{asmptn:regularity}.
Then there holds
\begin{align}\label{eqn:DualityEstimate}
\|\bm w_h- {\bm S}(\bm w_h)\|_{L^2(\Omega)}\le C_2 h^s\|\nabla \cdot \bm w_h\|_{L^2(\Omega)}\qquad \forall \bm w_h\in {\bm V}^0_h.
\end{align}
\end{lemma}
\begin{proof}
Let $(\bm \psi,r)\in {\bm H}^1_0(\Omega)\times L^2_0(\Omega)$ solve
the Stokes problem with source $\bm w_h - {\bm S}(\bm w_h)$ and unit viscosity:
\begin{alignat*}{2}
(\nabla \bm \psi,\nabla \bm z) - (\nabla \cdot \bm z,r) & = (\bm w_h -{\bm S}(\bm w_h),\bm z)\qquad &&\forall \bm z\in {\bm H}^1_0(\Omega),\\
(\nabla \cdot \bm \psi,q) & = 0\qquad &&\forall q\in L_0^2(\Omega).
\end{alignat*}
Testing the first equation with $\bm z = \bm w_h - {\bm S}(\bm w_h)$
and employing \eqref{eqn:DefS} leads to
\begin{align*}
\|\bm w_h - {\bm S}(\bm w_h)\|_{L^2(\Omega)}^2
&= (\nabla \bm \psi,\nabla (\bm w_h-{\bm S}(\bm w_h))) - (\nabla \cdot (\bm w_h-{\bm S}(\bm w_h)),r)\\
&= - (\nabla \cdot \bm w_h,r).
\end{align*}
Recall that $\pi_h r$ is the $L^2$--projection of $r$ defined by \eqref{eqn:L2ProjDef},
and note that
 it holds $(\nabla \cdot \bm w_h,\pi_h{r}_h)=0$
since $\bm w_h\in {\bm V}^0_h$.  Consequently, by \eqref{eqn:L2ProjApprox}, we have
\begin{align*}
\|\bm w_h - {\bm S}(\bm w_h)\|_{L^2(\Omega)}^2
&= -(\nabla \cdot \bm w_h,r-\pi_h{r}) \le \|\nabla \cdot \bm w_h\|_{L^2(\Omega)} \|r-\pi_h{r}\|_{L^2(\Omega)}\\
&\le C_{\pi_h,s} h^s \|\nabla \cdot \bm w_h\|_{L^2(\Omega)} \|r\|_{H^s(\Omega)}.
\end{align*}
Finally, the elliptic regularity Assumption~\ref{asmptn:regularity} implies
$\|r\|_{H^s(\Omega)}\le C_{\mathrm{ell},s} \|\bm w_h-{\bm S}(\bm w_h)\|_{L^2(\Omega)}$, and so
\begin{align*}
\|\bm w_h - {\bm S}(\bm w_h)\|_{L^2(\Omega)}^2
&\le C_{\mathrm{ell},s} C_{\pi_h, s} h^s \|\nabla \cdot \bm w_h\|_{L^2(\Omega)} \|\bm w_h-{\bm S}(\bm w_h)\|_{L^2(\Omega)}.
\end{align*}
Dividing the last inequality by $ \|\bm w_h-{\bm S}(\bm w_h)\|_{L^2(\Omega)}$ gets the desired result.
\end{proof}

\section{Quasi-optimal a priori error estimates for classical
finite element methods}\label{sec:estimates_classical}
This section derives a priori error estimates
for classical finite element methods that are not pressure-robust, i.e.\
do not satisfy Assumption~\ref{asmptn:reconstructionoperator} with ${\bm I}_h = \bm 1$ like
the Bernardi--Raugel, MINI or Taylor--Hood finite element methods. The proof of the
estimate bounds the error of the best-approximation
by the right-hand side data.

\begin{theorem}\label{thm:ConsistencyErrorEstimate_classical}
Suppose that the Stokes problem satisfies Assumption~\ref{asmptn:regularity}, the reconstruction operator
is taken to be the identity ${\bm I}_h = \bm 1$, and that ${\bm I}_h$ does not satisfy
Assumption~\ref{asmptn:reconstructionoperator}.
Then there holds
\begin{align*}
\|\nabla ({\bm v}_h-{\bm S}_h({\bm v}))\|_{L^2(\Omega)} &
 \le  C_2  h^s \| \nu^{-1} \mathbb{P} (\bm f) \|_{L^2(\Omega)}
+ \frac{C_{\pi_h,1} h}{\nu}  \|\bm f - \mathbb{P}(\bm f)\|_{L^2(\Omega)} \\
 & = C_2  h^s \| \mathbb{P} (-\Delta {\bm v}) \|_{L^2(\Omega)}
+ \frac{C_{\pi_h,1} h}{\nu}  \|\bm f - \mathbb{P}(\bm f)\|_{L^2(\Omega)}
\end{align*}
with $C_2>0$ given by \eqref{eqn:DualityEstimate}.
\end{theorem}
\begin{proof}
Write ${\bm e}_h := {\bm v}_h-{\bm S}_h({\bm v})$ and note that ${\bm e}_h\in {\bm V}^0_h$. Hence,
it follows from Lemmas~\ref{lem:StokesProjectoridentity} and  \ref{lem:Duality} that
\begin{align}
\|\nabla {\bm e}_h\|_{L^2(\Omega)}^2
&= (\nabla {\bm v}_h,\nabla {\bm e}_h) - (\nabla {\bm S}_h({\bm v}),\nabla {\bm e}_h) \nonumber \\
&= (\nabla {\bm v}_h,\nabla {\bm e}_h) - (\nabla {\bm v},\nabla {\bm S}({\bm e}_h)) \nonumber \\
& = \nu^{-1} (\bm f,{\bm e}_h-{\bm S}({\bm e}_h)) \nonumber \\
& = \nu^{-1} (\mathbb{P}(\bm f) + \nabla \alpha,{\bm e}_h-{\bm S}({\bm e}_h)) \nonumber \\
& = (\mathbb{P} (-\Delta {\bm v}),{\bm e}_h-{\bm S}({\bm e}_h))
- \nu^{-1} (\alpha - \pi_h \alpha,\mathrm{div}({\bm e}_h-{\bm S}({\bm e}_h))) \label{eq:rem:ref} \\
& \le
\left( C_2  h^s \| \mathbb{P} (-\Delta {\bm v}) \|_{L^2(\Omega)}
+ \nu ^{-1} \|\alpha - \pi_h \alpha \|_{L^2(\Omega)}
\right) \|\nabla {\bm e}_h\|_{L^2(\Omega)}, \nonumber
\end{align}
where $\alpha$ stems from the Helmholtz--Hodge decomposition
\eqref{eqn:HHLdecomposition} of $\bm f$.
The best approximation property of $\pi_h \alpha$
shows $\| \alpha - \pi_h \alpha \|_{L^2(\Omega)} \leq C_{\pi_h,1} h \| \nabla \alpha \|_{L^2(\Omega)} = C_{\pi_h,1} h \| \bm f - \mathbb{P}(\bm f) \|_{L^2(\Omega)}$.
This concludes the proof.
\end{proof}

\begin{remark}
Classical results for conforming
mixed methods \cite{GR86} show the a priori estimate
\[
\|\nabla ({\bm v}_h-{\bm S}_h({\bm v}))\|_{L^2(\Omega)} \leq
\frac{1}{\nu} \inf_{q_h\in Q_h} \|p-q_h\|_{L^2(\Omega)},
\]
which scales like $\nu^{-1} h^s$ under the given regularity assumptions.
Such an estimate is sometimes sharper than
Theorem \ref{thm:ConsistencyErrorEstimate_classical},
but can also be less sharp. \\
i) If it holds, e.g.,
$p \in Q_h$, then the error on the right hand side
of the classical estimate is zero.
This is also preserved in the computations for the new estimate
until \eqref{eq:rem:ref}, since in the special case
$\bm f = -\nu \Delta {\bm v}$ \eqref{eq:rem:ref}
can be shown to vanish identically. \\
ii) If it holds $p \not\in Q_h$ and if the solution
$({\bm v}, p) \in {\bm H}^{1+s}(\Omega) \times H^s(\Omega)$ has a low regularity with $s < 1$, then the new estimate can be sharper
e.g.\ for $\nu \ll 1$,
since it predicts an a priori error $\mathcal{O}(h^s + \nu^{-1} h)$, while the classical estimate predicts an error decay like
$\mathcal{O}(\nu^{-1} h^s)$.
We remark that
the pressure-dependent consistency error is influenced by two different contributions,
one determined by $-\Delta {\bm v}$ and another one
determined by
$\frac{1}{\nu} ({\bm f} - \mathbb{P}({\bm f}))$.
\end{remark}

\begin{theorem}[A priori error estimate]
  Under the assumptions of Theorem~\ref{thm:ConsistencyErrorEstimate_classical}, it holds
  \begin{multline*}
    \|\nabla({\bm v}-{\bm v}_h)\|^2_{L^2(\Omega)} \le (1+C_F) \inf_{\bm w_h\in {\bm X}_h} \|\nabla ({\bm v}-\bm w_h)\|^2_{L^2(\Omega)}\\
 +  \left(C_2  h^s \| \mathbb{P} (-\Delta {\bm v}) \|_{L^2(\Omega)}
+ \frac{C_{\pi_h,1} h}{\nu}  \|\bm f - \mathbb{P}(\bm f)\|_{L^2(\Omega)}\right)^2.
    \end{multline*}
\end{theorem}
\begin{proof}
  The proof starts with the Pythagoras theorem (using \eqref{eqn:DefSh})
  \begin{align*}
    \|\nabla ({\bm v}-{\bm v}_h)\|_{L^2(\Omega)}^2 = \|\nabla ({\bm v}-{\bm S}_h({\bm v}))\|_{L^2(\Omega)}^2
    + \|\nabla ({\bm v}_h-{\bm S}_h({\bm v}))\|_{L^2(\Omega)}^2.
  \end{align*}
  The second term can be estimated by Theorem~\ref{thm:ConsistencyErrorEstimate_classical}
  and the first term can be bounded by the best-approximation error in ${\bm X}_h$
  by the standard argument
\[
\|\nabla ({\bm v}-{\bm S}_h({\bm v}))\|_{L^2(\Omega)}\le \inf_{\bm w_h\in {\bm V}^0_h} \|\nabla ({\bm v}-\bm w_h)\|_{L^2(\Omega)}\le (1+C_F) \inf_{\bm w_h\in {\bm X}_h} \|\nabla ({\bm v}-\bm w_h)\|_{L^2(\Omega)},
\]
where $C_F \geq 1$ denotes the stability constant
of the Fortin operator of the mixed method, see e.g.\
\cite{jlmnr:sirev, GR86}.
\end{proof}

\section{Quasi-optimal pressure-robust a priori error estimates}
\label{sec:estimates_probust}

This section concerns novel quasi-optimal a priori error estimates
for conforming divergence-free and pressure-robustly modified finite element methods. Here, the distance between the discrete solution and the
discrete Stokes projector can be bounded by $\|\mathbb{P} (-\Delta {\bm v})\|_{L^2(\Omega)}$ which is in general much smaller than the bound in Theorem~\ref{thm:ConsistencyErrorEstimate_classical}, especially for small $\nu$.

\begin{theorem}\label{thm:ConsistencyErrorEstimate}
Suppose that the Stokes problem satisfies Assumption~\ref{asmptn:regularity} and that the reconstruction
operator ${\bm I}_h$ satisfies Assumption~\ref{asmptn:reconstructionoperator}.
Then there holds
\begin{align*}
\|\nabla ({\bm v}_h-{\bm S}_h({\bm v}))\|_{L^2(\Omega)} \le (C_1 h +C_2 h^s) \|\mathbb{P}(- \Delta {\bm v})\|_{L^2(\Omega)},
\end{align*}
with $C_1>0$ and $C_2>0$ given by \eqref{eqn:IhApprox} and \eqref{eqn:DualityEstimate}, respectively. Note, that there is
no dependency on $\nu^{-1}$.
\end{theorem}
\begin{proof}
Write ${\bm e}_h := {\bm v}_h-{\bm S}_h({\bm v})$ and note that ${\bm e}_h\in {\bm V}^0_h$. Hence,
\begin{align*}
(\nabla {\bm v}_h,\nabla {\bm e}_h)
&= \frac1{\nu} ({\bm f},{\bm I}_h {\bm e}_h) = (\mathbb{P}(-\Delta {\bm v}),{\bm I}_h {\bm e}_h)\\
& = (\mathbb{P}(-\Delta {\bm v}),{\bm I}_h {\bm e}_h-{\bm e}_h)+(\mathbb{P}(-\Delta {\bm v}),{\bm e}_h).
\end{align*}

The latter term is split up into (using also Lemma~\ref{lem:StokesProjectoridentity})
\begin{align*}
(\mathbb{P}(-\Delta {\bm v}),{\bm e}_h) & = (\mathbb{P}(-\Delta {\bm v}),{\bm e}_h-{\bm S}({\bm e}_h)) + (\mathbb{P}(-\Delta {\bm v}),{\bm S}({\bm e}_h))\\
& = (\mathbb{P}(-\Delta {\bm v}),{\bm e}_h-{\bm S}({\bm e}_h))+ (\nabla {\bm v},\nabla {\bm S}({\bm e}_h))\\
& = (\mathbb{P}(-\Delta {\bm v}),{\bm e}_h - {\bm S}({\bm e}_h)) + (\nabla {\bm S}_h({\bm v}),\nabla {\bm e}_h).
\end{align*}

It then follows from Lemma  \ref{lem:Duality}
and \eqref{eqn:IhApprox} that
\begin{align*}
\|\nabla {\bm e}_h\|_{L^2(\Omega)}^2
&= (\nabla {\bm v}_h,\nabla {\bm e}_h) - (\nabla {\bm S}_h({\bm v}),\nabla {\bm e}_h)\\
& = (\mathbb{P}(-\Delta {\bm v}),{\bm I}_h {\bm e}_h-{\bm e}_h)+(\mathbb{P}(-\Delta {\bm v}),{\bm e}_h-{\bm S}({\bm e}_h))\\
&\le \|\mathbb{P}(-\Delta {\bm v})\|_{L^2(\Omega)} \big( \|{\bm I}_h {\bm e}_h - {\bm e}_h\|_{L^2(\Omega)}+ \|{\bm e}_h-{\bm S}({\bm e}_h)\|_{L^2(\Omega)}\big)\\
&\le (C_1 h +C_2 h^s) \|\mathbb{P}(-\Delta {\bm v})\|_{L^2(\Omega)} \|\nabla {\bm e}_h\|_{L^2(\Omega)}.
\end{align*}
This concludes the proof.
\end{proof}

\begin{theorem}[A priori error estimate]
  Under the assumptions of Theorem~\ref{thm:ConsistencyErrorEstimate}, it holds
  \begin{align*}
    \|\nabla({\bm v}-{\bm v}_h)\|^2_{L^2(\Omega)} \le (1+C_F) \inf_{\bm w_h\in {\bm X}_h} \|\nabla ({\bm v}-\bm w_h)\|^2_{L^2(\Omega)} +  \left((C_1 h +C_2 h^s) \|\mathbb{P} (-\Delta {\bm v})\|_{L^2(\Omega)}\right)^2.
    \end{align*}
\end{theorem}
\begin{proof}
  The proof starts with the Pythagoras theorem (using \eqref{eqn:DefSh})
  \begin{align*}
    \|\nabla ({\bm v}-{\bm v}_h)\|_{L^2(\Omega)}^2 = \|\nabla ({\bm v}-{\bm S}_h({\bm v}))\|_{L^2(\Omega)}^2
    + \|\nabla ({\bm v}_h-{\bm S}_h({\bm v}))\|_{L^2(\Omega)}^2.
  \end{align*}
  The second term can be estimated by Theorem~\ref{thm:ConsistencyErrorEstimate}
  and the first term can be bounded by the best-approximation error in ${\bm X}_h$
  by the standard argument
\[
\|\nabla ({\bm v}-{\bm S}_h({\bm v}))\|_{L^2(\Omega)}\le \inf_{\bm w_h\in {\bm V}^0_h} \|\nabla ({\bm v}-\bm w_h)\|_{L^2(\Omega)}\le (1+C_F) \inf_{\bm w_h\in {\bm X}_h} \|\nabla ({\bm v}-\bm w_h)\|_{L^2(\Omega)},
\]
where $C_F \geq 1$ denotes the stability constant
of the Fortin operator of the mixed method,
see e.g.\ \cite{jlmnr:sirev, GR86}.
\end{proof}

\section{Estimates for the nonconforming Crouzeix--Raviart finite element method}
\label{sec:estimates_nonconforming}
In this section we
consider the space ${\bm X}_h \not\subset {\bm H}^1_0(\Omega)$ of nonconforming Crouzeix-Raviart functions, i.e., piecewise affine vector fields that are weakly
 continuous across edges  (2D)
or faces (3D) in the triangulation, see e.g.\ \cite{crouzeix:raviart:1973,crforty:brenner}.
To describe this space in detail we require some notation.
Recall that $\mathcal{T}_h$ is a conforming, shape--regular, and simplicial triangulation
of $\Omega$ parameterized by $h= \max_{T\in \mathcal{T}_h} {\rm diam}(T)$.
We denote by $\mathcal{E}_h$ the set of $(n-1)$--dimensional simplices in $\mathcal{T}_h$,
i.e., $\mathcal{E}_h$ is either the set of edges (2D) or faces (3D) in $\mathcal{T}_h$.
Let $P_m(T)$ denote the space of polynomials of degree $\le m$ on $T$,
and let ${\bm P}_m(T) = (P_m(T))^n$.
Then the Crouzeix-Raviart  space ${\bm X}_h$ consists of all functions $\bm w_h\in {\bm L}^2(\Omega)$
with the properties $\bm w_h|_T\in {\bm P}_1(T)$, $\int_E \bm w_h$
is single--valued for all $E\in \mathcal{E}_h$, and $\int_E \bm w_h =0$
for all boundary $E\in \mathcal{E}_h$.
The discrete pressure space $Q_h$ is the space of piecewise constants
with vanishing mean.  It is well--known that the pair ${\bm X}_h\times Q_h$
is inf--sup stable.

Note that Crouzeix-Raviart functions $\bm w_h \in {\bm V}^0_h$ are not divergence-free in a ${\bm H}(\mathrm{div})$-sense (as their normal traces are not continuous), but their piecewise divergence
vanishes.
Possible ${\bm H}(\mathrm{div})$-conforming reconstruction operators ${\bm I}_h$ for this
method are the lowest-order Raviart--Thomas or BDM  interpolation
operators, see \cite{BLMS15} for details.

In order to show the same quasi-optimal a priori error estimates for the Crouzeix--Raviart method some arguments have to be slightly modified.
First, the Stokes projectors ${\bm S}_h:{\bm H}^1_0(\Omega)\to {\bm V}^0_h$
and ${\bm S}:{\bm X}_h\to {\bm V}^0$ are now defined by
using the piecewise gradients $\nabla_h$, i.e.,
\begin{alignat}{2}
\label{eqn:DefShCR}
(\nabla_h {\bm S}_h({\bm v}),\nabla_h \bm w_h) &= (\nabla {\bm v},\nabla_h \bm w_h)\qquad &&\forall \bm w_h\in {\bm V}^0_h,\\
\label{eqn:DefSCR}
(\nabla {\bm S}({\bm v}_h),\nabla \bm w) & = (\nabla_h {\bm v}_h,\nabla \bm w)\qquad &&\forall \bm w\in {\bm V}^0.
\end{alignat}
Recall the Crouzeix--Raviart Fortin interpolation
\begin{align*}
  {\bm I}_\text{CR} {\bm v} \in {\bm X}_h
  \quad \text{defined by} \quad \int_E {\bm I}_\text{CR} {\bm v}  = \int_E {\bm v}  \quad \text{for all } E \in \mathcal{E}_h,
\end{align*}
which satisfies the approximation property
\begin{equation}\label{eqn:CRApprox}
\|\nabla_h ({\bm v}-{\bm I}_\text{CR} {\bm v})\|_{L^2(\Omega)}\le C_\text{CR} h^s \|{\bm v}\|_{H^{1+s}(\Omega)}
\end{equation}
for all $s \in [0,1]$.
This definition of the interpolant yields the well--known property \cite{crforty:brenner}
\begin{align*}
  \int_T \nabla ({\bm v} - {\bm I}_\text{CR} {\bm v}) = 0
  \quad \text{for all } T \in \mathcal{T}_h
\end{align*}
and in particular $\int_T \mathrm{div} ({\bm v} - {\bm I}_\text{CR} {\bm v})  = 0$
for any $T \in \mathcal{T}_h$. Since $\nabla_h \bm w_h$ is piecewise constant this also reveals that we have ${\bm S}_h = {\bm I}_\text{CR}$, i.e.,
the Crouzeix--Raviart interpolator is the discrete Stokes projector. Also note that the Stokes projector identity holds in the form
\[
(\nabla_h {\bm S}_h({\bm v}),\nabla_h \bm w_h) = (\nabla {\bm v},\nabla {\bm S}(\bm w_h))\qquad \forall {\bm v}\in {\bm V}^0,\ \bm w_h\in {\bm V}^0_h.
\]
However, in general ${\bm I}_\text{CR} {\bm v} \in {\bm X}_h$ does not imply ${\bm I}_\text{CR} {\bm v} \in {\bm H}(\mathrm{div},\Omega)$ and therefore
Lemma~\ref{lem:Duality} has to be modified as well.

The analysis also needs another mapping that projects a discretely divergence--free Crouzeix--Raviart function
to some $H^1$-conforming divergence-free function. Such an operator was introduced in \cite{MR3787384} and is based on rational bubble functions.

\begin{lemma}\label{lem:DualityCR}
Suppose that the Stokes problem satisfies Assumption~\ref{asmptn:regularity}.
Then there holds
\begin{align}\label{eqn:DualityEstimateCR}
\|\bm w_h- {\bm S}(\bm w_h)\|_{L^2(\Omega)}\le C_3 h^s\|\nabla_h \bm w_h\|_{L^2(\Omega)}\qquad \forall \bm w_h\in {\bm V}^0_h.
\end{align}
\end{lemma}
\begin{proof}
Consider the $H^1_0$-conforming and $H^1$-stable operator $\bm{E}_h$ from \cite{MR3787384} with
the properties
\begin{subequations}
\label{eqn:EhProperties}
\begin{alignat}{2}
  \nabla \cdot (\bm{E}_h \bm w_h) & = 0 \quad &&\text{for all } \bm w_h \in {\bm V}^0_h,\\  (\nabla_h \bm{u}_h, \nabla (\bm{E}_h \bm w_h - \bm w_h)) & = 0 \quad &&\text{for all } \bm w_h, \bm{u}_h \in {\bm X}_h,\\
 \|\nabla \bm{E}_h \bm w_h\|_{L^2(\Omega)}+ h^{-1}\|\bm{E}_h \bm w_h - \bm w_h\|_{L^2(\Omega)}&\le C_{E_h} \|\nabla_h \bm w_h\|_{L^2(\Omega)}\quad &&\text{for all }\bm w_h\in {\bm X}_h.
\end{alignat}
\end{subequations}
The second property follows from \cite[$I_2 = 0$ in proof of Theorem~5.1]{MR3787384}
As in Lemma~\ref{lem:Duality} we look at the solution
$(\bm \psi,r)\in {\bm H}^1_0(\Omega)\times L^2_0(\Omega)$ of
the Stokes problem with modified source $\bm{E}_h \bm w_h - {\bm S}(\bm w_h)$ and unit viscosity:
\begin{alignat*}{2}
(\nabla \bm \psi,\nabla \bm z) - (\nabla \cdot \bm z,r) & = (\bm{E}_h  \bm w_h -{\bm S}(\bm w_h),\bm z)\qquad &&\forall \bm z\in {\bm H}^1_0(\Omega),\\
(\nabla \cdot \bm \psi,q) & = 0\qquad &&\forall q\in L_0^2(\Omega).
\end{alignat*}
Testing the first equation with $\bm z = \bm{E}_h \bm w_h - {\bm S}(\bm w_h) \in {\bm H}^1_0(\Omega)$ and using \eqref{eqn:DefSCR}, \eqref{eqn:EhProperties}
and \eqref{eqn:CRApprox} leads to
\begin{align*}
\|\bm{E}_h \bm w_h - {\bm S}(\bm w_h)\|_{L^2(\Omega)}^2
&= (\nabla \bm \psi,\nabla (\bm{E}_h \bm w_h-{\bm S}(\bm w_h)))\\
&= (\nabla (\bm \psi- {\bm I}_{\text{CR}} \bm \psi),\nabla_h (\bm{E}_h \bm w_h - \bm w_h))\\
&\leq C_\text{CR} (1+C_{E_h}) h^s \| \bm \psi \|_{H^{1+s}(\Omega)} \| \nabla_h \bm w_h \|_{L^2(\Omega)}.
\end{align*}
The elliptic regularity assumption implies
$\| \bm \psi \|_{H^{1+s}(\Omega)} \le C_{\mathrm{ell},s} \|\bm{E}_h \bm w_h-{\bm S}(\bm w_h)\|_{L^2(\Omega)}$ and yields
\begin{align*}
  \|\bm{E}_h \bm w_h - {\bm S}(\bm w_h)\|_{L^2(\Omega)}
  \leq C_\text{CR} (1+C_{E_h}) C_{\mathrm{ell},s} h^s \| \nabla_h \bm w_h \|_{L^2(\Omega)}.
\end{align*}
Finally, a triangle inequality gives
\begin{align*}
  \|\bm w_h - {\bm S}(\bm w_h)\|_{L^2(\Omega)} & =
  \|\bm{E}_h \bm w_h - \bm w_h\|_{L^2(\Omega)}
  +\|\bm{E}_h \bm w_h - {\bm S}(\bm w_h)\|_{L^2(\Omega)}\\
  & \leq (C_{E_h} h + C_\text{CR} (1+C_{E_h}) C_{\mathrm{ell},s} h^s) \| \nabla_h \bm w_h \|_{L^2(\Omega)}.
\end{align*}
This concludes the proof.
\end{proof}

The previous result and similar arguments as in the conforming case enable us to prove the following theorem.

\begin{theorem}\label{thm:ConsistencyErrorEstimateCR}
Suppose that the Stokes problem satisfies Assumption~\ref{asmptn:regularity} and that the reconstruction
operator ${\bm I}_h$ satisfies Assumption~\ref{asmptn:reconstructionoperator}.
Then there holds
\begin{align*}
\|\nabla_h ({\bm v}_h-{\bm S}_h({\bm v}))\|_{L^2(\Omega)} \le (C_1 h +C_3 h^s) \|\mathbb{P} (-\Delta {\bm v})\|_{L^2(\Omega)},
\end{align*}
with $C_1>0$ and $C_3>0$ given by \eqref{eqn:IhApprox} and \eqref{eqn:DualityEstimateCR}, respectively.
Without Assumption~\ref{asmptn:reconstructionoperator}, a result similar to Theorem~\ref{thm:ConsistencyErrorEstimate_classical} is valid, i.e.,
\begin{align*}
\|\nabla_h ({\bm v}_h-{\bm S}_h({\bm v}))\|_{L^2(\Omega)} \le C_3  h^s \| \mathbb{P} (-\Delta {\bm v}) \|_{L^2(\Omega)}
+ \frac{C_{E_h} h}{\nu} \|\bm f - \mathbb{P} (\bm f) \|_{L^2(\Omega)}.
\end{align*}
\end{theorem}
\begin{proof}
The proof of the first result is nearly identical to the proof of Theorem~\ref{thm:ConsistencyErrorEstimate}
with slight changes concerning the application of $\nabla_h$ and the
replacement of Lemma~\ref{lem:Duality} by Lemma~\ref{lem:DualityCR}.
Likewise, the proof of the second result is almost identical to the proof of Theorem~\ref{thm:ConsistencyErrorEstimate_classical}. However, one term has to be estimated differently, as follows.
With $\mathrm{div} {\bm S}({\bm e}_h)=0$ and \eqref{eqn:EhProperties}, it holds
\begin{align*}
  \frac{1}{\nu}(\nabla \alpha, {\bm e}_h - {\bm S}({\bm e}_h))
  = \frac{1}{\nu}(\nabla \alpha, {\bm e}_h - \bm{E}({\bm e}_h))
  & \leq \frac{1}{\nu}\| \nabla \alpha \|_{L^2(\Omega)} \| {\bm e}_h - \bm{E}({\bm e}_h) \|_{L^2(\Omega)}\\
 & \leq \frac{C_{E_h} h}{\nu} \| \bm{f} - \mathbb{P} (\bm f) \|_{L^2(\Omega)}.
\end{align*}

\end{proof}

\section{Numerical Example}\label{sec:Numerics}

\begin{table}
\centering
\begin{tabular}{ccccc}
ndof & $ \!\| \nabla({\bm v} - {\bm v}_h) \|_{L^2(\Omega)} \!$ & order & $ \!\| \nabla({\bm v}_h - {\bm S}_h {\bm v}) \|_{L^2(\Omega)} \!$ & order\\
\hline
    379 & 1.4151e+00 &  -     & 5.0351e-02 & -     \\
   1414 & 9.7300e-01 &  0.542 & 3.0576e-02 & 0.722 \\
   5458 & 6.7235e-01 &  0.535 & 1.6366e-02 & 0.905 \\
  21442 & 4.6297e-01 &  0.540 & 8.4114e-03 & 0.964 \\
  84994 & 3.1819e-01 &  0.543 & 4.2567e-03 & 0.986 \\
 338434 & 2.1844e-01 &  0.545 & 2.1402e-03 & 0.996 \\
1350658 & 1.4988e-01 &  0.546 & 1.0729e-03 & 1.000 \\
\end{tabular}
\caption{\label{tab:ex1brlr_1c}Errors for the classical Bernardi--Raugel finite element method for $\bm f = \nabla (\sin(xy\pi))$ and $\nu = 1$.}
\end{table}

\begin{table}
\centering
\begin{tabular}{cccc}
ndof & $ \!\| \nabla({\bm v} - {\bm v}_h) \|_{L^2(\Omega)} \!$ & order & $ \!\| \nabla({\bm v}_h - {\bm S}_h {\bm v}) \|_{L^2(\Omega)} \!$\\
\hline
    379 & 1.4142e+00 &  -     & 8.5800e-11\\
   1414 & 9.7261e-01 &  0.542 & 1.2467e-13\\
   5458 & 6.7218e-01 &  0.535 & 1.9887e-14\\
  21442 & 4.6290e-01 &  0.540 & 4.3878e-14\\
  84994 & 3.1816e-01 &  0.543 & 9.8787e-14\\
 338434 & 2.1844e-01 &  0.545 & 2.2136e-13\\
1350658 & 1.4988e-01 &  0.546 & 4.4909e-13\\
\end{tabular}
\caption{\label{tab:ex1brlr_1m}Errors for the modified Bernardi--Raugel finite element method for $\bm f = \nabla (\sin(xy\pi))$ and $\nu = 1$.}
\end{table}

\begin{table}
\centering
\begin{tabular}{ccccc}
ndof & $\! \| \nabla({\bm v} - {\bm v}_h) \|_{L^2(\Omega)} \!$ & order & $ \!\| \nabla({\bm v}_h - {\bm S}_h {\bm v}) \|_{L^2(\Omega)} \!$ & order\\
\hline
    379 & 4.5794e+00 &  -     & 5.0351e+00 & -     \\
   1414 & 2.8168e+00 &  0.704 & 3.0576e+00 & 0.722 \\
   5458 & 1.5663e+00 &  0.850 & 1.6366e+00 & 0.905 \\
  21442 & 8.6350e-01 &  0.862 & 8.4114e-01 & 0.964 \\
  84994 & 4.8756e-01 &  0.828 & 4.2567e-01 & 0.986 \\
 338434 & 2.8682e-01 &  0.768 & 2.1402e-01 & 0.996 \\
1350658 & 1.7650e-01 &  0.703 & 1.0729e-01 & 1.000 \\
\end{tabular}
\caption{\label{tab:ex1brlr_2c}Errors for the classical Bernardi--Raugel finite element method for $\bm f = \nabla (\sin(xy\pi))$ and $\nu = 10^{-2}$.}
\end{table}

\begin{table}
\centering
\begin{tabular}{cccc}
ndof & $ \!\| \nabla({\bm v} - {\bm v}_h) \|_{L^2(\Omega)} \!$ & order & $ \!\| \nabla({\bm v}_h - {\bm S}_h {\bm v}) \|_{L^2(\Omega)} \!$\\
\hline
    379 & 1.4142e+00 &  -     & 8.5800e-09\\
   1414 & 9.7261e-01 &  0.542 & 1.2516e-11\\
   5458 & 6.7218e-01 &  0.535 & 6.5365e-13\\
  21442 & 4.6290e-01 &  0.540 & 1.3425e-12\\
  84994 & 3.1816e-01 &  0.543 & 2.7291e-12\\
 338434 & 2.1844e-01 &  0.545 & 5.5018e-12\\
1350658 & 1.4988e-01 &  0.546 & 1.1034e-11\\
\end{tabular}
\caption{\label{tab:ex1brlr_2m}Errors for the modified Bernardi--Raugel finite element method for $\bm f = \nabla (\sin(xy\pi))$ and $\nu = 10^{-2}$.}
\end{table}

\begin{table}
\centering
\begin{tabular}{ccccc}
ndof & $\! \| \nabla({\bm v} - {\bm v}_h) \|_{L^2(\Omega)} \!$ & order & $ \!\| \nabla({\bm v}_h - {\bm S}_h {\bm v}) \|_{L^2(\Omega)} \!$ & order\\
\hline
    379 & 4.3469e+02 &  -     & 5.0351e+02 & -     \\
   1414 & 2.6424e+02 &  0.721 & 3.0576e+02 & 0.722 \\
   5458 & 1.4134e+02 &  0.906 & 1.6366e+02 & 0.905 \\
  21442 & 7.2822e+01 &  0.960 & 8.4114e+01 & 0.964 \\
  84994 & 3.6908e+01 &  0.984 & 4.2567e+01 & 0.986 \\
 338434 & 1.8571e+01 &  0.995 & 2.1402e+01 & 0.996 \\
1350658 & 9.3137e+00 &  0.999 & 1.0729e+01 & 1.000 \\
\end{tabular}
\caption{\label{tab:ex1brlr_3c}Errors for the classical Bernardi--Raugel finite element method for $\bm f = \nabla (\sin(xy\pi))$ and $\nu = 10^{-4}$.}
\end{table}

\begin{table}
\centering
\begin{tabular}{cccc}
ndof & $ \!\| \nabla({\bm v} - {\bm v}_h) \|_{L^2(\Omega)} \!$ & order & $ \!\| \nabla({\bm v}_h - {\bm S}_h {\bm v}) \|_{L^2(\Omega)} \!$\\
\hline
    379 & 1.4142e+00 &  -     & 7.9830e-07\\
   1414 & 9.7261e-01 &  0.542 & 1.1388e-09\\
   5458 & 6.7218e-01 &  0.535 & 6.5440e-11\\
  21442 & 4.6290e-01 &  0.540 & 1.3407e-10\\
  84994 & 3.1816e-01 &  0.543 & 2.7257e-10\\
 338434 & 2.1844e-01 &  0.545 & 5.4872e-10\\
1350658 & 1.4988e-01 &  0.546 & 1.0965e-09\\
\end{tabular}
\caption{\label{tab:ex1brlr_3m}Errors for the modified Bernardi--Raugel finite element method for $\bm f = \nabla (\sin(xy\pi))$ and $\nu = 10^{-4}$.}
\end{table}

This sections gives a short numerical example to illustrate the theory.
We consider an L-shaped domain $\Omega:= (-1,1)^2 \setminus ((0,1) \times (-1,0))$ and the manufactured solution
\begin{align*}
 {\bm v}(r,\varphi)
& :=r^\gamma
\begin{pmatrix}
(\gamma+1)\sin(\varphi)\psi(\varphi) + \cos(\varphi)\psi^\prime(\varphi)
\\
-(\gamma+1)\cos(\varphi)\psi(\varphi) + \sin(\varphi)\psi^\prime(\varphi)
\end{pmatrix}^T,\\
 p_0(r,\varphi) &:= \nu r^{(\gamma-1)}((1+\gamma)^2 \psi^\prime(\varphi)+\psi^{\prime\prime\prime}(\varphi))/(1-\gamma)
\end{align*}
where
\begin{align*}
\psi(\varphi) &:=
\frac{1}{\gamma+1} \, \sin((\gamma+1)\varphi)\cos(\gamma\omega) - \cos((\gamma+1)\varphi)\\
&\qquad- \frac{1}{\gamma-1} \, \sin((\gamma-1)\varphi)\cos(\gamma\omega) + \cos((\gamma-1)\varphi)
\end{align*}
and $\gamma = 856399/1572864 \approx 0.54$, $\omega = 3\pi/2$
taken from \cite{Verfuerth1989}.
Note, that this yields $-\nu\Delta {\bm v} + \nabla p_0 = 0$.
To have a nonzero right-hand side we add $p_+ := \sin(xy\pi)$ to the pressure, i.e.\ $p := p_0 + p_+$ and $\bm f := \nabla(p_+)$.
Note that the exact solutions satisfy ${\bm v}\in {\bm H}^{1+s}(\Omega)$ and $p\in H^s(\Omega)$ for any $s<\gamma$.
Moreover, we set the viscosity parameter to either $\nu =1$, $\nu=10^{-2}$ or $\nu=10^{-4}$.

Tables~\ref{tab:ex1brlr_1c}-\ref{tab:ex1brlr_3m} compare the ${\bm H}^1$ errors of the classical
Bernardi--Raugel finite element method and its pressure-robust sibling on a series of unstructured uniformly red-refined meshes for $\nu = 1$ (Tables~\ref{tab:ex1brlr_1c} and \ref{tab:ex1brlr_1m}), $\nu=10^{-2}$ (Tables~\ref{tab:ex1brlr_2c} and \ref{tab:ex1brlr_2m})
and $\nu=10^{-4}$ (Tables~\ref{tab:ex1brlr_3c} and \ref{tab:ex1brlr_3m}).
For the classical method the distance between the discrete Stokes projector and the discrete solution is non-zero
and really scales with $\nu^{-1}$, but asymptotically converges with $h$
instead of $h^s$. At first glance this seems better than expected
in Theorem~\ref{thm:ConsistencyErrorEstimate_classical},
but the first term vanishes due to $\mathbb{P}(-\Delta {\bm v}) = \nu^{-1} \mathbb{P}(\nabla p_0) = 0$ in this example. This also pre-asymptotically leads to a slightly higher convergence order of the full error than in case of
$\nu=1$ at least for $\nu=10^{-2}$
and $\nu=10^{-4}$ where the $\mathcal{O}(h)$
error dominates at first.
The numbers of the modified pressure-robust variant convey that the discrete solution of the modified method and the discrete Stokes projector are identical
as predicted by Lemma~\ref{thm:ConsistencyErrorEstimate} (again due to $\mathbb{P}(-\Delta {\bm v}) = \nu^{-1} \mathbb{P}(\nabla p_0) = 0$). The numerical results
confirm that for pressure-robust methods, the
discrete velocity is independent of $\nu$.
However, this $\nu$-independence only holds up to
a quadrature error in the right-hand side,  which scales with
$\nu^{-1}$, and up to round-off errors.

\bibliographystyle{abbrv}


\section{Supplementary: Implementation Details}
\subsection{Generative Framework}
In this section, we describe the architecture that yields our reported results in the ZSL setting wherein there are no images from seen classes in test samples. For the SUN dataset, both the networks used for modeling mean and co-variance have linear (1800 and 2048 nodes), batch normalization and Relu layers. Additionally, the co-variance is restricted to be in the range of 0.5 to 1.5 for numerical stability via sigmoid activation. Both the networks are trained with ADAM optimizer \cite{kingma2014adam} using 0.001 and 0.1 as regularizer coefficients for means and covariance respectively.

For the AWA dataset, the generator networks have an architecture similar to SUN but consist of an additional dropout layer with probability 0.1. The parameters of the means are regularized with the coefficient $10^3$ while the parameters of covariance are regularized with the coefficient $10^4$.

For the CUB dataset, the generator networks have three linear layers (1200, 1800, 2048 nodes), 2 Relu, 2 batch normalization and 2 dropout layers. Their regularization coefficients are 0.01 and 0.1 for mean and covariance respectively. All the above networks are trained with a learning rate of $0.00001$. The training of these networks was additionally regularized via early stopping
\subsection{Adversarial Domain Adaption} As described above, ADA model comprises two discriminators ($D^{S,T}$), two classifiers ($C^{S,T}$) and two generators ($G^{S,T}$).

The discriminators $D^{S,T}$ are 5 layered neural networks comprising of 2 Linear layers of 1600 nodes and 1 node respectively, a single leaky Relu layer with a negative slope of 0.2 and batch normalization. The classifiers are single-layered networks with the number of nodes equal to the number of classes. $log(softmax(\mathbf{x}))$ is used as activation function for the classifiers $C^{S,T}$. The generators $G^{S,T}$ consist of three linear layers (1200,1200 and 2048 nodes), dropout layers, batch normalization and leaky Relu.

The overall objective is minimized using RMSprop (\cite{hinton2012neural}) optimizer with a learning rate of 0.00001. A manual seed of 100 has been used for all the ADA experiments.

\section{Supplementary: Training Procedures}
For brevity, we provide the training algorithms for our base ZSL model and the ADA ZSL model. We use the following loss definition for adversarial training as described in the paper
\begin{equation}
  \mathcal{L}= \mathcal{L}_{adv}^T+\mathcal{L}_{adv}^S+\chi \mathcal{L}_{cyc}+\xi\mathcal{L}_{clf}^T + \xi\mathcal{L}_{clf}^S
\end{equation}
where $\mathcal{L}_{adv}^{\{T,S\}} = \{L_G +L_D\}^{\{T,S\}}$ with
\begin{equation}
    L_G^{T} = \underset{c \sim {p_{c}}}{\mathbb{E}}[\beta \norm{G^T(x_{nc})-x_{nc}}_p - D^T_w \circ G^{T}(y_{nc})]

\end{equation}
\begin{equation}
    L_G^{S} = \underset{c \sim {p_{c}}}{\mathbb{E}}[\beta \norm{G^S(y_{nc})-y_{nc}}_p - D^S_w \circ G^{S}(x_{nc})]

\end{equation}
\begin{equation}
L^{T}_D =  \mathbb{E}_{c \sim {p_{c}}}[D^{T}_w \circ G^{T}(y_{nc})] - \mathbb{E}_{c \sim {p_{c}}}[D^T_w({x_{nc}})]
\end{equation}
\begin{equation}
L^{S}_D =  \mathbb{E}_{c \sim {p_{c}}}[D^{S}_w \circ G^{S}(x_{nc})] - \mathbb{E}_{c \sim {p_{c}}}[D^S_w(y_{nc})]
\end{equation}
Here, $D_w$ is the Wasserstein loss~\cite{arjovsky2017wasserstein} and $c$ denotes the class label.
\begin{align}
& \mathcal{L}_{\text{cyc}}(G^T, G^S)  =  \mathbb{E}_{c\sim p_{\text{c}}}[\norm{G^{S}\circ G^T(y_{nc})- x_{nc}}_p] \nonumber \\  &
+  \mathbb{E}_{c\sim p_{\text{c}}}[\norm{G^{T}\circ G^S(x_{nc})- y_{nc}}_p].\lbleq{cycle}
\end{align}
Here, $||\cdot||_p$ denotes the $L_p$ norm.
\begin{align}
L^{T}_{clf} = \mathbb{E}_{c \sim {p_{c}}}[L(C^{T}_{clf}\circ G^{T}(y_{nc}),Y^{T})] \nonumber \\
+ \mathbb{E}_{c \sim {p_{c}}}[L(C^{T}_{clf}(x_{nc}),\bar{Y}^U)]
\end{align}
\begin{align}
L^{S}_{clf} = \mathbb{E}_{{c} \sim {p_{c}}}[L(C^{S}_{clf}\circ G^{S}(x_{nc}),Y^{S})] \nonumber \\
+ \mathbb{E}_{c \sim {p_{c}}}[L(C^{S}_{clf}(y_{n c}),\bar{Y}^U)]
\end{align}

\begin{algorithm}[H]
	\caption{\textsc{Superficial Training Scheme} \label{alg:system}}
	\begin{algorithmic}[1]
		\State Train ResNet101 on Imagenet
		\State Randomly initialize model parameters
		\State Initialize dataset $\mathcal{D}_{ada}\gets \emptyset $
			\State Train the generative model to describe data from class $c$ by $p(\mathbf{x}|c,\Theta)$ $\forall c\in[1,...C]$
			\State $Source \longleftarrow$ data sampled from generative model
			\State $Target \longleftarrow$ unlabelled test data
			\State Initialize weights of $G^T:S\longrightarrow T$,$G^S:T\longrightarrow S$, $D^S$, $D^T$
			\State Augment class labels to $G^T$ and $G^S$
		\For{epoch $i=1:K$}
		\State Randomly sample $o_i^{s}\sim Source$, $o_i^{t}\sim {Target}$
		\State Perform ADA
		\EndFor

	\end{algorithmic}
\end{algorithm}
\begin{algorithm}[H]
	\begin{algorithmic}[1]
		\State {\bfseries Input:}Maximum iterations $N_{step}$, batch size $n$, the iteration number of discriminator in a loop $n_{d}$, RMSprop hyperparameters $\alpha$, class attributes $\{a_c\}_{c=1}^{S+U}$, ADAM hyperparameters $\alpha_2,\beta_1,\beta_2$

		\State // \textit{Pre-training}
		\For{$c = 1...S+U$}
		\State $\mathbf{\zeta_c} = f_{\mathbf{\Theta}}(\mathbf{a_c})$
         \EndFor
         \For{iter $= 1,..., N_{step}$}
            \State Sample minibatch of examples $x_1,x_2,...x_n$

         \State
         $L =\smashoperator[r]{\sum\limits_{c=1}^S}\smashoperator[r]{\sum\limits_{n:y_n=c}}\left(x_n - f_{\mathbf{\mu}}(a_c)\right)^T\left[f_{\Sigma}(a_c) \right] \left(x_n - f_{\mathbf{\mu}}(a_c)\right)$

          \State  $\mathbf{\Theta} \gets \text{Adam}(\bigtriangledown_{\mathbf{\Theta}} L,\mathbf{\Theta}, \alpha_2, \beta_1, \beta_2)$
         \EndFor
    \State // \textit{Adversarial Domain Adaptation}
	\State Initialize $\{G,D,C\}^{T,S}$
	\State $\lambda_d = Params(D^T,D^S,C^T,C^S)$
	\State $\lambda_g = Params(G^T,G^S)$

        \For{iter $= 1,..., N_{step}$}
            \State Sample minibatch $x_1,x_2,...x_n$ from test data for training $G^T$
            \State Sample minibatch $x'_1,x'_2,...x'_n$ from class conditional distributions for training $G^S$
        	\State Compute the overall loss $\mathcal{L}$ using Eq.1
        \State $\lambda_g \gets \text{RMSprop}(\bigtriangledown_{\lambda_g} \mathcal{L},\lambda_g, \alpha)$
		\For{$t = 1$, ..., $n_{d}$}
		 \State Sample minibatch of examples $x$
        \State Compute the overall loss $\mathcal{L}$ using Eq.1
        \State $\zeta_d \gets \text{RMSprop}(\bigtriangledown_{\lambda_d} \mathcal{L},\lambda_d, \alpha)$
        \EndFor

		\EndFor
	\end{algorithmic}
	\caption{All the notations have same meaning as that in the running text. \textit{Params}() return the parameters of the model in (.)}
	\label{alg_label}

\end{algorithm}

\section{Introduction}
In the conventional image classification tasks, examples from all classes are available during the training of the model. This assumption rarely holds in real-world problems, where we do not have the corresponding ubiquity of representative images from each class. Also, it is common knowledge that humans do not require prior visual evidence of a category to recognize an example from that category. Given that a child sees a picture of a horse and reads a description about zebra's appearance, he/she would more likely than not be able to easily recognize a zebra when an image is shown. The zero-shot learning (ZSL) problem~\cite{cmt,xian2018zero} in machine learning is motivated by similar considerations and seeks to exploit the existence of a labeled training set of `seen' classes and the knowledge about how each `unseen' class relates semantically to the seen classes.

The success of ZSL lies in learning an effective semantic representation (e.g. attributes / textual features) for the successful transfer of knowledge from the seen to the unseen classes. In Sec.~\ref{sec:relwork}, we provide a detailed overview of the prior work on ZSL, but in particular generative ZSL methods~\cite{xian2018feature,wang2017zero,verma2017simple,vermageneralized} have become quite popular recently, by the virtue of their ability to \emph{generate} labeled examples for the unseen classes. However, a key requirement in such methods is the reliable estimation of the class distribution of seen and unseen classes. Even then, zero-shot learning suffers from hubness problem \cite{dem} mostly because of the use of nearest neighbor classifiers exploiting different distance metrics. It can be mitigated by using neural nets or any classifier which does not explicitly compare the inter-class distances in high dimensional data for label prediction. Hence, a generative model makes it plausible to train deep classifiers on synthesized data from the unseen classes.

\begin{figure*}[!htbp]
    \centering
    \includegraphics[width=0.82\textwidth]{Overall.png}
    \caption{The overall architecture of the proposed approach. All the notations are consistent with that described in Section 2. $\mathbf{A_c}$ denotes the class attributes for all classes i.e. $\{\mathbf{a_c}\}_{c=1}^{S+U}$. }
    \label{fig:architecture}
\end{figure*}

 A simple, yet principled, way to construct generative models for ZSL is to learn the class distributions for the seen and unseen classes~\cite{verma2017simple,wang2017zero}. While this is straightforward for seen classes (for which we have access to labeled data), it can't be done for the unseen classes. In recent work, ~\cite{verma2017simple} used exponential family to model the distributions of the class conditionals in terms of learnable parameters. This is an effective model; however, their approach does not extend to non-exponential family distributions. Moreover, they used offline learning techniques to learn the parameters of seen classes, and rely on kernel-based regression to estimate the class parameters, given the class attributes. The model also requires careful tuning of hyperparameters. Our model, on the other hand, exploits the advantages of neural nets and end-to-end training to provide stability during the learning phase and remains less susceptible to hyperparameter variations.

However, such a model alone is not sufficient as there may be a domain shift between the original unseen class data and the synthesized unseen class data. The presence of acute domain shift between the seen and unseen classes hinders the performance of ZSL models~\cite{UDA}. Since the predictions for the unseen classes rely on the transfer of knowledge learnt from the seen classes, we might have poor performance on unseen classes due to the domain shift. We note that by enforcing domain adaptation to tackle problem settings where the train and test distributions are far apart, the model's performance can be greatly improved. An earlier approach~\cite{UDA} used the idea of \emph{joint} sparse coding for minimizing domain shift between the seen and unseen class data, however, since then there have been developments in adversarial domain adaptation that enable robust detection and resolution of domain shift~\cite{adda,cycada}. Adversarial learning and adaptation methods have found applicability in a wide range of fields from robotics navigation~\cite{icra} to recommender systems~\cite{recgan,irgan}. Several adversarial adaptation techniques like ADDA\cite{adda} require explicit source-target pairs of data points. Such a luxury is not present in Zero-Shot transductive setting where the test data is unlabelled. Similarly, unsupervised domain matching methods like CycleGAN\cite{cyclegan} use cyclic consistency to find the data point most similar to the source sample and then minimize the gap between these two. Though this is effective in maintaining the inherent clusters, it can match the unrelated class clusters together in the source and target domain if the classes are close enough.

Motivated by these desiderata, in this work, we develop an Adversarial Domain Adaptation framework for ZSL that leverages a generative ZSL model to improve upon the classification for unseen classes. Our model can transform the synthesized samples for unseen classes into the true test/unseen class domain while maintaining the data clusters associations. We first learn a generative model for the class conditional distribution of the seen and unseen classes by utilizing labeled data from the seen classes. Then, by domain adaptation, we explicitly bring closer the learnt distribution and the true distribution of the unseen class conditionals. We employ a scheme of cyclically consistent adversarial domain adaptation~\cite{cycada} to minimize domain shift without assuming any particular parametric form of the source and target distributions.

To the best of our knowledge, there is no adversarial framework for semi-supervised domain matching where explicit pairs of data points are not given but an external agent associates noisy labels to the samples. In addition, since we leverage neural nets for classification, we overcome the hubness issue, by virtue of not classifying based on explicit distances (unlike classical KNN type algorithms).

\section{The Proposed Approach}
\label{sec:model}

Our approach consists of a pre-training phase followed by adversarial domain adaptation (ADA). We first describe the generative model and then elaborate on the ADA setting. A detailed illustration of our method is shown in Figure 1.

\subsection{The Generative Model}
We model the data from class $c$ by a class conditional probability $p(\mathbf{x}|c,\mathbf{\Theta})$ where $\mathbf{\Theta}$ denotes the global parameters of the model. We do not have any restriction on the the type of distribution chosen. Let us denote the total number of classes whose examples are seen during training by $S$, and the classes, none of whose examples are seen during training by $U$. For the sake of defining the prediction rule formally, assume the unseen classes are known. Then, for an observation $\mathbf{x}$ from either a seen or unseen class $c$, where $c \in [1, S+U]$, we have $y_n = c$, and, assume the input to be generated as
$
 x_n \sim p(\mathbf{x}|c,\mathbf{\Theta})

$

Under this framework, given test example $\mathbf{x_+}$, the predicted class $\hat{y}_+$ can be given by computing the most-probable class as follows
$
    \hat{y}_{+} = \argmax_c p(c|\mathbf{x_+},\mathbf{\Theta})
$ and using Baye's Rule we have,
\begin{equation}
    p(c|\mathbf{x_+},\mathbf{\Theta})  = \frac{p(\mathbf{x_+}|c,\mathbf{\Theta}) p(c|\mathbf{\Theta})}{\sum_{c \in [1,S+U]}{p(\mathbf{x_+}|c,\mathbf{\Theta}) p(c|\mathbf{\Theta}) }}
\end{equation}
Thus
\begin{equation}
    \hat{y}_+ = \underset{c}{\argmax}\  p(\mathbf{x_+}|c,\mathbf{\Theta})p(c|\mathbf{\Theta})
\end{equation}

For the sake of simplicity, we ignore the estimation of class prior probabilities and choose to treat them as equal for all the classes. However, their correct estimation can in principle provide better results. The prediction rule then becomes:
\begin{equation}
    \hat{y}_+ = \underset{\mathbf{c}}{\argmax}\  p(\mathbf{x_+}|c,\mathbf{\Theta})
\end{equation}
If labeled training data for all the classes are available, then standard inference techniques like Maximum Likelihood Estimation (MLE), Maximum-a-Posteriori (MAP) Estimation, or fully Bayesian inference can be used to determine the class conditional distributions. However, since the unseen classes do not have labeled training examples, we need a way to ``extrapolate'' the seen class distribution parameters to unseen class distribution parameters. This will be done via the class attribute vectors as described ahead (\ref{sec:attmapping}).

First, assuming $\mathbf{X}$ and $\mathbf{C}$ ($c_k \in S\cup U$) denote the inputs and the associated output class labels respectively, a standard generative approach seeks to  maximize their joint distribution $\mathbb{P}[\mathbf{X^{S\cup U},\mathbf{C^{S\cup U}}}|\mathbf{\Theta}]$.

Assuming i.i.d. observations, we have
\begin{equation}
    \mathbb{P}[\mathbf{X^{S\cup U},\mathbf{C^{S\cup U}}}|\mathbf{\Theta}] = \mathbb{P}[\mathbf{X}^S,\mathbf{C^{S}}|\mathbf{\Theta}] \mathbb{P}[\mathbf{X^U},\mathbf{C^{U}}|\mathbf{\Theta}]
\end{equation}
Since, $\mathbf{X}^U,\mathbf{C^{U}}$ are unavailable for $\Theta$ estimation,
usually $\mathbb{P}[\mathbf{X^{S},\mathbf{C^{S}}}|\mathbf{\Theta}]$ is maximized instead, expecting the learnt $\Theta$ to behave as a proxy to true value.
\begin{equation}
    \mathbb{P}[\mathbf{X^{S},\mathbf{C^{S}}}|\mathbf{\Theta}]=\prod_{\mathbf{x},c \sim S} p(\mathbf{x},c|\mathbf{\Theta})
\end{equation}
\begin{align}
   \Rightarrow log(\mathbb{P}[\mathbf{X^{S},\mathbf{C^{S}}}|\mathbf{\Theta}])=\underset{\mathbf{x},c \sim S}{\sum}log(p(\mathbf{x},c|\mathbf{\Theta})) \nonumber \\
   =\underset{\mathbf{x},c \sim S}{\sum}log(p(\mathbf{x}|c,\mathbf{\Theta}))+log(p(c|\mathbf{\Theta}))
\end{align}
Since we are not modelling the class probability distribution $p(c|\mathbf{\Theta})$, the objective becomes
\begin{equation}
    \underset{\mathbf{\Theta}}{\argmax} \underset{\mathbf{x},c \sim S}{\mathbb{E}}[log(p(\mathbf{x}|c,\mathbf{\Theta}))]
\end{equation}

This sub-optimal $\Theta$ produces an inherent domain shift between the true unseen class distribution and the learnt distribution. We mitigate this by using adversarial domain adaptation to bring the unseen distribution and learnt distributions closer (refer \ref{sec:ADA}).

\subsubsection{Mapping Class Attributes to Class Parameters}
\label{sec:attmapping}
Since each class is described in terms of attribute vectors $\mathbf{a}_c$, we condition the class distribution on their respective attribute vector $\mathbf{a}_c$. Let these class-specific parameters be $\zeta_c$ which can be uniquely determined from the class attribute vector $\mathbf{a_c}$ and global parameters $\mathbf{\Theta}$ by a functional mapping $f$. This mapping for most purposes will be a complicated relationship and using a linear mapping (e.g., as done in~\cite{verma2017simple}) here would severely affect the generation quality of the network. We model this function $f:\{\mathbf{a}_c\} \rightarrow \{\mathbf{\zeta}_c\}$ using neural networks with trainable weights $\mathbf{\Theta}$  bringing extensive expressiveness and hierarchical relationships among attribute features. Thus, the class parameters can be written as
\begin{equation}
    \mathbf{\zeta_c} = f_{\Theta}(\mathbf{a_c})
\end{equation}
\begin{figure}[H]
\centering
\begin{minipage}[c]{0.45\textwidth}
\centering
    \includegraphics[width=\textwidth,height=0.45\textwidth]{PGM.png}
    \caption{Samples of seen/unseen classes $\textbf{x}_n$ generated by the class conditional distribution defined by parameters $\zeta_c$  which in turn are the outputs of neural network $f_\Theta$ with $\textbf{a}_c$ as input}
    \label{fig:PGM}
\end{minipage}
\end{figure}

Such an approach leads to a stable training procedure w.r.t hyperparameters and enables us to perform the joint learning of $f_{\boldsymbol{\Theta}}$ and consequently, class parameters $\{\zeta_c\}$ in an end to end fashion. This is an important difference between our approach and the generative approach used by~\cite{verma2017simple}. Their approach first learns the class conditional parameters and then learns the attribute to class parameter mappings. We provide empirical justification for the stability of our approach in the Results section.

For simplicity, we take $p(\mathbf{x}|c,\mathbf{\Theta})$ to be a Gaussian distribution with parameters mean and co-variance, $\zeta_c =  \{\mathbf{\mu}_c,\Sigma_c \}$ where $c \in S$.
We model $\mathbf{\mu_c}$ and $\mathbf{\Sigma}^{-1}_c$ as non-linear functions of the attribute vector $\mathbf{a_c}$ with neural networks of weights $\mathbf{\Theta}=\{ \mathbf{\theta^{\mu},\theta^{\Sigma}} \}$ in the following manner,
\begin{equation*}
    \mathbf{\mu_c} = f_{\theta^{\mathbf{\mu}}}(\mathbf{a_c}),\
        \mathbf{\Sigma_c}^{-1} = diag(f_{\theta^{\mathbf{\Sigma}}}(\mathbf{a_c})),\      \mathbf{x} \sim N \left(\mu_c, \mathbf{\Sigma_c}\right)
\end{equation*}
To ensure the condition of the covariance matrix ($\Sigma_c$) being a positive semi-definite matrix we model the inverse covariance to a diagonal matrix with positive diagonal entries. Thus $f_{\theta^{\Sigma}}$ outputs a vector in $\mathbb{R}^d_{>0}$ where d is the dimension of mean vector (equivalently the dimension of semantic space). The overall objective function becomes:
\begin{align}
\small
\underset{\mathbf{\theta^\mu,\theta^\Sigma}}{\argmax} \underset{(\mathbf{x},c) \sim S}{\mathbb{E}}\left[log(\mathbf{\Sigma_c}^{-1} )
    -\left(\mathbf{x} - \mathbf{\mu_c}\right)^T  \mathbf{\Sigma_c}^{-1}\left(\mathbf{x} - \mathbf{\mu_c}\right)\right]
\end{align}

We again emphasize the fact that choosing Gaussian distribution is only for expositional purposes and one can also try other non-exponential family distributions as a part of inductive bias. The model does not restrict the choice of class conditional distribution to exponential family distribution
We take $\mathbf{x_n}$ as the features extracted from dataset images by resnet-101\cite{resnet} pre-trained on the Imagenet dataset \cite{imagenet2015}.  For the rest of the paper, $\{\mathbf{x}_n\}_c$ denotes the entire test data comprising of samples from all the classes. Similarly, $\{y_n\}_c$ denotes the samples generated from the generative model (as defined here) for all the classes. For the rest of the paper, we refer to the generative model defined above as the base ZSL model.

\subsection{Adversarial Domain Adaptation}
\label{sec:ADA}

The procedure described in the previous section only leverages the data from seen classes to estimate the class conditional distributions of all the classes. However, if there is a domain shift between the seen and unseen classes, then the estimated unseen class conditional distribution would also suffer from this domain shift due to reliance on the seen class data. Hence, to mitigate the issue of domain shift,  we propose to incorporate the unlabelled data from the unseen classes. In our overall architecture, we denote the process of learning the ZSL model parameters $\mathbf{\Theta}$ as `pre-training.' Based on the generative framework learned during pre-training, we can sample the class-conditional distribution for unseen classes to generate the unseen samples. We then minimize the domain gap between the generated distribution of the unseen classes and the actual distribution of the unseen classes.

In this section, we denote the source domain as $S$ and the target domain as $T$. Through adversarial adaptation, we aim to bring the target distribution of $\{y_n \}_c$ (referred as $ y_{nc}$) closer to the source distribution  of $\{x_n \}_c$ (referred as $x_{nc}$); hence we learn a function $G^T(y_{nc})$ that maps class conditionals from the generated distribution $y_n$ to the real test distribution $x_n$ for all unseen classes $\{c\}_{c=1}^{U}$. Hence, $G^T: S\rightarrow T$ is a mapping from source $S$ to target $T$. Similarly, we define another function $G^S: T\rightarrow S$ that maps the class conditionals from the real test distribution to the same latent space as the class conditionals from the generated distribution. $D^T$ and $D^S$ are the corresponding discriminators.

Our design is inspired by CycleGAN \cite{cyclegan} and we make modifications to its base architecture for supporting zero shot learning. We consider a cyclic consistency loss instead of the vanilla adversarial loss (and variants) primarily because we want to learn as constrained a latent space for the Generator as possible. Cyclic consistency is an additional constraint on top of the adversarial loss that acts as an appropriate regularizer for transfer learning, as motivated in the original paper~\cite{cyclegan} . We also justify the cyclic consistency loss empirically in the Ablation Study section 4.4

\subsubsection{Label Augmentation}
Inspired by conditional GAN\cite{conGAN}, we augment the input to the generators $G^T$ and $G^S$ with the respective class labels, which facilitates the preservation of relationships between the synthesized data and their correct class labels. For $G^S$ the input data (test data) is unlabeled, hence we use the predictions from our pre-trained ZSL model as the guiding labels. Note that these labels are noisy labels and both the generators and discriminators should be capable of handling data corruption during the training phase. This is yet again a problem with the conventional GAN architectures.
\subsubsection{Classifiers}
 We handled label recovery by adding two classifiers ($C^T,C^S$) in parallel with the discriminators. The parameters of the classifiers are trained jointly with the corresponding discriminators. Recently, parallel to our work, \cite{robustGAN} gave theoretical support to the use of external classifier with conditional GAN architecture to counter noisy data labels. The classifiers provide an additional benefit of enforcing linear separability for the generator but the impact is reduced if the classifier is multilayered. We justify clustering in the next section.

\subsubsection{Optimization Function}
Let the loss defined in CycleGAN which consists of cyclic consistency loss ($\mathcal{L}_{cyc}$) and the adversarial loss ($\mathcal{L}_{adv}$) for domains T and S be $\mathcal{L}$
\begin{equation}
  \mathcal{L}= \mathcal{L}_{cyc}+\mathcal{L}_{adv}^T+\mathcal{L}_{adv}^S
\end{equation}
For our case, $L_1$ norm worked the best for cyclic consistency loss\cite{cyclegan} $\mathcal{L}_{cyc}$, while Wasserstein loss \cite{arjovsky2017wasserstein} was found suitable for $\mathcal{L}_{adv}$. Additionally, identity regularizer with $l1$ norm (see eq 12) was added to the generator to ensure that the output domains for each generator remain unmodified if given as input. We add the classification loss ($\mathcal{L}_{clf}$) of the real data (not generated by G) to the discriminator loss during adversarial training. We ensured that the classification loss is not added at the beginning for data transformed by generators in accordance with mismatch loss addition for only real images~\cite{conGAN}. We evaluated cross entropy loss for both the correct and mismatched pairs of label-image. This enforces a stronger clustering than considering the loss term for only mismatched pairs, as in \cite{conGAN}. However once the GAN training has converged and the classification accuracy over the generated data becomes close or greater than the accuracy of pseudo-labels, we do a corruption recovery by training the classifiers over both the transformed samples $G^T(y_{nc})$ and true data samples $x_{nc}$ (\textit{refer eq (\ref{eq: classfifier})})

With $\chi, \xi,\beta$ as tune-able hyper parameters, the overall loss function then becomes
\begin{equation}
  \mathcal{L}= \mathcal{L}_{adv}^T+\mathcal{L}_{adv}^S+\chi \mathcal{L}_{cyc}+\xi\mathcal{L}_{clf}^T + \xi\mathcal{L}_{clf}^S
\end{equation}
where $\mathcal{L}_{adv}^{\{T,S\}} = \{L_G +L_D\}^{\{T,S\}}$ with
\begin{equation}
    L_G^{T} = \underset{c \sim {p_{c}}}{\mathbb{E}}[\beta \norm{G^T(x_{nc})-x_{nc}}_p - D^T_w \circ G^{T}(y_{nc})]
\end{equation}
\begin{equation}
L^{T}_D =  \underset{c \sim {p_{c}}}{\mathbb{E}}[D^{T}_w \circ G^{T}(y_{nc})] - \underset{c \sim {p_{c}}}{\mathbb{E}}[D^T_w(x_{nc})]
\end{equation}
Here, $D_w$ is the Wasserstein loss~\cite{arjovsky2017wasserstein} and $c$ denotes the class label.
\begin{align}
& \mathcal{L}_{\text{cyc}}(G^T, G^S)  =  \mathbb{E}_{c\sim p_{\text{c}}}[\norm{G^{S}\circ G^T(y_{nc})- x_{nc}}_p] \nonumber \\  &
+  \mathbb{E}_{c\sim p_{\text{c}}}[\norm{G^{T}\circ G^S(x_{nc})- y_{nc}}_p].\lbleq{cycle}
\end{align}
Here, $||\cdot||_p$ denotes the $L_p$ norm.
\begin{align}
L^{T}_{clf} = \mathbb{E}_{c \sim {p_{c}}}[L(C^{T}_{clf}\circ G^{T}(y_{nc}),Y^{T})] \nonumber \\
+ \mathbb{E}_{c \sim {p_{c}}}[L(C^{T}_{clf}(x_{nc}),\bar{Y}^U)]
\label{eq: classfifier}
\end{align}
Similarly, we can define $L^S_{clf} , L^S_D , L^S_G$. Please refer to supplementary  material for exact equations and training algorithm.

\section{Related Work}
\label{sec:relwork}
Due to its ability to overcome the drawbacks of conventional classification problems, ZSL has attained tremendous recent interest for a wide range of AI problems, including those in computer vision. Earlier works~\cite{DAP,IAP} on ZSL were based on directly or indirectly mapping the instances of specific examples to their class-attributes. The learned mapping was then used during inference; this mapping works by first projecting the unseen data to the class-attribute space and then using the nearest neighbor search to classify the unseen image.
Another approach for ZSL focuses on learning the map of bi-linear compatibility between the visual space and the semantic space. ALE \cite{akata2013label}, DEVISE \cite{frome2013devise}, SJE \cite{SJE}, ESZSL \cite{ESZSL2015}, and SAE \cite{SAE2017} are based on the approach of measuring the bi-linear compatibility.

Generative models~\cite{mishra2017generative,Chen_2018_CVPR,xian2018feature,verma2017simple,guo2017synthesizing,BucherZSL,wang2017zero} have shown promising results for both ZSL and GZSL setups. Another advantage of the generative approach is that by using synthesized samples, we can convert the ZSL problem to the conventional supervised learning problem that can handle the biases towards the seen classes. The \cite{verma2017simple} used a simple generative model based on the exponential family framework while \cite{guo2017synthesizing} synthesize the classifier. While recent generative approaches for the ZSL are deep generative models based on the VAE \cite{VAE} and GAN \cite{GAN}. The approach \cite{vermageneralized,BucherZSL,mishra2017generative} is based on the VAE architecture while \cite{xian2018feature,Chen_2018_CVPR,lu2017zero} used the adversarial sample generation based on the class conditioned attribute.

In ZSL, the train and test classes are disjoint and hence there is a high probability of domain shift for the unseen classes. This is another challenge in the ZSL setup and needs to be handled. Previously, very few works have handled the domain shift problem and worked on both the transductive as well as inductive settings. \cite{verma2017simple} adapted to the new domain by simple Gaussian mixture model updates. \cite{song2018transductive} used the unbiased embedding in the transductive setting. \cite{kodirov2015unsupervisedDA,ye2017zero_DOMAIN} proposed unsupervised domain adaption for the ZSL. \cite{saligram2016learningJoint} used the structural SVM formulation for domain adaption.

In this paper, we propose the design of a deep generative model that has many differences as compared to the previously proposed VAE/GAN based deep generative models for ZSL. The VAE based architecture minimizes the ELBO \cite{VAE} by using a scheme of approximate optimization, making it less robust in handling domain shift. This also applies to the latent class distributions learned by VAE. While we explicitly estimate the class conditional distributions, VAE based methods learn these distributions as latent variables via approximate inference. Thus the complexity of our model in representing the class conditionals is on par with VAE based models but on the other hand, we reap the benefits of direct optimization. The GAN based generative approach is difficult to train, requiring a lot of seen class examples during training. Moreover, they need the attribute vectors of unseen classes at the beginning of the procedure while our model can handle on the fly addition of new classes.
To this end, we propose a simple CNN based architecture that can learn any parametric distribution with exact optimization, and unlike the GAN based approach, has stable training. This makes it especially suitable for domain shift minimization by adapting the distribution of the unseen classes.

\section{Experiments and Results}
\label{sec:model}
\begin{figure*}[t]
    \centering
    \begin{tabular}{ccc}
        \includegraphics[width=0.32\textwidth,height=0.30\textwidth]{pretrained.png} &
    \includegraphics[width=0.32\textwidth,height=0.30\textwidth]{adversarial_wgan_iterative.png} &
    \includegraphics[width=0.32\textwidth,height=0.30\textwidth]{Awa_HyperParam.png}\\
    (a) & (b) &(c)
    \end{tabular}

    \caption{\textbf{(a)} shows the t-SNE plot for the output of the generative model as compared to the test data. Crosses represent the test data while dots represent the generated data. The domain shift is visible in this plot. \textbf{(b)} shows the t-SNE plot after domain shift minimization with our model. The scale for the axes in \textbf{(a)} and \textbf{(b)} is kept constant for comparison. The model can allot the clusters correctly except those where the prediction of pseudo-labels suffered a lot and recovery was difficult. \textbf{(c)} shows the stability of the generative model wrt regularization coefficients on the AWA2 dataset. The x-labels and y-labels are the weight decay in Adam optimizer for learning the NN parameters predicting the Mean and Sigma of the class-conditional distributions respectively. The shaded grid values represent the top-1 accuracy obtained for the given configuration of hyperparameters. Note that even on a logarithmic scale the changes in accuracy are about 1-3\%}
    \label{fig:Clustering}
\end{figure*}
To demonstrate the effectiveness of the proposed approach we performed extensive experimentation on the three standard datasets for ZSL, namely AWA2 \cite{xian2018zero}, CUB-200 \cite{welinder2010caltech} and SUN \cite{xiao2010sun}. In all the experiments, we follow the newly proposed train test split suggested by \cite{xian2018zero}. Since we are using the pre-trained resnet-101 model, therefore, we first sought to make sure that any class that belongs to the test classes is not present in the ImageNet \cite{imagenet2015} training samples. This was already rectified in the split proposed by \cite{xian2018zero} for ZSL. For reference, the network architecture and training procedure is provided in the supplementary material.
\begin{table}[H]
\small
\centering
 \addtolength{\tabcolsep}{1.0pt}
 \begin{tabular}{||c|c|c|c||}
 \hline
 Dataset & Attribute/Dim & \#Image & Seen/Unseen Class \\ [0.5ex]
 \hline\hline
AWA2 & A/85 & 37322 & 40/10 \\
CUB & A/1024 & 11788 & 150/50 \\
SUN & A/102 & 14340 & 645/72 \\
 \hline
 \hline
 \end{tabular}
 \caption{Datasets used in our experiments, and their statistics} \label{tab:data}
 \vspace{-2em}
\end{table}
\textbf{Animal with Attribute (AWA2):} The dataset has 50 classes of animals, with 40 classes used for training data and the rest 10 as test data. Each class also has a human annotated 85-dimensional attribute vector associated with it.

\textbf{CUB-200:} This is a fine-grained dataset, containing 200 classes of birds, with 150/50 as the train/test ZSL class split. It has $11788$ data points and 1024-dimensional human-annotated attribute vectors for each class. The attribute vector comprises of 312-dimensional CUB vector appended with word vectors of class names as proposed by \cite{xian2018zero}.

\textbf{SUN Seen Recognition:} There are a total of 14340 images from 717 classes. Hence, every class has nearly 20 samples. Each class is associated with a 102-dimensional human-annotated attribute vector.

\subsection{Zero-Shot Learning (ZSL)}

We report per class accuracy as is the convention in standard ZSL. It is a better metric to report the accuracy of the model as compared to the overall (across classes) accuracy when the classes are unbalanced. We use the newly proposed splits \cite{xian2018zero} for dividing the train and test examples to ensure that the Imagenet classes (used for training ResNet) and the test classes are disjoint. We use the corresponding attribute vector provided against each dataset. Please refer to table-\ref{tab:data} for details on the dataset.

The results of the ZSL setting are shown in the table-\ref{tab:zsl}. As apparent from the table, the proposed approach shows a significant improvement over the previous state-of-art approaches. On the SUN dataset and the AWA2 dataset, we have our top-1 accuracies $63.3\% $ and $70.4\%$ respectively, which are better than its close competitor \cite{verma2017simple}. Also, their top-1 accuracy on the fine-grained CUB dataset is \textbf{significantly} reduced to $49.2\%$,  compared to our model's top-1 accuracy of $70.9\%$. Our model thus performs consistently well and beats other models on all the three benchmark datasets.

Additionally, our approach is more stable to hyperparameter variations as compared to the other competing generative approaches like GFZSL\cite{verma2017simple}. We get only 2-4\% drops in accuracy on a logarithmic scale, unlike GFZSL\cite{verma2017simple} (\textit{figure \ref{fig:Clustering},c}). ADA model has three hyperparameters $\xi,\chi,\beta$ and we chose $\chi=10,\beta =0.5\chi$  as used by the original Cycle GAN \cite{cyclegan}. Random search was used for $\xi=0.0001$.

\begin{table}[h!]
\small
 \centering
 \addtolength{\tabcolsep}{5pt}
 \begin{tabular}{|l| c |c | c | c |}
 \hline
  & {\textbf{SUN}} & {\textbf{CUB}} &  {\textbf{AWA2}} \\ \hline
 \textbf{Method} &  \textbf{PS} &   \textbf{PS} &  \textbf{PS} \\
  \hline
 \textbf{DAP}\cite{DAP} & 39.9  & 40.0    & 46.1 \\
 \textbf{IAP} \cite{IAP}  & 19.4  & 24.0    & 35.9\\
 \textbf{CONSE} \cite{norouzi2013zero} &  38.8 &  34.3   & 44.5 \\
 \textbf{CMT} \cite{cmt}  & 39.9  & 34.6 &  37.9 \\
  \textbf{SSE} \cite{saligram2016learningJoint}  & 51.5  & 43.9 & 61.1\\
  \textbf{LATEM} \cite{latem}  & 55.3  & 49.3    & 55.8\\

 \textbf{DEVISE} \cite{frome2013devise}  & 56.5 & 52.0 & 59.7 \\
 \textbf{SJE} \cite{SJE}  & 53.7  & 53.9 &  61.9 \\
 \textbf{ESZSL} \cite{ESZSL2015} & 54.5  & 53.9   & 58.6 \\
 \textbf{SYNC}\cite{changpinyo2016synthesized}  & 56.3  & 55.6  & 46.6\\
 \textbf{SAE} \cite{SAE2017}  & 40.3  & 33.3  & 54.1 \\
 \textbf{DEM} \cite{dem}  & 61.9  & 51.7  &  67.1 \\

\textbf{GFZSL}\cite{verma2017simple}  & 63.1  & 49.2   & 67.0 \\
 \textbf{CVAE-ZSL}\cite{mishra2017generative}  & 61.7  & 52.1  &  65.8 \\

 \hline
 \textbf{W/O ADA (Ours)} & \textbf{63.3} & \textbf{70.9}  &  \textbf{70.4} \\

 \hline
 \end{tabular}
 \caption{Zero Shot Learning Accuracy on the SUN, CUB, and AWA2 dataset. Here PS is the proposed split recently adopted in the ZSL community after \cite{xian2018zero}. The results reported in the table for the other approaches were taken from the Table 3 of \cite{xian2018zero} }
 \label{tab:zsl}
 \vspace{-1em}
\end{table}
\begin{table}[h!]
\small
 \centering
 \addtolength{\tabcolsep}{7pt}
 \begin{tabular}{|l| c |c | c | c |}
 \hline
 \textbf{Method}  & {\textbf{SUN}} & {\textbf{CUB}} &  {\textbf{AWA2}} \\
  \hline
 \textbf{DSRL}\cite{ye2017zero_DOMAIN} & 56.8 & 48.7 & 72.8 \\
 \textbf{ALE} \cite{akata2013label}  & 55.7 & 54.5 & 70.7\\
 \textbf{GFZSL} \cite{verma2017simple} &  64.2 &  50.5   & \textbf{78.6} \\
  \hline
 \textbf{With ADA (Ours)} & \textbf{65.5}  & \textbf{74.2} &  \textbf{78.6} \\
 \hline
 \end{tabular}
 \caption{Transductive Zero-Shot Learning results on the SUN, CUB, and AWA2 dataset. Transductive setting for our model corresponds to ADA. We note that the compared results are reported using the same ResNet101 feature and same train-test split. The results are taken from~\cite{xian2018zero} paper which has evaluated the models with ResNet101 features. }
 \label{tab:adversarialzsl}
 \vspace{-1.5em}
\end{table}

\subsection{Domain Adaption}
In ZSL, since $S\cap U=\phi$, there is a high probability that the seen and
unseen data do not come from the same underlying domain. This implies that the estimated parameters for the unseen classes, based on the training data of the seen class are likely to deviate from their optimal values. To this end, we propose an Adversarial Domain Adaptation (ADA) method (refer section \ref{sec:ADA}) to explicitly handle the domain shift problem.

In Table-\ref{tab:adversarialzsl}  we show the results of the proposed ADA method and compare against the previous transductive setting approaches. The result of ALE \cite{akata2013label} and GFZSL \cite{verma2017simple} are taken from the Figure 8 of \cite{xian2018zero}. Here we observe that using the domain adaption method boosts the generative model's performance. In the case of the AWA2 dataset without domain adaption, the top-1 accuracy was 70.4\% while with the domain adaption it rises to 78.6\%. A similar pattern is observed for the CUB (3.3\% improvement) and SUN dataset also. The domain shift in SUN dataset is ameliorated by the presence of a large number of training and testing classes and hence we see a smaller increment after ADA.

Moreover, our ADA method can minimize the domain shift (apparent in the Figure \ref{fig:Clustering} (a),(b)) in accordance with the clusters allotted by the base ZSL model. We can see that the model associates wrong clusters for only two classes owing to the low prediction accuracy of the base ZSL model for these classes which, itself is due to a strong overlap in test clusters of these classes. Thus, a reduction in label corruption will further improve the domain matching.

\subsection{Ablation Study}

In this section, we compare variants of our proposed approach through an ablation study to empirically analyze the benefits of each component. In particular, we check whether enforcing cyclic consistency leads to better performance than the vanilla adversarial loss, whether incorporating deep classifiers in architecture leads to improved performance, and whether adversarial domain adaptation is required for domain shift minimization for training the deep classifiers.
\subsubsection{Experimental Setup}
We have kept the pre-trained base ZSL model the same for consistent ablation results. The different variants of our model for the ablation study are described below:

\begin{itemize}
    \item \textbf{Std DA}: To test the relative importance of adversarial domain adaptation and hence domain shift minimization, we trained a deep classifier (with the same architecture as other variants) on the labeled samples synthesized from our generative model (base ZSL model) and the unlabelled test data with its pseudo labels.
    \item \textbf{Vanilla ADA}: This domain adaptation model comprises of a single generator and discriminator augmented with a classifier, where the generator maps the source domain to the target domain. For effective comparison, we used the same architecture of generators, discriminators and classifiers for ablation and experimental evaluations.
    \item \textbf{CycleGAN w/o}: In this variant, we removed the classifiers $C^T$ and $C^S$ associated with our proposed ADA model. Hence, the adversarial architecture is similar to CycleGAN which comprises of two generators and associated two discriminators.
    \item \textbf{Ours}: This is our proposed ADA model, defined in section \ref{sec:ADA}

\end{itemize}

We employ two different techniques for predicting the class labels. For the above-defined variants which have a trainable classifier in them like \textit{vanilla ADA} and \textit{Ours}, we report the class averaged top-1 accuracy of the predictions from the classifier attached to the discriminator (referred as \textbf{M1} in \textit{table \ref{tab:ablationzsl}}). For the approach \textit{Ours}, classifier $C^T$ ( mapping the source domain to the target domain ) is used for class label predictions.

We also report the 1 nearest neighbor classification accuracies using the Gaussian distance between the class conditionals mapped to the target domain by the generator and the test data feature (referred to as \textbf{M2} in Table \ref{tab:ablationzsl}). This method predicts the most probable class via the mixture of class conditionals, in a similar way like the base ZSL model in the inductive setting. To generate the mapped cluster prototypes for each class conditional, we sample the data points from its class conditional distribution, transform them into another target domain (test domain) via the generator of ADA and then extract the required statistics (mean of the cluster for our case) from the new distribution. Like ADA experimental setup, the learned covariance matrix is not changed after domain adaptation.

For the variant without a trainable classifier, \textit{CycleGAN w/o}, we only use the later method (method M2)  for evaluating the accuracies. Also, note that due to the absence of any adversarial generator in the variant \textit{Std DA}, the \textbf{M2} accuracy is computed in the exact same way as in our base inductive ZSL model.
\begin{table}[h!]
\scriptsize
 \centering
 \addtolength{\tabcolsep}{0.2pt}
 \begin{tabular}{|l| c | c | c | c | c | c |}
 \hline
   & \multicolumn{2}{c|}{\textbf{SUN}} & \multicolumn{2}{c|}{\textbf{CUB}} &  \multicolumn{2}{c|}{\textbf{AWA2}} \\
  \hline
  \textbf{Variant} & \textbf{M1}& \textbf{M2}& \textbf{M1} & \textbf{M2} & \textbf{M1} & \textbf{M2} \\
  \hline

 \textbf{Std DA}  & 64.8  & NA & 72.2 & NA & 71.3 & NA\\
 \hline
 \textbf{Vanilla ADA}  & 64.9  & 47.1 & 71.5 & 57.8 & 77.3 & 56.1\\
 \hline
 \textbf{CycleGAN w/o}  & NA  &  \textbf{57.2}  & NA & \textbf{68.4} & NA & \textbf{75.8}\\
  \hline
 \textbf{Ours} & \textbf{65.5}  & {55.8} & \textbf{74.2}  & 67.5 & \textbf{78.6} & 74.9\\
 \hline
 \end{tabular}
 \caption{Ablation study on ZSL with splits proposed in \cite{xian2018zero}}
 \label{tab:ablationzsl}
 \vspace{-3.5em}
\end{table}

\subsubsection{Analysis}
When we compare the performance of \textit{Std DA} with the base inductive ZSL model (results in Table \ref{tab:zsl}), we only see a marginal performance increase. During the experiments, the classifiers initially came close to our benchmark model (\textit{Ours}), but soon converged to a sub-optimum where they mimicked the accuracy of pseudo-labels provided by the base ZSL model. Owing to the domain shift, the classifier was not able to transfer the supervision from generated samples to the test data. This supports the claim that adversarial networks reduce domain mismatch, precluding the classifiers from converging at pseudo-labels.

The addition of trainable classifiers with ADA gave a heavy accuracy boost. This is mostly because of the higher expressivity and generalizability of such neural net classifiers as compared to nearest neighbor based classifiers. This is empirically suggested by diminished performance of about 3-10\% on various datasets in \textit{CycleGan w/o} wrt \textbf{M1} accuracy of \textit{Ours}. The addition of classification loss term does reduce the linear separability (reduction in \textbf{M2} score of \textit{Ours} vs \textit{CycleGAN w/o}) but the performance gain from classifiers overshadows this degradation.

Cyclic consistency further restrains the output space of the generator which drastically improves the linear separability of the generated data points (\textbf{M2} score of \textit{Ours}). This causes the proposed model to perform better than standard adversarial architecture using a similar classifier. This is apparent when we compare \textbf{M2} accuracies of \textit{vanilla ADA} with \textit{ours}. Even though the \textbf{M1} accuracies of these two models differ by about 1-2\%, the drop in \textbf{M2} accuracies are severe. Since nearest neighbor models rely on linear separability they suffer with as large as 10-30\% drop.This is also apparent in \textit{figure \ref{fig:Clustering}} t-SNE plots.

We can safely conclude, adding adversarial domain adaptation to the generative ZSL framework allows us to leverage the expressivity of neural net classifiers to classify novel classes while being trained only using the labels from seen classes. The adversarial adaptation minimizes the domain shift which is a crucial requirement for classifiers to transfer knowledge from the synthesized data and hence helps to train incisive classifiers that do not face the hubness issue, unlike distance-based nearest neighbor classifiers.

\section{Conclusion}
In this paper, we address the issue of domain shift between the distributions of the seen and unseen classes in zero-shot learning. We adopt an end-to-end approach for generative modeling that captures non-linear dynamics better as compared to previous state-of-the-art approaches. The proposed approach first learns the class conditional distributions for both the seen and unseen classes by leveraging the data from only the seen classes. Following this, we explicitly minimize the domain shift between the estimated unseen class distributions and the true unseen class distributions by using a cyclic consistency based adversarial scheme. We show through detailed experimentation, that our proposed generative model, although much simpler than GAN/VAE based frameworks, outperform existing models in the ZSL setting. Also, we show that our scheme of minimizing domain shift significantly improves performance, as compared to the transductive setting methods adopted by previous approaches. The generative framework can in principle assume any form, some of the popular ones being GAN and VAE based models. However, they lack explainability and they require further sampling to extract statistics like class variance. A larger intra-class variance would be an outcome of larger variations in the visual appearances of class attributes and hence the samples would be harder to classify together. An interesting future direction can be to use these statistics to model selective attention mechanisms or training with hard negative mining.

\textbf{Acknowledgements:} VKV acknowledges support from Visvesvaraya PhD Fellowship and PR acknowledges support from Visvesvaraya Young Faculty Fellowship.

{\small
\bibliographystyle{ieee}