\section{Introduction}

The last two decades have seen the development of several complexity dichotomies for pattern counting problems in graphs, including full classifications for counting \emph{subgraphs}, \emph{induced subgraphs}, and \emph{homomorphisms} from fixed computable pattern classes $\mathcal H$.
The input to such problems is a \emph{pattern} graph $H \in \mathcal H$ and an unrestricted \emph{host} graph $G$;
the task is to count the relevant occurrences of $H$ in $G$.
Depending on $\mathcal H$, these problems are known to be either polynomial-time solvable or $\cc{\#W[1]}$-hard when parameterized by $|V(H)|$.
The latter rules out polynomial-time algorithms under the complexity assumption $\cc{FPT} \neq \cc{\#W[1]}$.

In this paper, we focus on counting \emph{subgraphs} from any fixed graph class $\mathcal H$.
On the positive side, given a pattern graph~$H\in\mathcal H$ whose smallest vertex-cover has size $\operatorname{vc}(H)$ and an $n$-vertex host graph~$G$, there are known $O(n^{\operatorname{vc}(H)+1})$ time  algorithms~\cite{DBLP:journals/siamcomp/WilliamsW13,DBLP:journals/siamdm/KowalukLL13,DBLP:conf/focs/CurticapeanM14} to count subgraphs of $G$ that are isomorphic to $H$:
First, find a minimum vertex-cover~$C$ of~$H$ using exhaustive search. Then, iterate over all possible embeddings~$f$ of~$H[C]$ into~$G$ and count the possible extensions of $G[f(C)]$ to a full copy of~$H$.
Complementing this algorithm, an almost matching running time lower bound of $n^{\Omega(\operatorname{vc}(H) / \log \operatorname{vc}(H))}$ under the exponential-time hypothesis~($\ensuremath{{\mathrm{ETH}}}$) is also known~\cite{DBLP:conf/stoc/CurticapeanDM17}.
Thus, assuming $\ensuremath{{\mathrm{ETH}}}$ or $\cc{FPT} \neq \cc{\#W[1]}$, the problem $\#\ensuremath{{\mathrm{Sub}}}(\mathcal H)$ of counting subgraphs from a fixed class~$\mathcal H$ is polynomial-time solvable if and only if the vertex-cover numbers (or equivalently, the maximum matching sizes) of the graphs in~$\mathcal H$ are bounded by a constant.
The rightmost column of \cref{fig: overview} visualizes this situation.

Turning from counting to the problem $\ensuremath{{\mathrm{Sub}}}(\mathcal H)$ of \emph{detecting} subgraphs from fixed classes~$\mathcal H$, the picture is less clear.
Evidence points at three strata of complexity:
Define the \emph{matching-split number} of $H$ to be the minimum number of vertices whose deletion turns $H$ into a matching, that is, a graph of maximum degree $1$.
Jansen and Marx~\cite{JansenMarx} show that, if this number is bounded in a graph class~$\mathcal H$, then $\ensuremath{{\mathrm{Sub}}}(\mathcal H)$ is polynomial-time solvable.
For classes~$\mathcal H$ of bounded tree-width, it is known~\cite{DBLP:conf/wg/PlehnV90,DBLP:journals/algorithmica/AlonYZ97,DBLP:journals/jcss/FominLRSR12} that the problem $\ensuremath{{\mathrm{Sub}}}(\mathcal H)$ is fixed-parameter tractable when parameterized by~$|V(H)|$.
For pattern classes $\mathcal H$ of unbounded tree-width, it is conjectured that $\ensuremath{{\mathrm{Sub}}}(\mathcal H)$ is W[1]-hard---so far, this hardness has only been established for cliques, bicliques~\cite{DBLP:journals/jacm/Lin18}, grids~\cite{DBLP:conf/wg/ChenGL17}, and less natural graph classes.
The leftmost column of \cref{fig: overview} visualizes the situation.

We propose to study an intermediate setting between decision and counting, namely, counting subgraph patterns \emph{modulo fixed integers} $q \in \mathbf{N}$.
Modular counting has a tradition in classical complexity theory, where the complexity classes $\ensuremath{{\mathrm{Mod}}}_q \cc{P}$ for $q \in \mathbf{N}$ capture problems that ask to count accepting paths of polynomially time-bounded non-deterministic Turing machines modulo~$q$.
In particular, the class $\ensuremath{{\mathrm{Mod}}}_2 \cc{P}$ (better known as~$\cc{\oplus P}$) plays a central role in the proof of Toda's theorem~\cite{Toda}.
Several (partial) classification results for frameworks of modular counting problems are known; this includes homomorphisms to fixed graphs~\cite{DBLP:journals/toc/FabenJ15,DBLP:journals/toct/0001GR14,DBLP:journals/toct/0001GR16,DBLP:conf/mfcs/KazeminiaB19,doi:10.1137/1.9781611976465.137}, constraint satisfaction problems~\cite{DBLP:journals/corr/abs-0809-1836,DBLP:conf/stacs/GuoHLX11}, and Holant problems~\cite{DBLP:conf/focs/CurticapeanX15}.

\cref{fig: overview} summarizes our understanding.
If the vertex-cover number is bounded, the polynomial-time algorithms (regions 7 and 8) follow from the algorithm for $\#\ensuremath{{\mathrm{Sub}}}(\mathcal H)$ described above and require no further attention.
Our paper is concerned with the remaining regions~\mbox{1--6}.

As argued above, matchings play a central role in decision and counting, so it is natural that they reprise their role in modular subgraph counting:
On the positive side, there are known polynomial-time algorithms for counting matchings of a given size modulo fixed powers of two. (For bipartite graphs and counting modulo $2$, this essentially follows from the fact that determinant and permanent coincide modulo~$2$.)
On the negative side, if $q$ is not a power of two, counting matchings modulo $q$ is known to be $\ensuremath{{\mathrm{Mod}}}_p \cc{P}$-complete for any odd prime $p$ dividing $q$.
We establish a parameterized analogue of this fact:
Let $\ensuremath{{\mathrm{Mod}}}_q\cc{W[1]}$ be the class of parameterized problems that are fpt-reducible to counting $k$-cliques modulo~$q$.
We show that counting $k$-matchings (that is, sets of $k$ pairwise disjoint edges) in graphs modulo fixed odd primes $q \in \mathbf{N}$ is $\ensuremath{{\mathrm{Mod}}}_q\cc{W[1]}$-hard.
In our proof, modular counting allows us to sidestep the algebraic machinery from previous works~\cite{DBLP:conf/focs/CurticapeanM14,DBLP:phd/dnb/Curticapean15,DBLP:conf/stoc/CurticapeanDM17}, resulting in a surprisingly simple and self-contained argument.
\begin{theorem}\label{thm: hardness-matchings}
  For any integer $q \in \mathbf{N}$ containing an odd prime factor $p$,
  counting $k$-matchings modulo $q$ is $\ensuremath{{\mathrm{Mod}}}_p\cc{W[1]}$-hard under Turing fpt-reductions and admits no $n^{o(k/\log k)}$ time algorithm under $\ensuremath{{\mathrm{ETH}}}$.
\end{theorem}

\begin{figure}[t]
  \begin{tikzpicture}[xscale=2.4]
    \draw [decorate,decoration={brace,amplitude=5pt}] (1,4.2) --  node [midway, yshift=1.2em] {This paper} (3,4.2);

    \fill [green!20] (0,0) rectangle (4,1);
    \fill [green!20] (0,1) rectangle (2,2);
    \fill [red!20] (1,3) rectangle (4,4);
    \fill [red!20] (3,1) rectangle (4,4);
    \fill [yellow!20] (0,2) rectangle (1,3);

    \draw (0,0) rectangle (1,4);
    \draw (1,0) rectangle (2,4);
    \draw (2,0) rectangle (3,4);
    \draw (3,0) rectangle (4,4);
    \draw [dashed] (0, 1) --(4,1);
    \draw [dashed] (0, 2) --(4,2);
    \draw [dashed] (0, 3) --(4,3);

    \node at (-.1,.5) [anchor = east, align=center] {vertex-cover\\number};
    \node at (-.1,1.5) [anchor = east, align=center] {matching-split\\number};
    \node at (-.1,2.5) [anchor = east] {tree-width};
    \node at (.5, -.1) [anchor = north] {decision};
    \node at (1.5, -.1) [anchor = north,align=center, text width = 2cm] {counting mod $2^t$};
    \node at (2.5, -.1) [anchor = north, align=center,text width = 2cm] {counting mod odd $q$};
    \node at (3.5, -.1) [anchor = north] {counting};

    \node (pt) at (1.5, 2.4) [inner sep=1pt, black!50!red] {\small\emph{paths}};
    \node (p) at (1.2, 2.85) [inner sep=0, red!20] {$\bullet$};
    \draw [black!50!red] (p) -- (pt);
    \node (mt) at (2.5, 1.4) [inner sep=1pt, black!50!red] {\small\emph{matchings}};
    \node (m) at (2.95, 1.35) [inner sep=0, red!20] {$\bullet$};
    \draw [black!50!red] (m) -- (mt);
    \node (ct) at (0.5, 3.7) [inner sep=1pt, black!50!red] {\small\emph{grids}};
    \node (c) at (0.95, 3.8) [inner sep=0, red!20] {$\bullet$};
    \draw [black!50!red] (c) -- (ct);
    \node (gt) at (0.5, 3.3) [inner sep=1pt, black!50!red] {\small\emph{bicliques}};
    \node (g) at (0.95, 3.2) [inner sep=0, red!20] {$\bullet$};
    \draw [black!50!red] (g) -- (gt);

    \node at (3.9, 3) [black!50!red, fill=red!20, anchor=east] {$W[1]$-hard};
    \node at (.1, 1) [black!50!green, fill=green!20, anchor=west, inner sep=0pt] {Polynomial time};
    \node at (.1, 2.5) [black!50!yellow, anchor=west] {FPT};

    \foreach \position/\text in { (1,0)/7, (2,0)/8, (1,1)/5, (2,1)/6, (1,2)/3, (2,2)/4, (1,3)/1, (2,3)/2 }
      \node at \position [xshift = 2.4cm, yshift = .8cm, anchor=east] {\small\text};
  \end{tikzpicture}

  \caption{\label{fig: overview}%
  An overview over known results and our new results.
  The columns correspond, from left to right, to the problem types $\ensuremath{{\mathrm{Sub}}}(\mathcal H)$, $\ModSub{\mathcal H}{2^t}$, $\ModSub{\mathcal H}{q}$ for $q\ne 2^t$, or $\#\ensuremath{{\mathrm{Sub}}}(\mathcal H)$; our results are depicted in the two middle columns.
  The rows correspond, from bottom to top, to requiring $\mathcal H$ to have bounded vertex-cover number, matching-split number, tree-width, or no requirement at all.
  The complexity along each row is monotone: By \cref{lem:decision-to-parity}, decision is no harder than modular counting, and modular counting trivially is no harder than counting.
  %
  Our results are depicted in the middle two columns:
  Regions~1 and~2 are \cref{cor: hardness-tw}.
  Region~5 is \cref{thm: algorithm-split}.
  Regions~7 and~8 already follow from~\cite{DBLP:journals/siamcomp/WilliamsW13}.
  The point in region~6 is \cref{thm: hardness-matchings}, and the point in region~3 is \cref{thm: hardness-paths}.
  We view the hardness of these points as evidence to conjecture their enclosing regions to be hard, see \cref{conjecture}.
  }
\end{figure}

Known arguments from Ramsey theory (see~\cite[Section~5]{DBLP:journals/corr/CurticapeanM14}) extend \cref{thm: hardness-matchings} from matchings to $\ModSub{\mathcal H}{q}$ for any hereditary class~$\mathcal H$ of unbounded vertex-cover number.
This suggests that modular subgraph counting may only become tractable when the modulus is a power of two.
Indeed, we show that patterns of matching-split number $s$ can be counted modulo $q = 2^t$ in time $n^{O(t4^s)}$.
To prove this, we follow the general idea of the bounded vertex-cover number algorithm for $\ensuremath{{\mathrm{Sub}}}(\mathcal H)$ outlined before, and we reduce to counting matchings modulo powers of two.
This however requires us to overcome technical complications to avoid unwanted cancellations.
Overall, we obtain:
\begin{theorem}\label{thm: algorithm-split}
  There is an algorithm that, given a graph $H$ of matching-split number $s \in \mathbf{N}$ and an $n$-vertex graph $G$,
  computes the number of $H$-isomorphic subgraphs of~$G$ modulo $2^t$ in time $n^{O(t4^s)}$.
\end{theorem}
We complement this result in two ways:
First, we observe that $\oplus\ensuremath{{\mathrm{Sub}}}(\mathcal H)$ is $\oplus\W$-complete for pattern classes $\mathcal H$ of unbounded tree-width;
this follows directly from previous hardness proofs for $\#\ensuremath{{\mathrm{Sub}}}(\mathcal H)$.
More interestingly, we establish the $\oplus\W$-completeness of counting $k$-paths modulo~$2$ in undirected graphs, thus solving an open problem from~\cite{DBLP:conf/icalp/BjorklundDH15}, where this problem was considered in the context of Hamiltonian cycle detection, following~\cite{BjorklundHusfeldt}.
\begin{theorem}\label{thm: hardness-paths}
  Counting $k$-paths modulo $2$ is $\oplus\W$-complete.
\end{theorem}
This result adds to a rich range of previous work on the $k$-path problem, and is of interest outside our framework.
Bodlaender~\cite{BODLAENDER19931} and Monien~\cite{MONIEN1985239} showed that \emph{finding} a $k$-path is fixed-parameter tractable.
In contrast, Flum and Grohe~\cite{DBLP:journals/siamcomp/FlumG04} showed that \emph{exactly counting} $k$-paths is $\cc{\#W[1]}$-hard.
Nevertheless, Arvind and Raman~\cite{ArvindRaman} showed that \emph{approximately counting} $k$-paths, which corresponds to computing the most significant bit(s) of the number of~$k$-paths, is fixed-parameter tractable.
Our \cref{thm: hardness-paths} suggests that the \emph{least} significant bit of the number of $k$-paths is hard to compute.
This is surprising, because some of the most influential fpt-algorithms for finding a $k$-path work over characteristic~$2$,
based on the group algebra framework introduced by Koutis~\cite{Koutis08}.

\medskip

Let us conclude with a general remark on the techniques used in this paper:
Recent works successfully exploited a connection between subgraph counts and (linear combinations of) homomorphism counts to obtain algorithms and hardness results~\cite{DBLP:conf/stoc/CurticapeanDM17,DBLP:conf/esa/Roth17,DBLP:conf/mfcs/DorflerRSW19,DBLP:journals/algorithmica/RothS20,DBLP:conf/focs/Roth0W20}.
For example, the number of $k$-matchings in a graph $G$ is a linear combination of homomorphism counts from $f(k)$ fixed graphs.
Insights on the complexity of counting the homomorphisms occurring in this linear combination then lead to complexity results for counting $k$-matchings.
This connection however does not readily transfer to modular counting, as the relevant linear combinations (which involve rational coefficients) may be \emph{undefined} modulo~$p$. We therefore prove \cref{thm: hardness-matchings,thm: algorithm-split,thm: hardness-paths} using more combinatorial approaches.

\section{Preliminaries}\label{sec: preliminaries}
Unless otherwise stated, we consider finite, undirected, simple graphs without self-loops.

\paragraph*{Subgraph problems}

A \emph{homomorphism} from graph $H$ to graph $G$ is a mapping $\varphi\colon V(H)\to V(G)$ such that $\{\varphi(u),\varphi(v)\} \in E(G)$ for each $\{u,v\}\in E(H)$.
An \emph{embedding} is an injective homomorphism, and we let $\ensuremath{{\mathrm{Emb}}}(H,G)$ denote the set of embeddings from~$H$ to~$G$.
An \emph{isomorphism} is a bijective homomorphism,
and an \emph{automorphism} is an isomorphism from $H$ to itself.
The set of all automorphisms of $H$ is called $\ensuremath{{\mathrm{Aut}}}(H)$, and forms a group when endowed with function composition~$\circ$.

We let $\ensuremath{{\mathrm{Sub}}}(H,G)$ be the set of all $H$-subgraphs of~$G$, that is, the set of all $H'$ with $V(H')\subseteq V(G)$ and $E(H')\subseteq E(G)$ such that $H'$ is isomorphic to~$H$.
This terminology fixes the possible confusion about isomorphic copies of subgraphs:
For example, there is exactly one $K_k$-subgraph in $K_k$,
but there are $k!$ embeddings.
The \emph{subgraph problem $\ensuremath{{\mathrm{Sub}}}$} is given a pair $(H,G)$ to decide whether $G$ has at least one $H$-subgraph.
The \emph{subgraph counting problem $\#\ensuremath{{\mathrm{Sub}}}$} is given a pair $(H,G)$ to determine the number of $H$-subgraphs in $G$.

For a graph class~$\ensuremath{{\mathcal{H}}}$, we write $\#\ensuremath{{\mathrm{Sub}}}(\ensuremath{{\mathcal{H}}})$ for the restricted
problem where the input $(H,G)$ is promised to satisfy $H\in \ensuremath{{\mathcal{H}}}$.
For $q\in\mathbf{Z}_{\ge 2}$, the \emph{modular subgraph counting problem $\ModSub{\ensuremath{{\mathcal{H}}}}{q}$} is the problem to compute the number of $H$-subgraphs modulo~$q$.
In the special case with~$q=2$, we write $\oplus\ensuremath{{\mathrm{Sub}}}$.

It will be useful to consider \emph{colorful} subgraph problems, where~$G$ is $H$-colored, that
is, there is a given homomorphism~$c:V(G)\to V(H)$.
Due to the homomorphism property, we allow edges $\{u,v\}\in E(G)$ only if the
corresponding colors satisfy $\{c(u), c(v)\} \in E(H)$.
A subgraph~$H'$ of an $H$-colored graph~$G$ is \emph{vertex-colorful} if $c$ is
bijective on $V(H')$.
Let $\ensuremath{{\mathrm{VertexColorfulSub}}}(H,G)$ be the set of vertex-colorful subgraphs~$H'$ for which $c$ is
an isomorphism from $H'$ to $H$.
The corresponding computational problems are defined analogously to the
uncolored case; the input consists of a graph~$G$ together with an
$H$-coloring~$c$.

\paragraph*{Background from complexity theory}
A \emph{parameterized counting problem} is a pair $(f,\kappa)$ of functions $f,\kappa:\set{0,1}^\ast\to\mathbf{N}$ where $\kappa$ is computable.
A \emph{parsimonious fpt-reduction} from a parameterized counting problem~$(f,\kappa)$ to a parameterized counting problem~$(g,\iota)$ is a function~$R$ with the following properties:
(i) $f(x)=g(R(x))$ for all $x\in\set{0,1}^\ast$,
(ii) $\iota(R(x))$ is bounded by a computable function in~$\kappa(x)$, and
(iii) the reduction is computable in time $h(\kappa(x))\op{poly}(\abs{x})$ for some computable function~$h$.
A~\emph{Turing fpt-reduction} may query the oracle multiple times for instances whose parameter is bounded by a function of the input parameter, and combine the query answers in fpt-time to produce the correct output.
Moreover, reductions can also be \emph{randomized}, in which case we require that their error probability is bounded by a small constant.

The \emph{exponential-time hypothesis} (ETH) postulates the existence of some~$\varepsilon>0$ such that no algorithm solves $n$-variable~$3$-CNF formulas in time~$O(2^{\varepsilon n})$. We write for short that $3$-CNF-SAT does not have $2^{o(n)}$-time algorithms, and we also disallow bounded-error randomized algorithms.

\paragraph*{Modular counting}
For our purposes, we define the class $\ensuremath{{\mathrm{Mod}}}_q\cc{W[1]}$ as the class of all parameterized problems $(f,\kappa)$ with $f\colon \Sigma^\ast\to\{0,\dots,q-1\}$ such that $(f,\kappa)$ has a parsimonious fpt-reduction to the problem of counting $k$-cliques modulo~$q$.
For $q=2$, it was shown in~\cite{DBLP:conf/icalp/BjorklundDH15} that
all problems in $\cc{W[1]}$ admit randomized fpt-reductions to problems in $\oplus\W$.
Another result~\cite[Lemma~2.1]{DBLP:conf/soda/WilliamsWWY15} yields the corresponding generalization for all $q>2$.
We use the following analogous proposition for the vertex-colorful subgraph problem, proven in the appendix.

\begin{lemma}\label{lem:decision-to-parity}%
  For any integer $q \geq 2$,
  there is a randomized Turing fpt-reduction from the problem~$\ensuremath{{\mathrm{VertexColorfulSub}}}$ to the problem~$\#\ensuremath{{\mathrm{VertexColorfulSub}}}\bmod q$.
  On input $(H,G)$, the reduction only queries instances with the same pattern $H$.
\end{lemma}

Our work relies on the following hardness result for parameterized modular subgraph counting, which follows easily from known results on the colorful subgraph decision problem~\cite{DBLP:journals/tcs/DalmauJ04,DBLP:journals/toc/Marx10}. See the appendix for a proof.

\begin{lemma}\label{cor: hardness-tw}
  Let $\ensuremath{{\mathcal{H}}}$ be a graph family of unbounded tree-width and let~$q$ be an integer with $q\ge 2$.
  Then $\#\ensuremath{{\mathrm{VertexColorfulSub}}}(\ensuremath{{\mathcal{H}}})\bmod q$ parameterized by $k=\abs{E(H)}$ is $\ensuremath{{\mathrm{Mod}}}_q\cc{W[1]}$-hard under parsimonious fpt-reductions.
  Moreover, if ETH is true, then the problem does not have an algorithm running
  in time $n^{o(k/\log k)}$, where~$n=\abs{V(G)}$.
\end{lemma}

\section{Hardness of counting $k$-matchings}\label{sec: hardness-matchings}

In this section, we prove \cref{thm: hardness-matchings}.
We first establish $\ensuremath{{\mathrm{Mod}}}_q\cc{W[1]}$-hardness of the problem $\#\pp{ColMatch}\bmod q$ for odd $q \geq 3$:
Given a graph $G$ with an edge-coloring $c: E(G) \to \mathcal C$ for some set of colors $\mathcal C$ with $\abs{\mathcal C}=k$,
this problem asks to count modulo $q$ the edge-colorful matchings in $G$.
These are the matchings that use each color in $\mathcal C$ exactly once.

\begin{lemma}\label{lem:reduction-to-colmatch}
  For any fixed integer $p$ with odd prime factor $q$,
  the problem ${\#\pp{ColMatch}\bmod p}$ is $\ensuremath{{\mathrm{Mod}}}_q \cc{W[1]}$-hard under parsimonious fpt-reductions and has no $n^{o(k/\log k)}$ time algorithm under $\ensuremath{{\mathrm{ETH}}}$.
\end{lemma}
\begin{proof}
  The class $\mathcal H$ of all $3$-regular graphs has unbounded tree-width, and hence by \cref{cor: hardness-tw}, the problem $\#\ensuremath{{\mathrm{VertexColorfulSub}}}(\ensuremath{{\mathcal{H}}})\bmod q$ is $\ensuremath{{\mathrm{Mod}}}_q \cc{W[1]}$-hard and hard under $\ensuremath{{\mathrm{ETH}}}$.
  We reduce it to $\#\pp{ColMatch}\bmod q$, implying the hardness of $\#\pp{ColMatch}\bmod p$.
  Let $H\in\ensuremath{{\mathcal{H}}}$ and~$G$ be the input for the reduction with $k=\abs{V(H)}$ and $H$-colored $G$,
  and let $\setc{V_a}{a \in V(H)}$ be the color classes of $G$,
  with edge-sets $E_{a,b}(G)$ for $ab \in E(H)$.
  Using the gadgets from \cref{fig:matchings mod p}, we construct a graph~$G'$:

  \begin{enumerate}
    \item Each vertex $u \in V(G)$ is replaced by three vertices $u_1$, $u_2$, and $u_3$.
    We insert a \emph{consistency gadget} $Q_u$ at these vertices by adding $q$ gadget vertices, connecting $u_2$ and $u_3$ to all gadget vertices, and $u_1$ to the first $(q+1)/2$ gadget vertices.
    For $S \subseteq \set{u_1,u_2,u_3}$, let $m_S$ count the matchings in $Q_u$ that match precisely $S$;
    it can be checked that
    \begin{equation}
    \label{eq: consistency-gadget-extensions}
      m_S
      \equiv_q
      \begin{cases}
        1 & \text{if $S$ is $\emptyset$ or $\set{u_1,u_2,u_3}$,}
        \\
        0 & \text{if $S$ is $\set{u_2}$, $\set{u_3}$, or $\set{u_2,u_3}$.}
      \end{cases}
    \end{equation}
    We explicitly ignore the other three cases for $S$, as they will not be relevant.

  \item For $\{a,b\} \in E(H)$, suppose that $a$ is the $j$th vertex incident to $b$, and that $b$ is the $i$th vertex incident to~$a$.
    For each edge $\{u,v\} \in E_{a,b}(G)$ with $u\in V_a$ and $v\in V_b$, we insert an AND-gadget $A_{uv}$ at $\set{u_i,v_j}$. Then for any set $S\subsetneq \set{u_i,v_j}$, the number of edges in $A_{uv} - S$ is divisible by $q$, whereas $A_{uv} - \set{u_i,v_j}$ has exactly one edge.
    \item The edge-colors of $G'$ are defined as follows:
    For~$u\in V_a$ with $a \in V(H)$, we assign color $(\mathtt{CONS},a,i)$ to all edges of $Q_u$ incident to vertex~$u_i$, for $i\in\set{1,2,3}$.
    For each $ab \in E(H)$, we assign color $(\mathtt{AND},ab)$ to all edges in AND-gadgets between edges in $E_{a,b}(G)$.
    Overall, we have $k'=3k+3k/2$ colors.
  \end{enumerate}

  \begin{figure}[tp]

  	    \centering
  	\begin{subfigure}[t]{0.4\textwidth}
  		\centering
  		\input{gadgets_and.pdf_tex}
  		\subcaption{The AND-gadget, shown here for $q=5$.
  			In general, the upper $u_i,v_j$-path always has $3$ edges; both \emph{external vertices} $u_i$ and~$v_j$ have $q-2$ neighboring leaves.
  			If exactly~$0$ or~$1$ external vertices are
  			removed, this graph has $0$ edges modulo~$q$; if both vertices are removed,
  			the graph has $1$ edge.}
  	\end{subfigure}%
  	\hspace{1cm}
  	\begin{subfigure}[t]{0.4\textwidth}
  		\centering
  		\input{gadgets_cons.pdf_tex}
  		\subcaption{The consistency gadget contains $q$ gadget vertices, shown here for $q=5$. The number of matchings of size~$0$ and $3$ equals~$1$ modulo~$q$, and if $u_1$ is deleted, the number of non-empty edge-colorful matchings equals~$0$ modulo $q$.}
  	\end{subfigure}
  	\caption{\label{fig:matchings mod p}%
    The two gadgets used in the proof of \cref{lem:reduction-to-colmatch}.}
  \end{figure}

  Every $H$-copy $F$ in $G$ induces a set $\mathcal{M}_F$ of colorful matchings in $G'$. We describe this set in the following, show that $|\mathcal{M}_F| \equiv_q 1$, and that $\mathcal{M}_F$ and  $\mathcal{M}_{F'}$ are disjoint for $F \neq F'$.
  \begin{itemize}
    \item For each $v \in V(F)$, match all of $\set{v_1,v_2,v_3}$ within $Q_v$.
    For fixed $v$, the number of possible matchings in $Q_v$ is $m_{\set{v_1,v_2,v_3}} \equiv_q 1$ by~\eqref{eq: consistency-gadget-extensions}.
    Let $\mathcal{Q}_F$ denote the set of all matchings that can be obtained by the previous step.
    Since matchings can be chosen independently for distinct $Q_v$, we obtain $|\mathcal{Q}_F| \equiv_q 1^{|V(F)|} \equiv_q 1$.
    %
  \item Any $M \in \mathcal{Q}_F$ can be extended to several colorful matchings by choosing one edge from each color $(\mathtt{AND},ab)$ for $\{a,b\} \in E(H)$.
    For each $\{u,v\} \in E(F)$, the AND-gadget $A_{uv}$ has exactly one such edge,
    while the other AND-gadgets of color $(\mathtt{AND},ab)$ have $0$ such edges modulo $q$.
    Hence, the number edges of color $(\mathtt{AND},ab)$ that can extend $M$ is $1$ modulo $q$.
    This implies that the overall number $r_M$ of matchings extending $M$ into a colorful matching is also $r_M \equiv_q 1^{|E(F)|} \equiv_q 1$.
  \end{itemize}
  Overall, every $H$-copy $F$ induces $\sum_{M \in \mathcal{Q}_F} r_M \,\equiv_q\, \sum_{M \in \mathcal{M}_F} 1 \,\equiv_q\, 1$ colorful matchings, so we indeed have $|\mathcal{M}_F| \equiv_q 1$. We also observe from the construction that $\mathcal{M}_F \cap \mathcal{M}_{F'} = \emptyset$ for distinct $H$-copies $F$ and $F'$. In the appendix, we use properties of the gadgets to prove that colorful matchings $M \notin \bigcup_{F} \mathcal{M}_F$ cancel modulo $q$.
\begin{claim}\label{claim: irrelevant-matchings}
  The number of colorful matchings $M$ that are not contained in $\mathcal{M}_F$ for any $H$-copy $F$ is divisible by $q$.
\end{claim}
Overall, we have shown that the number of $H$-copies in $G$ and the number of colorful matchings in $G'$ agree modulo $q$. As $G'$ can be computed in polynomial time and the parameter is increased only by a constant factor, the claimed hardness results follow.
\end{proof}

To prove \cref{thm: hardness-matchings}, it suffices to give an fpt-reduction from ${\#\pp{ColMatch}\bmod q}$ to counting $k$-matchings modulo $q$.
This is achieved by a standard inclusion-exclusion argument that can be found in the appendix.

\section{Counting matching-splittable subgraphs modulo $2^t$}\label{sec: algorithm}

In this section, we prove \cref{thm: algorithm-split} by describing an $n^{O(t4^s)}$-time algorithm for counting modulo $2^t$ the subgraphs of matching-split number $s$.
Our algorithm builds upon known algorithms for the decision and counting versions of subgraph problems;
we first review their underlying ideas and sketch our algorithm for \cref{thm: algorithm-split}.

\subsubsection*{Counting subgraphs of bounded vertex-cover number.}
The basic structure of our algorithm is similar to a known $O(n^{s+1})$ time algorithm~\cite{DBLP:journals/siamcomp/WilliamsW13,DBLP:journals/siamdm/KowalukLL13,DBLP:conf/focs/CurticapeanM14} for counting embeddings from~$H$ to~$G$ if $H$ has a vertex-cover $S\subseteq V(H)$ of size $s\in\mathbf{N}$. Counting embeddings is sufficient for counting subgraph copies, as we can first compute the number $\#\ensuremath{{\mathrm{Aut}}}(H)$ of automorphisms on~$H$ as $\#\ensuremath{{\mathrm{Emb}}}(H,H)$, and then use
\begin{equation}\label{eq: sub-emb-aut}%
\#\ensuremath{{\mathrm{Sub}}}(H,G) = \frac{\#\ensuremath{{\mathrm{Emb}}}(H,G)}{\#\ensuremath{{\mathrm{Aut}}}(H)}\,.
\end{equation}

Given $(H,G)$ with $h=\abs{V(H)}$ and $n=\abs{V(G)}$, the algorithm for computing $\#\ensuremath{{\mathrm{Emb}}}(H,G)$ first finds a minimum vertex-cover~$S$ of~$H$ in time~$h^{O(s)}$; then $I:=V(H)\setminus S$ is an independent set.
Then the algorithm enumerates all partial embeddings~$f$ from $H[S]$ to $G$, which takes time at most~$n^{O(s)}$.
Finally, for each~$f$, it remains to count all functions $g\colon I \to V(G)$ that extend~$f$ to a full embedding from $H$ to~$G$.
We observe that $g$ extends $f$ to a full embedding if and only if every vertex $u\in I$ maps via~$g$ to a vertex~$v=g(u)\in V(G)\setminus f(S)$ that satisfies the \emph{neighborhood constraint} $N_G(v)\cap f(S)\supseteq f(N_H(u))$.
Counting functions~$g$ with this property can be achieved (in a not completely obvious way) with dynamic programming;
we only need to know the number of vertices $v\in V(G)\setminus f(S)$ that have a specific neighborhood $N_G(v)\cap f(S)$, and for each~$f$, there are at most $2^s$ different possible such neighborhoods.
Overall, in~$n^{O(s)}$ time, we can compute the number $\#\ensuremath{{\mathrm{Emb}}}(H,G)$.

\subsubsection*{Detecting subgraphs of bounded matching-split number.}
Jansen and Marx~\cite{JansenMarx} extend the above approach and obtain an $n^{O(s)}$ time algorithm for the \emph{decision} problem $\ensuremath{{\mathrm{Sub}}}(H,G)$
when~$H$ has matching-split number~$s$.
In this case, we consider a \emph{splitting set} $S$ of size~$s$ instead of a vertex-cover, that is, the graph $M=H-S$ may have isolated edges besides isolated vertices.
Now the idea is to not only classify the vertices $v\in V(G)\setminus f(S)$ by their neighborhoods $N_v=N_G(v)\cap f(S)$, but to also classify the edges $\set{u,v}\in E(G-f(S))$ by their neighborhoods $\set{N_u,N_v}$.
It then remains to find a matching in~$G-f(S)$ that has as many isolated vertices and isolated edges as $H-S$, such that these vertices and edges satisfy the neighborhood constraints in~$f(S)$.
Jansen and Marx achieve this by reduction to a colored matching problem.

\subsubsection*{Our algorithm.}
In our algorithm for \cref{thm: algorithm-split}, we need to overcome two challenges:
\begin{enumerate}[(a)]
  \item\label{challenge-modulo} Since counting embeddings is algorithmically more straight-forward than counting subgraphs, we would like to count embeddings and divide by the number of automorphisms $\#\ensuremath{{\mathrm{Aut}}}(H)$ as in \cref{eq: sub-emb-aut}.
  However, since we are counting modulo~$2^t$, the number $\#\ensuremath{{\mathrm{Aut}}}(H)\bmod 2^t$ may be~$0$, and so the division in \cref{eq: sub-emb-aut} is impossible. (In fact, even even numbers $\#\ensuremath{{\mathrm{Aut}}}(H)$ have no inverse modulo $2^t$.)
  \item\label{challenge-count} When mimicking Jansen and Marx's detection algorithm, we cannot just \emph{count} the relevant matchings in $G-f(S)$, since counting perfect matchings is~$\cc{\#P}$-hard.
\end{enumerate}
Most of our effort focuses on overcoming~\eqref{challenge-modulo}:
In \cref{sub: rigidize-split}, we show that every graph~$H$ of matching-split number $s$ has a splitting set $R$ of size $O(s^2)$ that remains rigid under automorphisms, i.e., any automorphism $f$ of $H$ must
satisfy $f(R) = R$.
In \cref{sub: fixed split}, we show how to compute $\#\ensuremath{{\mathrm{Sub}}}(H,G)$ if
such a rigid splitting set~$R$ for~$H$ is given.
Rather than counting $H$-embeddings and attempting a division by $\#\ensuremath{{\mathrm{Aut}}}(H)$, we use the rigidity of $R$ to keep track of the automorphisms of~$H$ in a more explicit way.

To overcome~\eqref{challenge-count}, we use a determinant-based algorithm~\cite{HiraiNamba} to compute the Hafnian over a polynomial ring modulo~$2^t$ in \cref{sec:Hafnian}.
We then reduce our constrained matching counting problem to computing such Hafnians in \cref{sec:color-demands}.

\subsection{Rigidizing the splitting set}%
\label{sub: rigidize-split}

Let $H$ be a graph with a splitting set $S$ of size $s$, and let
$M=H-S$ be the remaining graph of maximum degree $1$; we speak of $M$ as a matching, even though it may contain isolated vertices.
An automorphism~$f$ of $H$ may map a vertex $v\in S$ in the splitting set to $f(v) \notin S$.
We show that if a splitting set of size $s$ exists then there is also a \emph{rigid} splitting set~$R$ of size~$O(s^2)$, i.e.,
such that every $f \in \ensuremath{{\mathrm{Aut}}}(H)$ satisfies $f(R) = R$.
In fact, the following algorithm can find such a set~$R$.

\begin{algorithm}{Rigidize}{$H$}{Given a graph $H$ of matching-split number $s$, this algorithm computes a rigid splitting set $R\subseteq V(H)$ of size $O(s^2)$.}
\item[R1] (Find small splitting set.)
  Using brute-force, compute a set~$S\subseteq V(H)$ of size~$s$ such that $H-S$ is a matching.
\item[R2] (Extend it to neighbors of low-degree vertices.)
  Let $D\subseteq V(H)$ be the set of all vertices whose degree in~$H$ is at most~$s+1$.
  Set $T:=S$.
  While there is an edge $\{u, v\}$ with $u\in T\cap D$ and $v\in\overline T$, add $v$ to~$T$.
\item[R3] (Refine it.)
  Set $R:=T$.
  For each component~$C$ of~$H[T\cap D]$ with at most two vertices,
  remove~$V(C)$ from~$R$.
\end{algorithm}
The following lemma captures useful properties of Rigidize. See the appendix for a proof.

\begin{lemma}\label{lem:rigidize}
  The algorithm Rigidize runs in time $h^{O(s)}$ where $h=\abs{V(H)}$, and the output set $R\subseteq V(H)$ has the properties that $\abs{R}\le O(s^2)$, that $H-R$ is a matching, and that every $f \in \ensuremath{{\mathrm{Aut}}}(H)$ satisfies $f(R)=R$.
\end{lemma}

\subsection{Counting subgraphs with rigid splitting sets}\label{sub: fixed split}

We use the rigid splitting set~$R$ from \cref{lem:rigidize} to
compute the number of times~$H$ occurs as a subgraph modulo a power of two.
As a subroutine, we use an algorithm for counting colored matchings modulo a power of two in a setting involving particular ``color demands''.
\begin{definition}
Let~$G$ be a graph, let~$C$ be a finite set of colors, and let $c\colon V(G)\cup E(G)\to 2^C$ be a function that labels each vertex and edge with a subset of~$C$.
For any matching~$M$, let $I(M)$ be the set of its isolated vertices.
  For a coloring $c_M:I(M)\cup E(M)\to C$, the colored matching $(M,c_M)$ is \emph{permissible} if $c_M(t)\in c(t)$ holds for all~$t\in I(M)\cup E(M)$.

\emph{Color demands} are functions~$D_I,D_E\colon C \to \mathbf{N}$.
The pair $(M,c_M)$ \emph{satisfies the demands~$D_I,D_E$} if, for each~$i\in C$, the graph~$M$ contains exactly $D_I(i)$ isolated vertices~$v$ with $c_M(v)=i$ and exactly $D_E(i)$ edges with $c_M(e)=i$.
Let $\mathcal M(G,c,D_I,D_E)$ be the set of all permissible matchings~$(M,c_M)$ that satisfy the demand~$D$.
\end{definition}

As shown in the appendix, we obtain the following algorithm as a corollary to Hirai and Namba's algorithm~\cite{HiraiNamba} for computing the Hafnian over polynomial rings modulo~$2^t$.
\begin{lemma}\label{lem:color-demands}
  Given a graph~$G$, permissible colors~$c\colon V(G)\cup E(G)\to 2^C$, color demands $D_I,D_E\colon C\to\mathbf{N}$, and $t\in\mathbf{N}_{\ge1}$, there is an algorithm that computes the number $\abs{\mathcal M(G,c,D_I,D_E)}\bmod 2^t$ in time $n^{O(t\abs{C})}$.
\end{lemma}

Before we state the main algorithm, we introduce some basic group-theoretic notation.
Let~$R$ be a splitting set of~$H$ that satisfies $f(R)=R$ for all $f\in\ensuremath{{\mathrm{Aut}}}(H)$.
Let $G$ be a graph and let $S\subseteq V(G)$ be a set with $\abs{S}=\abs{R}$.
For an embedding $\sigma\in\ensuremath{{\mathrm{Emb}}}(H[R],G[S])$ and an
automorphism~$\varphi\in\ensuremath{{\mathrm{Aut}}}(H)$, we note that the function
$\sigma\circ(\varphi\restrict[R])$ is again an embedding in $\ensuremath{{\mathrm{Emb}}}(H[R],G[S])$.
Indeed, we view this operation as a right-action of the group $\ensuremath{{\mathrm{Aut}}}(H)$ on the
set $\ensuremath{{\mathrm{Emb}}}(H[R],G[S])$.
We call two embeddings $\sigma,\sigma'\in\ensuremath{{\mathrm{Emb}}}(H[R],G[S])$ \emph{equivalent} if
there exists~$\varphi\in\ensuremath{{\mathrm{Aut}}}(H)$ such that
$\sigma'=\sigma\circ(\varphi\restrict[R])$; this clearly defines an equivalence
relation.
The equivalence class $\sigma\ensuremath{{\mathrm{Aut}}}(H)$ is called the \emph{orbit} of $\sigma$
under $\ensuremath{{\mathrm{Aut}}}(H)$.
All orbits have the same size.
Let $E_S$ be a set of representatives for each orbit, that is, a maximal set of
mutually non-equivalent embeddings in~$\ensuremath{{\mathrm{Emb}}}(H[R],G[S])$.

We are ready to state the modular counting algorithm for $s$-matching-splittable subgraphs.

\begin{algorithm}{ModCount}{$H,G,t$}{Given an $s$-matching-splittable graph~$H$, a
    host graph~$G$, and an integer~$t\ge 2$, this algorithm computes the number
    $\#\ensuremath{{\mathrm{Sub}}}(H,G)\bmod 2^t$.}
\item[C1] (Compute rigid splitting set.)
  Call Rigidize$(H)$ to compute the set~$R$.
\item[C2] (Reduce to counting colored matchings.)
  For each $S\subseteq V(G)$ with $\abs{S}=\abs{R}$ (that is, a possible image of $R$)
  and each representative embedding $\sigma\in E_S$ from $H[R]$ to $G[S]$,
  we construct an instance~$(G-S,c_\sigma,D_I,D_E)$ of colored matching with demands and
  use \cref{lem:color-demands} to obtain the number $\abs{\mathcal M(G-S,c_\sigma,D_I,D_E)}\bmod 2^t$:
  \begin{description}
    \item[a.]
      (Set permitted colors.)
      Let~$C=2^R\cup \binom{2^R}{1}\cup\binom{2^R}{2}$.
      For each vertex ${v\in V(G)\setminus S}$,
      let $N_v\subseteq R$ be the vertices of $R$ that hit $N_G(v)\cap S$ under $\sigma$, that is,
      $N_v=\sigma^{-1}(N_G(v)\cap S)$.
      Define~$c_\sigma(v)=\setc{N}{N\subseteq N_v}$.
      Moreover, for each $\set{u,v}\in E(G-S)$, define $c_\sigma(\set{u,v})=\setc{\set{N,N'}}{N\subseteq N_u,\; N'\subseteq N_v}$.
    \item[b.]
      (Make demands.)
      The demands $D_I,D_E:C\to\mathbf{N}$ depend only on~$H$ and~$R$.
      For each $N\subseteq R$, we let
      $D_I(N)$ be the number of isolated vertices~$v$ in $H-R$ whose neighborhood in~$H$
      satisfies $N_H(v)\cap R=N$.
      Moreover, for all $N,N'\subseteq R$, we let
      $D(\set{N,N'})$ be the number of edges~$\set{u,v}\in E(H-R)$ with $\set{N_H(u)\cap R, N_H(v)\cap R}=\set{N,N'}$.
  \end{description}
\item[C3] (Sum up.)
  Output the sum modulo $2^t$ of all integers returned by the queries in C2.
\end{algorithm}

  In the appendix, we prove that ModCount satisfies the properties stated in \cref{thm: algorithm-split}.

\section{Hardness of counting paths modulo two}\label{sec: hardness-paths}

In this section, we prove \cref{thm: hardness-paths}, that counting $k$-paths modulo $2$ is $\oplus\W$-hard.
We first formally introduce this and some intermediate problems.

The \emph{length} of a path is the number of its edges.
For a graph~$G$ and vertices~$s,t\in V(G)$, an \emph{$s,t$-path} is a simple path from~$s$ to~$t$.
For a computable, strictly increasing function~$f\colon\mathbf{N}\to\mathbf{N}$, we define \pp{$f$-Flexible Path} to be the problem that is given~$(G,s,t,k)$ to decide whether there exists any~$s,t$-path in~$G$ whose length~$\ell$ satisfies $k\le\ell\le f(k)$.
When $\ensuremath{{\mathrm{id}}}$ denotes the identity function, then \pp{Path} (also known as \pp{$k$-Path} or \textsc{Longest Path}) is defined as \pp{$\ensuremath{{\mathrm{id}}}$-Flexible Path}.
We similarly define \pp{Directed $f$-Flexible Path} and \pp{Directed Path} for directed graphs, and we define the counting and parity versions of these problems in the canonical manner.

We start our reduction at the vertex-colorful subgraph problem $\oplus\ensuremath{{\mathrm{VertexColorfulSub}}}(\mathcal H)$ for a class $\mathcal H$ of unbounded tree-width, which is $\oplus\W$-hard by \cref{cor: hardness-tw}.
The core part of the reduction is in \cref{lem:vcsubs-to-dipaths}, where we reduce to counting paths of somewhat flexible length in a directed graph (modulo~$2$).
From there, we reduce to the familiar $k$-path problem in undirected graphs using standard tricks.

\begin{lemma}\label{lem:vcsubs-to-dipaths}%
  For any class $\mathcal H$ of connected $4$-regular graphs without non-trivial automorphisms,
  there is a computable strictly increasing function~$f$
  such that $\oplus\ensuremath{{\mathrm{VertexColorfulSub}}}(\mathcal H)$ admits a parsimonious polynomial-time fpt-reduction to \pp{$\oplus$Directed $f$-Flexible Path}.
\end{lemma}

\begin{proof}
  Let $H\in\mathcal H$ and $G$ be given as input,
  where $G$ is given with color classes $V_u$ for~$u\in V(H)$.
  Let $k=\abs{E(H)}$.
  Since $H$ is connected and $4$-regular, it has an Eulerian tour
  $u_0,u_1,\dots,u_{k-1},u_{k}$ such that $E(H)=\setc{ \set{u_i,u_{i+1}} }{ i\in\set{0,\dots,k-1} }$.
  Every vertex of~$H$ appears exactly twice on the tour, except
  for $u_0=u_{k}$, which appears three times.
  We construct a graph $G'$ that visits every color class of~$G$ two (or three)
  times according to the Eulerian tour in~$H$.

  Before we give a formal construction, we give an overview.
  The graph~$G'$ will essentially be a sequence of bipartite
  graphs~$B_0,\dots,B_\ell$ whose edges are all directed from left to right.
  For a bipartite graph $B$, we write $L(B)$ and $R(B)$ for
  its left and right part, respectively.
  We will have $R(B_j)=L(B_{j+1})$ for all $j\in\set{0,\dots,\ell}$.
  Each $B_j$ will either be a perfect matching~$M_i$ or a graph $G_i$ that is a
  copy of $G[V_{u_{i-1}}\cup V_{u_{i}}]$. (Note that $G_i$ is indeed bipartite, since $G$ is $H$-colored and $H$ contains no self-loops.)
  Pictorially, the sequence of bipartite graphs is $M_0 G_1 M_1 \dots M_{k-1}
  G_k M_k M_{k+1}$.
  We also add some additional gadget edges between non-neighboring
  layers.

  We now describe the construction of $G'$ in detail:
  \begin{enumerate}
    \item
      \textbf{Graph edges.}
      For each $i\in\set{1,\dots,k}$, let $G_i$ be a fresh copy of
      $G[V_{u_{i-1}} \cup V_{u_{i}}]$, renamed so that $L(G_i)=\set{i}\times
      V_{u_{i-1}}$ and $R(G_i)=\set{i}\times V_{u_i}$.

    \item
      \textbf{Matching edges.}
      For each $i\in\set{1,\dots,k-1}$, we connect $G_i$ with $G_{i+1}$ by
      adding the canonical perfect matching~$M_i$ from $L(M_i)\doteq
      R(G_i)=\set{i}\times V_{u_i}$ to $R(M_i)\doteq L(G_{i+1})=\set{i+1}\times
      V_{u_i}$.

      At the last graph~$G_{k}$, we add a dangling matching $M_{k}$ with
      $L(M_k)=R(G_k)=\set{k}\times V_{u_k}$ and $R(M_k)=\set{k+1}\times
      V_{u_k}$.
      We add a further dangling matching $M_{k+1}$ with $L(M_{k+1})=R(M_k)$ and
      $R(M_{k+1})=\set{k+1}\times V_{u_k}$.
      And a final dangling matching~$M_0$ to the left of $G_1$ so that
      $R(M_0)=L(G_1)$ and $L(M_0)=\set{0}\times V_{u_0}$.

    \item
      \textbf{Gadget edges.}
      For all $i,j\in\set{1,\dots,k}$ with $i < j$ and $u_i = u_j$, note that
      $L(M_i)=\set{i}\times V_{u_i}$ and $L(M_j)=\set{j}\times V_{u_i}$.
      We add the canonical \emph{bidirected} matching between $L(M_i)$ and
      $L(M_j)$.
      Similarly, we add the canonical bidirected matching between $R(M_i)$ and
      $R(M_j)$.
      We also add those matchings between $L(M_0)$ and $L(M_{k+1})$ and between
      $R(M_0)$ and $R(M_{k+1})$.

    \item
      \textbf{Source/sink.}
      Let $s$ be a new vertex and add all edges $(s,v)$ for $v\in L(M_0)$.
      Let $t$ be a new vertex and add all edges $(v',t)$ for $v'\in R(M_{k+1})$.

    \item
      \textbf{Parameters.}
      Finally, we set $k' = 2k + 4$ and $f(k') = k + 3(k+2) + 2$ so that we are
      counting all $s,t$-paths whose length is between $k'$ and $f(k')$.
  \end{enumerate}
  This finishes the construction of~$G'$.
  Note that $G'\setminus\set{s,t}$ is indeed a sequence $B_0\dots B_\ell$ of
  bipartite graphs with some additional gadget edges, where $\ell=2k+2$.
  We define the \emph{$j$-th layer of $G'$} as the set $L_j=L(B_j)$ for
  $j\le\ell$ and $L_{\ell+1}=R(B_\ell)$ for $j=\ell+1$.
  Recall that $L_j=R(B_{j-1})$ for $j>0$.

  \medskip
  We describe the \emph{canonical} solutions in the output of the reduction.

  Thus, let $H'$ be an $H$-subgraph of $G$ that is colorful.
  Recall that this means that $H$ is isomorphic to $H'$ and that the `coloring' homomorphism $c\colon V(G)\rightarrow V(H)$ is bijective on $H'$.
  Since $H'$ is finite, $c$ is in fact an isomorphism.
  (It must map non-edges to non-edges because $E(H)$ and $E(H')$ have the same size.)
  Let $\phi\colon H'\rightarrow H$ be the inverse of $c$ and note that $\phi$ is a homomorphism.

  We define the canonical $s,t$-path $p=p_{H'}$ in $G'$ corresponding to $H'$ by picking one
  vertex from each layer.
  More precisely, $p$ is the unique $s,t$-path that does not use any gadget
  edges and that satisfies the following property for all
  $j\in\set{0,\dots,\ell+1}$: Setting $i$ and $u$ such that $L_j=\set{i}\times
  V_u$, we have $V(p) = \set{(i,\phi(u))}$.
  Note that this determines a unique vertex set of~$p$.
  We claim that $V(p)$ indeed induces a path on the graph and matching edges,
  which must be unique if it exists.

  To show that $p$ is a path, let $j\in\{0,\ldots, \ell\}$.
  We claim that the two vertices in $V(p)\cap V(B_j)$ are adjacent in $B_j$.
  If $B_j$ is one of the matching graphs, then $U_j=\set{i}\times V_u$ for some
  $i$ and $u$, and $U_{j+1}=\set{i+1}\times V_u$.
  By choice of the matching, there is indeed an edge from $(i,\phi(u))$ to
  $(i+1,\phi(u))$ in~$B_j$.
  Otherwise, $B_j$ is one of the graph copies, say $L_j=\set{i}\times
  V_u$ and $L_{j+1}=\set{i}\times V_{v}$; note that by construction of $G'$, we
  only do this if $\{u,v\}$ is an edge of $H$ since $v$ is the successor of $u$ in an
  Eulerian tour of $H$.
  Since $\phi$ is a graph homomorphism respecting the coloring, we have that
  $\phi(u) \phi(v)$ is an edge in $G[V_u\cup V_v]$.
  Since $B_j$ was a copy of this graph by construction, there is an edge from
  $(i,\phi(u))$ to $(i,\phi(v))$.
  Overall, we get that $spt$ is an $s,t$-path in $G'$, and its length is
  $\ell+2=2k+4=k'$; this is the canonical path corresponding to~$\phi$.

  In summary, every vertex-colorful $H$-subgraph $H'$ defines a unique canonical $s,t$-path.
  Conversely, let $p$ be a path whose internal vertex set is of the form
  $V(p) = \set{(i_j,x_j)}$ where $x_j\in V_u$ for layer  $L_j=\set{i}\times V_u$, and that is consistent in the sense that it always picks the same vertex from each color class.
  (Formally, if $u_i=u_j$ then $x_i=x_j$.)
  Then $p$ describes a set of $|V(H)|$ many vertices and $k$ edges in $G$ that make up a colorful $H$-subgraph $H'$ of $G$.
  Since $H$ has no nontrivial automorphisms, this graph is uniquely defined.

  \medskip

  Let $\mathcal P$ be the set of all $s,t$-paths whose length~$\ell$ satisfies $k'\le\ell\le f(k')$.
  The central claim is that the number of paths in~$\mathcal P$ that are not canonical solutions is even.
  For this, we construct a fixed-point free involution~$\pi$ on non-canonical
  solutions.

  First we concentrate on solutions that are \emph{jumpy}:
  A path $p\in\mathcal P$ is jumpy if
  there exist $i,j\in\set{1,\dots,k}$ with $i\ne j$ so that we added gadget
  edges between $M_i$ and $M_j$, and
  \begin{enumerate}[(i)]
    \item
      $p$ uses the edge from $(i,x)\in L(M_i)$ to $(j,x)\in L(M_j)$ but not the
      edge from $(j+1,x)\in R(M_j)$ to $(i+1,x)\in R(M_i)$,
      or
    \item
      $p$ uses the edge from $(i+1,x)\in R(M_i)$ to $(j+1,x)\in R(M_j)$ but not
      the edge from $(j,x)\in L(M_j)$ and $(i,x)\in L(M_i)$.
  \end{enumerate}
  If $p$ is jumpy, we define $\pi(p)$ as follows: First we identify the
  lexicographically first position where $p$ is jumpy.
  Then we exchange state (i) with state (ii) at that position.
  Note that (i) implies that $p$ uses the edge from $(j,x)$ to $(j+1,x)$ in
  $M_j$ and (ii) implies that $p$ uses the edge from $(i,x)$ to $(i+1,x)$ in
  $M_i$; we swap these edges too, when applying $\pi$.
  Now $\pi$ is a fixed-point free involution on jumpy paths, and note that
  $\pi(p)$ has the same length as $p$.

  Next we consider paths that are \emph{bad}.
  A path $p\in \mathcal P$ is bad if it is not jumpy, and yet there exist
  $i,j\in\set{1,\dots,k}$ with $i\ne j$ so that we added gadget edges between
  $M_i$ and $M_j$, and
  \begin{enumerate}[(i)]
    \item
      $p$ uses the $M_i$-edge from $(i,x)$ to $(i+1,x)$ and not the $M_j$-edge
      from $(j,x)$ to $(j+1,x)$, or
    \item
      $p$ uses the three edges from $(i,x)$ to $(j,x)$ to $(j+1,x)$ to
      $(i+1,x)$.
  \end{enumerate}
  We define $\pi(p)$ for a bad path~$p$ by finding the first bad position, and
  switching between (i) and (ii).
  Say, $p$ was in state (i) at that position, then $\pi(p)$ is in state (ii) at
  that position, which makes the overall path exactly two edges longer.
  We claim that all non-jumpy paths pair up in this manner, without having to
  consider arbitrarily long paths.
  Indeed, since $p$ is not jumpy, when we look at the vertices that~$p$
  traverses, we are essentially traversing the layers in a monotone order, except for potential short two-vertex detours in the case of gadgets in
  a bad state~(ii).
  More precisely, when $p$ reaches a vertex in a layer~$L_j$, either the next
  vertex is in $L_{j+1}$ or the third vertex after it is in $L_{j+1}$.
  This means that the paths moves to the right by one layer at least once ever
  $4$ vertices, and thus the longest path that is not jumpy has length at most
  $k+3(k+2)+2=k'$, accounting for $k$ graph edges, at most~$3$ matching/gadget edges
  for each of the $k+2$ gadgets, and two edges at the source and sink.
\end{proof}

For completeness, we include two simple reductions:
First from the flexible-length to the fixed-length problem in directed graphs, then from the directed to undirected problem.
\begin{lemma}\label{lem:dipath-flexibility}
  Let $f:\mathbf{N}\to\mathbf{N}$ be a computable, strictly increasing function.
  There is a parsimonious polynomial-time fpt-reduction from \pp{\#Directed $f$-Flexible Paths} to \pp{\#Directed
  Paths}.
\end{lemma}

\begin{lemma}\label{lem:dipath-to-path}
  There is a parsimonious polynomial-time fpt-reduction from \pp{\#Directed Paths} to \pp{\#Paths}.
\end{lemma}

With all these prerequisites collected, we can complete the proof.

\begin{proof}[Proof of \cref{thm: hardness-paths}]
  Let $\mathcal H$ be the class of all graphs~$H$ that are connected, $4$-regular, and whose automorphism group has size one.
  We use the probabilistic method to argue that the tree-width of graphs in~$\mathcal H$ is not bounded.
  With probability $1-o(1)$ as $h\to\infty$, random $4$-regular graphs with~$h$ vertices are connected~\cite[Theorem~2.10]{wormald1999models},
  they have no nontrivial automorphisms~\cite{DBLP:journals/rsa/KimSV02},
  and they are almost Ramanujan~\cite[Theorem~7.10]{hoory2006expander}, that is, their second-largest eigenvalue in absolute value satisfies $\lambda\le 2\sqrt{3}+o(1)<3.5$.
  By a union bound, $H$ has all three properties simultaneously with probability $1-o(1)$.
  By Cheeger's inequality~\cite[Theorem 4.11]{hoory2006expander}, we have
  \[ \min_{S\subseteq V(H), |S|\le\frac12 h } \frac{\abs{E(S,\overline S)}}{\abs{S}}\ge \frac{1}{2}(4-\lambda)>0.1\,,\] that is, the edge expansion is bounded away from zero, which implies that the tree-width of~$H$ is at least linear in~$h$ (see, e.g.,~\cite[Exercise~7.34]{DBLP:books/sp/CyganFKLMPPS15}).
  Thus, $\mathcal H$ contains graphs of arbitrarily large tree-width.

  By \cref{cor: hardness-tw}, $\oplus\ensuremath{{\mathrm{VertexColorfulSub}}}(\mathcal H)$ is $\oplus\W$-hard under parsimonious fpt-reductions.
  If there is a parsimonious fpt-reduction from problem $\#A$ to problem $\#B$ then in particular the parity version $\oplus A$ reduces to $\oplus B$.
  Writing $\oplus A\leq \oplus B$ we can summarize the chain of reductions in \cref{lem:vcsubs-to-dipaths,lem:dipath-flexibility,lem:dipath-to-path} as
  \[
    \oplus\ensuremath{{\mathrm{VertexColorfulSub}}}(\mathcal H)
    \le
    \oplus\pp{Directed $f$-Flexible Paths}\\
   \le
    \oplus\pp{Directed Paths}
    \le
    \oplus\pp{Paths}\,.
  \]
  This proves the $\oplus\W$-hardness of $\oplus\pp{Paths}$.
  The containment follows from the standard fpt-reduction from $\#\pp{Paths}$ to $\#\pp{Clique}$, which is parsimonious.
  Overall, the claim follows.
\end{proof}

\section{Conclusion}

We conducted an initial investigation of modular subgraph counting, leading to the partial classification depicted in \cref{fig: overview}. To obtain a complete picture, the following conjecture needs to be addressed.

\begin{conjecture}\label{conjecture}
  For any computable pattern class $\mathcal H$:
  \begin{itemize}
    \item If $\mathcal H$ has unbounded matching-split number, then the problem $\oplus\ensuremath{{\mathrm{Sub}}}(\mathcal H)$ is $\oplus\W$-complete.
    \item If $\mathcal H$ has unbounded vertex-cover number, then $\ModSub{\mathcal H}{q}$ for fixed $q \in \mathbf{N}$ is $\ensuremath{{\mathrm{Mod}}}_p \cc{W[1]}$-complete for any odd divisor $p$ of $q$.
  \end{itemize}
\end{conjecture}

An appropriate transfer of the subgraph-homomorphism framework to modular counting is likely to help in settling this conjecture.

\bibliographystyle{plainurl}


\section{Introduction}
Likelihood-based generative modeling is a central task in machine learning that is the basis for a wide range of applications ranging from speech synthesis~\citep{oord2016wavenet}, to translation~\citep{sutskever2014sequence}, to compression~\citep{mackay2003information}, to many others. Autoregressive models have long been the dominant model class on this task due to their tractable likelihood and expressivity, as shown in Figure~\ref{fig:sota}. Diffusion models have recently shown impressive results in image \citep{ho2020denoising, song2020score,nichol2021improved} and audio generation \citep{kong2020diffwave, chen2020wavegrad} in terms of perceptual quality, but have yet to match autoregressive models on density estimation benchmarks. In this paper we make several technical contributions that allow diffusion models to challenge the dominance of autoregressive models in this domain.
\begin{figure}[bht]
\begin{subfigure}{.485\textwidth}
  \centering
  %
  \includegraphics[width=.99\linewidth]{media/images/cifar10_sota_over_time}
  \caption{CIFAR-10 without data augmentation}
\end{subfigure}
\begin{subfigure}{.495\textwidth}
  \centering
  %
  \includegraphics[width=.99\linewidth]{media/images/imagenet64_sota_over_time}
  \caption{ImageNet 64x64}
\end{subfigure}
\caption{Autoregressive generative models were long dominant in standard image density estimation benchmarks. In contrast, we propose a family of diffusion-based generative models, \emph{Variational Diffusion Models} (VDMs), that outperforms contemporary autoregressive models in these benchmarks. See Table \ref{table:results} for more results and comparisons.}
\label{fig:sota}
\end{figure}
Our main contributions are as follows:
\begin{itemize}
\item We introduce a flexible family of diffusion-based generative models that achieve new state-of-the-art log-likelihoods on standard image density estimation benchmarks (CIFAR-10 and ImageNet). This is enabled by incorporating Fourier features into the diffusion model and using a learnable specification of the diffusion process, among other modeling innovations.
\item We improve our theoretical understanding of density modeling using diffusion models by analyzing their variational lower bound (VLB), deriving a remarkably simple expression in terms of the signal-to-noise ratio of the diffusion process. This result delivers new insight into the model class: for the continuous-time (infinite-depth) setting we prove a novel invariance of the generative model and its VLB to the specification of the diffusion process, and we show that various diffusion models from the literature are equivalent up to a trivial time-dependent rescaling of the data.
\end{itemize}

\section{Related work}
\label{sec:relatedwork}

Our work builds on diffusion probabilistic models (DPMs)~\citep{sohl2015deep}, or \emph{diffusion models} in short. DPMs can be viewed as a type of variational autoencoder (VAE)~\citep{kingma2013auto,rezende2014stochastic}, whose structure and loss function allows for efficient training of arbitrarily deep models. Interest in diffusion models has recently reignited due to their impressive image generation results \citep{ho2020denoising, song2020improved}.

\citet{ho2020denoising} introduced a number of model innovations to the original DPM, with impressive results on image generation quality benchmarks. They showed that the VLB objective, for a diffusion model with discrete time and diffusion variances shared across input dimensions, is equivalent to multi-scale denoising score matching, up to particular weightings per noise scale. Further improvements were proposed by \citet{nichol2021improved}, resulting in better log-likelihood scores. \citet{gao2020learning} show how diffusion can also be used to efficiently optimize energy-based models (EBMs) towards a close approximation of the log-likelihood objective, resulting in high-fidelity samples even after long MCMC chains.

\citet{song2019generative} first proposed learning generative models through a multi-scale denoising score matching objective, with improved methods in \citet{song2020improved}. This was later extended to continuous-time diffusion with novel sampling algorithms based on reversing the diffusion process~\citep{song2020score}.

Concurrent to our work, \cite{song2021maximum}, \cite{huang2021variational}, and \cite{vahdat2021score} also derived variational lower bounds to the data likelihood under a continuous-time diffusion model. Where we consider the infinitely deep limit of a standard VAE, \cite{song2021maximum} and \cite{vahdat2021score} present different derivations based on stochastic differential equations. \cite{huang2021variational} considers both perspectives and discusses the similarities between the two approaches. An advantage of our analysis compared to these other works is that we present an intuitive expression of the VLB in terms of the signal-to-noise ratio of the diffused data, leading to much simplified expressions of the discrete-time and continuous-time loss, allowing for simple and numerically stable implementation. This also leads to new results on the invariance of the generative model and its VLB to the specification of the diffusion process. We empirically compare to these works, as well as others, in Table~\ref{table:results}.

Previous approaches to diffusion probabilistic models fixed the diffusion process; in contrast optimize the diffusion process parameters jointly with the rest of the model. This turns the model into a type of VAE~\citep{kingma2013auto,rezende2014stochastic}. This is enabled by directly parameterizing the mean and variance of the marginal $q({\mathbf{z}}_t|{\mathbf{z}}_0)$, where previous approaches instead parameterized the individual diffusion steps $q({\mathbf{z}}_{t+\epsilon}|{\mathbf{z}}_t)$. In addition, our denoising models include several architecture changes, the most important of which is the use of Fourier features, which enable us to reach much better likelihoods than previous diffusion probabilistic models.

\section{Model}
\label{sec:model}
We will focus on the most basic case of generative modeling, where we have a dataset of observations of ${\mathbf{x}}$, and the task is to estimate the marginal distribution $p({\mathbf{x}})$. As with most generative models, the described methods can be extended to the case of multiple observed variables, and/or the task of estimating conditional densities $p({\mathbf{x}}|{\mathbf{y}})$. The proposed latent-variable model consists of a diffusion process (Section ~\ref{sec:inferencemodel}) that we invert to obtain a hierarchical generative model (Section ~\ref{sec:generativemodel}). As we will show, the model choices below result in a surprisingly simple variational lower bound (VLB) of the marginal likelihood, which we use for optimization of the parameters.

\subsection{Forward time diffusion process}\label{sec:inferencemodel}
Our starting point is a Gaussian diffusion process that begins with the data ${\mathbf{x}}$, and defines a sequence of increasingly noisy versions of ${\mathbf{x}}$ which we call the \emph{latent variables} ${\mathbf{z}}_t$, where $t$ runs from $t=0$ (least noisy) to $t=1$ (most noisy).
The distribution of latent variable ${\mathbf{z}}_t$ conditioned on ${\mathbf{x}}$, for any $t \in [0,1]$ is given by:
\begin{align}
q({\mathbf{z}}_t|{\mathbf{x}}) &= \mathcal{N}\left(\alpha_t {\mathbf{x}}, \sigma^2_t \mathbf{I}\right),
\label{eq:zt_given_x}
\end{align}
where  $\alpha_t$ and $\sigma^{2}_t$ are strictly positive scalar-valued functions of $t$. Furthermore, let us define the \emph{signal-to-noise ratio} (SNR): \begin{align}
    \text{SNR}(t) = \alpha^2_t / \sigma^2_t.
\end{align}
We assume that the $\text{SNR}(t)$ is strictly monotonically decreasing in time, i.e. that $\text{SNR}(t) < \text{SNR}(s)$ for any $t > s$. This formalizes the notion that the ${\mathbf{z}}_t$ is increasingly noisy as we go forward in time. We also assume that both $\alpha_{t}$ and $\sigma^{2}_t$ are smooth, such that their derivatives with respect to time $t$ are finite.  This diffusion process specification includes the \emph{variance-preserving} diffusion process as used by ~\citep{sohl2015deep,ho2020denoising} as a special case, where $\alpha_t = \sqrt{1 - \sigma^2_t}$. Another special case is the \emph{variance-exploding} diffusion process as used by \citep{song2019generative,song2020score}, where $\alpha^2_t = 1$. In experiments, we use the variance-preserving version.

The distributions $q({\mathbf{z}}_t|{\mathbf{z}}_s)$ for any $t > s$ are also Gaussian, and given in Appendix \ref{sec:distributions}. The joint distribution of latent variables $({\mathbf{z}}_s, {\mathbf{z}}_t, {\mathbf{z}}_u)$ at any subsequent timesteps $0 \leq s < t < u \leq 1$ is Markov: $q({\mathbf{z}}_u | {\mathbf{z}}_t, {\mathbf{z}}_s) = q({\mathbf{z}}_u | {\mathbf{z}}_t)$. Given the distributions above, it is relatively straightforward to verify through Bayes rule that $q({\mathbf{z}}_s|{\mathbf{z}}_t,{\mathbf{x}})$, for any $0 \leq s < t \leq 1$, is also Gaussian. This distribution is also given in Appendix \ref{sec:distributions}.

\subsection{Noise schedule}
\label{sec:learning_noise}

In previous work, the noise schedule has a fixed form (see \cref{app:monotonicnn}, \cref{fig:snr_schedule}). In contrast, we learn this schedule through the parameterization
\begin{align}
\sigma^{2}_t = \text{sigmoid}(\gamma_{{\boldsymbol{\eta}}}(t))
\end{align}
where $\gamma_{{\boldsymbol{\eta}}}(t)$ is a monotonic neural network with parameters ${{\boldsymbol{\eta}}}$, as detailed in Appendix \ref{app:monotonicnn}.

Motivated by the equivalence discussed in Section \ref{sec:schedule_invariance}, we use $\alpha_t = \sqrt{1 - \sigma^2_t}$ in our experiments for both the discrete-time and continuous-time models, i.e. variance-preserving diffusion processes. It is straightforward to verify that $\alpha^2_t$ and $\text{SNR}(t)$, as a function of $\gamma_{{\boldsymbol{\eta}}}(t)$, then simplify to:
\begin{align}
\alpha^{2}_t = \text{sigmoid}(-\gamma_{{\boldsymbol{\eta}}}(t)) \\
\text{SNR}(t) = \exp(-\gamma_{{\boldsymbol{\eta}}}(t))
\end{align}

\subsection{Reverse time generative model}\label{sec:generativemodel}
We define our generative model by inverting the diffusion process of Section~\ref{sec:inferencemodel}, yielding a hierarchical generative model that samples a sequence of latents ${\mathbf{z}}_t$, with time running backward from $t=1$ to $t=0$.
We consider both the case where this sequence consists of a finite number of steps $T$, as well as a continuous time model corresponding to $T\rightarrow\infty$. We start by presenting the discrete-time case.

Given finite $T$, we discretize time uniformly into $T$ timesteps (segments) of width $\tau = 1/T$. Defining $s(i) = (i-1)/T$ and $t(i) = i/T$, our hierarchical generative model for data ${\mathbf{x}}$ is then given by:
\begin{align}
p({\mathbf{x}}) &= \int_{{\mathbf{z}}} p({\mathbf{z}}_{1})p({\mathbf{x}}|{\mathbf{z}}_{0}) \prod_{i=1}^T p({\mathbf{z}}_{s(i)}|{\mathbf{z}}_{t(i)}).
\end{align}
With the variance preserving diffusion specification and sufficiently small $\text{SNR}(1)$, we have that $q({\mathbf{z}}_1|{\mathbf{x}}) \approx \mathcal{N}({\mathbf{z}}_1; 0,\mathbf{I})$. We therefore model the marginal distribution of ${\mathbf{z}}_1$ as a spherical Gaussian:
\begin{align}
    p({\mathbf{z}}_1) = \mathcal{N}({\mathbf{z}}_1; 0, \mathbf{I}).
\end{align}
We wish to choose a model $p({\mathbf{x}}|{\mathbf{z}}_0)$ that is close to the unknown $q({\mathbf{x}}|{\mathbf{z}}_0)$. Let $x_i$ and $z_{0,i}$ be the $i$-th elements of ${\mathbf{x}},  {\mathbf{z}}_0$, respectively. We then use a factorized distribution of the form:
\begin{align}
p({\mathbf{x}}|{\mathbf{z}}_0) = \prod_i p(x_i|z_{0,i})\label{eq:recon},
\end{align}
where we choose $p(x_i |z_{0,i}) \propto q(z_{0,i}|x_i)$, which is normalized by summing over all possible discrete values of $x_i$ (256 in the case of 8-bit image data). With sufficiently large $\text{SNR}(0)$, this becomes a very close approximation to the true $q({\mathbf{x}}|{\mathbf{z}}_0)$, as the influence of the unknown data distribution $q({\mathbf{x}})$ is overwhelmed by the likelihood $q({\mathbf{z}}_0|{\mathbf{x}})$.
Finally, we choose the conditional model distributions as \begin{align}p({\mathbf{z}}_s|{\mathbf{z}}_t) = q({\mathbf{z}}_s|{\mathbf{z}}_t, {\mathbf{x}}=\hat{{\mathbf{x}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}_t; t)),
\end{align}
i.e.\ the same as $q({\mathbf{z}}_s|{\mathbf{z}}_t,\mathbf{x})$, but with the original data ${\mathbf{x}}$ replaced by the output of a \emph{denoising model} $\hat{{\mathbf{x}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}_t; t)$ that predicts ${\mathbf{x}}$ from its noisy version ${\mathbf{z}}_t$. Note that in practice we parameterize the denoising model as a function of a \emph{noise prediction model} (Section \ref{sec:noisepred}), bridging the gap with previous work on diffusion models \citep{ho2020denoising}. The means and variances of $p({\mathbf{z}}_s|{\mathbf{z}}_t)$ simplify to a remarkable degree; see Appendix \ref{sec:distributions}.

\subsection{Noise prediction model and Fourier features}
\label{sec:noisepred}

We parameterize the denoising model in terms of a \emph{noise prediction model} $\hat{{\boldsymbol{\epsilon}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}_t; t)$:
\begin{align}
    \hat{{\mathbf{x}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}_t; t) = ({\mathbf{z}}_t - \sigma_t \hat{{\boldsymbol{\epsilon}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}_t; t))/\alpha_t,
\label{eq:noisepred}
\end{align}
where $\hat{{\boldsymbol{\epsilon}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}_t; t)$ is parameterized as a neural network. The noise prediction models we use in experiments closely follow~\cite{ho2020denoising}, except that they process the data solely at the original resolution. The exact parameterization of the noise prediction model and noise schedule is discussed in Appendix \ref{sec:hps}.

Prior work on diffusion models has mainly focused on the perceptual quality of generated samples, which emphasizes coarse scale patterns and global consistency of generated images. Here, we optimize for likelihood, which is sensitive to fine scale details and exact values of individual pixels.  To capture the fine scale details of the data, we propose adding a set of \emph{Fourier features} to the input of our noise prediction model. Let ${\mathbf{x}}$ be the original data, scaled to the range $[-1,1]$, and let ${\mathbf{z}}$ be the resulting latent variable, with similar magnitudes. We then append channels $\sin(2^n\pi{\mathbf{z}})$ and $\cos(2^n\pi{\mathbf{z}})$, where $n$ runs over a range of integers $\{n_{min}, ..., n_{max}\}$. These features are high frequency periodic functions that amplify small changes in the input data~${\mathbf{z}}_t$; see Appendix \ref{sec:fourier} for further details.
Including these features in the input of our denoising model leads to large improvements in likelihood as demonstrated in Section~\ref{sec:experiments} and Figure \ref{fig:fourier-ablation}, especially when combined with a learnable SNR function. We did not observe such improvements when incorporating Fourier features into autoregressive models.

\subsection{Variational lower bound}
We optimize the parameters towards the variational lower bound (VLB) of the marginal likelihood, which is given by
\begin{align}
    - \log p({\mathbf{x}}) \leq %
    - \text{VLB}({\mathbf{x}}) =
    \underbrace{%
    D_{KL}(q({\mathbf{z}}_1|{\mathbf{x}})||p({\mathbf{z}}_1))
    %
    }_{\text{Prior loss}}
    + \underbrace{%
    \mathbb{E}_{q({\mathbf{z}}_0|{\mathbf{x}})}\left[ - \log p({\mathbf{x}}|{\mathbf{z}}_0) \right]
    }_{\text{Reconstruction loss}}
    +
    \underbrace{
    \mathcal{L}_T({\mathbf{x}}).
    }_{\text{Diffusion loss}}
    \label{eq:finite_elbo_1}
\end{align}
The prior loss and reconstruction loss can be (stochastically and differentiably) estimated using standard techniques; see \citep{kingma2013auto}. The diffusion loss, $\mathcal{L}_T({\mathbf{x}})$, is more complicated, and depends on the hyperparameter $T$ that determines the depth of the generative model.

\section{Discrete-time model}

In the case of finite $T$, using $s(i) = (i-1)/T$, $t(i) = i/T$, the diffusion loss is:
\begin{align}
    \mathcal{L}_T({\mathbf{x}})
    &=
    \sum_{i=1}^T
    \mathbb{E}_{q({\mathbf{z}}_{t(i)}|{\mathbf{x}})}
    D_{KL}[q({\mathbf{z}}_{s(i)}|{\mathbf{z}}_{t(i)},{\mathbf{x}})||p({\mathbf{z}}_{s(i)}|{\mathbf{z}}_{t(i)})].
    \label{eq:main_finite_ELBO}
\end{align}
In appendix \ref{sec:objective} we show that this expression simplifies considerably, yielding:
\begin{align}
\mathcal{L}_T({\mathbf{x}}) =
\frac{T}{2}
\mathbb{E}_{{\boldsymbol{\epsilon}} \sim \mathcal{N}(0,\mathbf{I}), i \sim U\{1, T\}}
\left[
\left(\text{SNR}(s)-\text{SNR}(t)\right)
||{\mathbf{x}} - \hat{{\mathbf{x}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}_t;t) ||_2^2
\right],
\label{eq:lT_main}\end{align}
where $U\{1, T\}$ is the uniform distribution on the integers $\{1,\ldots,T\}$, and ${\mathbf{z}}_t = \alpha_t {\mathbf{x}} + \sigma_t {\boldsymbol{\epsilon}}$. This is the general discrete-time loss for any choice of forward diffusion parameters $(\sigma_t,  \alpha_t)$.  When plugging in the specifications of $\sigma_t$, $\alpha_t$ and $\hat{{\mathbf{x}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}_t; t)$ that we use in experiments, given in Sections \ref{sec:learning_noise} and \ref{sec:noisepred}, the loss simplifies to:
\begin{align}
\mathcal{L}_T({\mathbf{x}}) =
\frac{T}{2}
\mathbb{E}_{{\boldsymbol{\epsilon}} \sim \mathcal{N}(0,\mathbf{I}), i \sim U\{1, T\}}
\left[
(\exp(\gamma_{{\boldsymbol{\eta}}}(t)-\gamma_{{\boldsymbol{\eta}}}(s))-1)
\left\lVert{\boldsymbol{\epsilon}} - \hat{{\boldsymbol{\epsilon}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}_t; t) \right\lVert_2^2
\right]
\label{eq:discrete}\end{align}
where ${\mathbf{z}}_t = \text{sigmoid}(-\gamma_{{\boldsymbol{\eta}}}(t)) {\mathbf{x}} + \text{sigmoid}(\gamma_{{\boldsymbol{\eta}}}(t)) {\boldsymbol{\epsilon}}$.
In the discrete-time case, we simply jointly optimize ${\boldsymbol{\eta}}$ and ${\boldsymbol{\theta}}$ by maximizing the VLB through a Monte Carlo estimator of \Eqref{eq:discrete}.

Note that $\exp(.)-1$ has a numerically stable primitive $expm1(.)$ in common numerical computing packages; see figure \ref{fig:expm1}. \Eqref{eq:discrete} allows for numerically stable implementation in 32-bit or lower-precision floating point, in contrast with previous implementations of discrete-time diffusion models (e.g. ~\citep{ho2020denoising}), which had to resort to 64-bit floating point.

\subsection{More steps leads to a lower loss}
\label{sec:moresteps}
A natural question to ask is what the number of timesteps $T$ should be, and whether more timesteps is always better in terms of the VLB. In Appendix~\ref{appendix:more_steps} we analyze the difference between the diffusion loss with $T$ timesteps, $\mathcal{L}_T({\mathbf{x}})$, and the diffusion loss with double the timesteps, $\mathcal{L}_{2T}({\mathbf{x}})$, while keeping the SNR function fixed.
We then find that if our trained denoising model $\hat{{\mathbf{x}}}_{{\boldsymbol{\theta}}}$ is sufficiently good, we have that $\mathcal{L}_{2T}({\mathbf{x}}) < \mathcal{L}_T({\mathbf{x}})$, i.e.\ that our VLB will be better for a larger number of timesteps. Intuitively, the discrete time diffusion loss is an upper Riemann sum  approximation of an integral of a strictly decreasing function, meaning that a finer approximation yields a lower diffusion loss. This result is illustrated in Figure~\ref{fig:more_steps_is_better}.
\begin{figure}[t]
\includegraphics[width=0.9\textwidth]{media/images/vdm-riemann-figure.pdf}
  \caption{
Illustration of the diffusion loss with few segments $T$ (left), more segments $T$ (middle), and infinite segments $T$ (continuous time, right). The continuous-time loss (\Eqref{eq:elbo_v}) is an integral of mean squared error (MSE) over SNR, here visualized as a black curve. The black curve is strictly decreasing when the model is sufficiently well trained, so the discrete-time loss (\Eqref{eq:lT_main}) is an upper bound (an upper Riemann sum approximation) of this integral that becomes better when segments are added.}
  \label{fig:more_steps_is_better}
\end{figure}

\section{Continuous-time model: $T\rightarrow\infty$}

Since taking more time steps leads to a better VLB, we now take $T\rightarrow\infty$, effectively treating time $t$ as continuous rather than discrete. The model for $p({\mathbf{z}}_t)$ can in this case be described as a continuous time diffusion process \citep{song2020score} governed by a stochastic differential equation; see Appendix \ref{sec:sde}. In Appendix \ref{sec:objective} we show that in this limit the diffusion loss $\mathcal{L}_T({\mathbf{x}})$ simplifies further. Letting $\text{SNR}'(t) = d \text{SNR}(t)/dt$, we have, with ${\mathbf{z}}_t = \alpha_t {\mathbf{x}} + \sigma_t {\boldsymbol{\epsilon}}$:
\begin{align}
\mathcal{L}_{\infty}({\mathbf{x}})
&= -\frac{1}{2}\mathbb{E}_{{\boldsymbol{\epsilon}}\sim\mathcal{N}(0,\mathbf{I})} \int_{0}^{1} \text{SNR}'(t) \left\rVert {\mathbf{x}} - \hat{{\mathbf{x}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}_t;t) \right\lVert_{2}^{2} dt, \label{eq:fusion_integral}
\\
&= -\frac{1}{2}\mathbb{E}_{{\boldsymbol{\epsilon}}\sim\mathcal{N}(0,\mathbf{I}), t \sim \mathcal{U}(0,1)}
\left[ \text{SNR}'(t) \left\rVert {\mathbf{x}} - \hat{{\mathbf{x}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}_t;t) \right\lVert_{2}^{2} \right].
\label{eq:fusion}\end{align}
This is the general continuous-time loss for any choice of forward diffusion parameters $(\sigma_t,  \alpha_t)$. When plugging in the specifications of $\sigma_t$, $\alpha_t$ and $\hat{{\mathbf{x}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}_t; t)$ that we use in experiments, given in Sections \ref{sec:learning_noise} and \ref{sec:noisepred}, the loss simplifies to:
\begin{align}
\mathcal{L}_{\infty}({\mathbf{x}})
&= \frac{1}{2}\mathbb{E}_{{\boldsymbol{\epsilon}}\sim\mathcal{N}(0,\mathbf{I}), t \sim \mathcal{U}(0,1)}
\left[ \gamma'_{{\boldsymbol{\eta}}}(t) \left\rVert {\boldsymbol{\epsilon}} - \hat{{\boldsymbol{\epsilon}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}_t;t) \right\lVert_{2}^{2} \right],
\end{align}
where $\gamma'_{{\boldsymbol{\eta}}}(t) = d \gamma_{{\boldsymbol{\eta}}}(t)/dt$. We use the Monte Carlo estimator of this loss for evaluation and optimization.

\subsection{Equivalence of diffusion models in continuous time}
\label{sec:schedule_invariance}
The signal-to-noise function $\text{SNR}(t)$ is invertible due to the monotonicity assumption in Section \ref{sec:inferencemodel}. Due to this invertibility, we can perform a change of variables, and make everything a function of $v \equiv \text{SNR}(t)$ instead of $t$, such that $t = \text{SNR}^{-1}(v)$. Let $\alpha_v$ and $\sigma_v$ be the functions $\alpha_t$ and $\sigma_t$ evaluated at $t = \text{SNR}^{-1}(v)$, and correspondingly let ${\mathbf{z}}_v = \alpha_v {\mathbf{x}} + \sigma_v {\boldsymbol{\epsilon}}$. Similarly, we rewrite our noise prediction model as $\tilde{{\mathbf{x}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}, v) \equiv \hat{{\mathbf{x}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}, \text{SNR}^{-1}(v))$. With this change of variables, our continuous-time loss in \Eqref{eq:fusion_integral} can equivalently be written as:
\begin{align}
\mathcal{L}_{\infty}({\mathbf{x}}) = \frac{1}{2}\mathbb{E}_{{\boldsymbol{\epsilon}}\sim\mathcal{N}(0,\mathbf{I})} \int_{\text{SNR}_{\text{min}}}^{\text{SNR}_{\text{max}}} \left\rVert {\mathbf{x}} - \tilde{{\mathbf{x}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}_v, v) \right\lVert_{2}^{2} dv,
\label{eq:elbo_v}
\end{align}
where instead of integrating w.r.t. time $t$ we now integrate w.r.t. the signal-to-noise ratio $v$, and where $\text{SNR}_\text{min}=\text{SNR}(1)$ and $\text{SNR}_\text{max}=\text{SNR}(0)$.

What this equation shows us is that the only effect the functions $\alpha(t)$ and $\sigma(t)$ have on the diffusion loss is through the values $\text{SNR}(t) = \alpha^{2}_t/\sigma^{2}_t$ at endpoints $t=0$ and $t=1$. Given these values $\text{SNR}_\text{max}$ and $\text{SNR}_\text{min}$, the diffusion loss is invariant to the shape of function $\text{SNR}(t)$ between $t=0$ and $t=1$. The VLB is thus only impacted by the function $\text{SNR}(t)$ through its endpoints $\text{SNR}_{\text{min}}$ and $\text{SNR}_{\text{max}}$.

Furthermore, we find that the distribution $p({\mathbf{x}})$ defined by our generative model is also invariant to the specification of the diffusion process. Specifically, let $p^{A}({\mathbf{x}})$ denote the distribution defined by the combination of a diffusion specification and denoising function $\{\alpha^A_v, \sigma^A_v, \tilde{{\mathbf{x}}}^A_{{\boldsymbol{\theta}}}\}$, and similarly let $p^{B}({\mathbf{x}})$ be the distribution defined through a different specification $\{\alpha^B_v, \sigma^B_v, \tilde{{\mathbf{x}}}^B_{{\boldsymbol{\theta}}}\}$, where both specifications have equal $\text{SNR}_\text{min}, \text{SNR}_\text{max}$; as shown in Appendix~\ref{sec:app_equivalence}, we then have that $p^{A}({\mathbf{x}}) = p^{B}({\mathbf{x}})$ if $\tilde{{\mathbf{x}}}^B_{{\boldsymbol{\theta}}}({\mathbf{z}}, v) \equiv \tilde{{\mathbf{x}}}^A_{{\boldsymbol{\theta}}}((\alpha^A_v/\alpha^B_v){\mathbf{z}}, v)$. The distribution on all latents ${\mathbf{z}}_v$ is then also the same under both specifications, up to a trivial rescaling. This means that any two diffusion models satisfying the mild constraints set in 3.1 (which includes e.g.\ the \emph{variance exploding} and \emph{variance preserving} specifications considered by~\cite{song2020score}), can thus be seen as equivalent in continuous time.

\subsection{Weighted diffusion loss}
This equivalence between diffusion specifications continues to hold even if, instead of the VLB, these models optimize a \emph{weighted} diffusion loss of the form:
\begin{align}
\mathcal{L}_{\infty}({\mathbf{x}}, w) = \frac{1}{2}\mathbb{E}_{{\boldsymbol{\epsilon}}\sim\mathcal{N}(0,\mathbf{I})} \int_{\text{SNR}_{\text{min}}}^{\text{SNR}_{\text{max}}} w(v) \left\rVert {\mathbf{x}} - \tilde{{\mathbf{x}}}_{{\boldsymbol{\theta}}}({\mathbf{z}}_v, v) \right\lVert_{2}^{2} dv,
\label{eq:weighted_diff}
\end{align}
which e.g.\ captures all the different objectives discussed by \cite{song2020score}, see Appendix~\ref{sec:ncsn_ddpm_weights}. Here, $w(v)$ is a weighting function that generally puts increased emphasis on the noisier data compared to the VLB, and which thereby can sometimes improve perceptual generation quality as measured by certain metrics like FID and Inception Score. For the models presented in this paper, we further use $w(v)=1$ as corresponding to the (unweighted) VLB.

\subsection{Variance minimization}

Lowering the variance of the Monte Carlo estimator of the continuous-time loss generally improves the efficiency of optimization. We found that using a low-discrepancy sampler for $t$, as explained in Appendix \ref{sec:lowdisc}, leads to a significant reduction in variance. In addition, due to the invariance shown in Section~\ref{sec:schedule_invariance} for the continous-time case, we can optimize the schedule \emph{between} its endpoints w.r.t. to minimize the variance of our estimator of loss, as detailed in Appendix \ref{app:variance_minimization}. The endpoints of the noise schedule are simply optimized w.r.t. the VLB.

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

\begin{table}[t]
\footnotesize
\begin{center}
\setlength{\tabcolsep}{4pt}
\begin{tabular}{lccccc}
\textbf{Model} & Type & CIFAR10 & CIFAR10 & ImageNet & ImageNet %
\\
(Bits per dim on test set)&  & no data aug. & data aug. & 32x32 & 64x64%
\\
\midrule
\textit{Previous work}\\
ResNet VAE with IAF~\citep{kingma2016improving} & VAE & 3.11 & \\
Very Deep VAE~\citep{child2020very} & VAE & 2.87 & & 3.80 & 3.52 \\
NVAE~\citep{vahdat2020nvae} & VAE & 2.91 & & 3.92 & \\
Glow~\citep{kingma2018glow} & Flow & & $3.35^{(B)}$ & 4.09 & 3.81 \\
Flow++~\citep{ho2019flow++} & Flow & 3.08 & & 3.86 & 3.69 \\
PixelCNN~\citep{van2016pixel} & AR & 3.03 & & 3.83 & 3.57\\
PixelCNN++~\citep{salimans2017pixelcnn++} & AR & 2.92 &\\
Image Transformer~\citep{parmar2018image} & AR & 2.90 & & 3.77\\
SPN ~\citep{menick2018generating} & AR & & & & 3.52 \\
Sparse Transformer~\citep{child2019generating}& AR &2.80 & & & 3.44 %
\\
Routing Transformer~\citep{roy2021efficient} & AR & & & & 3.43\\
Sparse Transformer + DistAug~\citep{jun2020distribution} & AR & & $2.53^{(A)}$ & & \\
DDPM~\citep{ho2020denoising} & Diff & & $3.69^{(C)}$ & & %
\\
EBM-DRL~\citep{gao2020learning} & Diff & & $3.18^{(C)}$ & & \\
Score SDE~\citep{song2020score} & Diff & 2.99 & & & %
\\
Improved DDPM~\citep{nichol2021improved} & Diff & 2.94 & & & 3.54 %
\\
\midrule
\textit{Concurrent work}\\
CR-NVAE \citep{sinha2021consistency} & VAE & & $2.51^{(A)}$ \\
LSGM~\citep{vahdat2021score} & Diff & 2.87\\
ScoreFlow~\citep{song2021maximum} (variational bound) & Diff & & $2.90^{(C)}$ & 3.86 & \\
ScoreFlow~\citep{song2021maximum} (cont. norm. flow) & Diff & 2.83 & $2.80^{(C)}$ & 3.76 & \\
\midrule
\textit{Our work}\\
\textbf{VDM (variational bound)} & Diff & \textbf{2.65} & $\textbf{2.49}^{(A)}$ & \textbf{3.72} & \textbf{3.40} %
\\
\bottomrule
\end{tabular}
\end{center}
\caption{Summary of our findings for density modeling tasks, in terms of bits per dimension (BPD) on the test set. Model types are autoregressive (AR), normalizing flows (Flow), variational autoencoders (VAE), or diffusion models (Diff). Our results were obtained using the continuous-time formulation of our model. CIFAR-10 data augmentations are: (A) extensive, (B) small translations, or (C) horizontal flips. The numbers for VDM are variational bounds, and can likely be improved by estimating the marginal likelihood through importance sampling, or through evaluation of the corresponding continuous normalizing flow as done by \cite{song2021maximum}.}
\label{table:results}
\end{table}

We demonstrate our proposed class of diffusion models, which we call \emph{Variational Diffusion Models} (VDMs), on the CIFAR-10~\citep{krizhevsky2009learning} dataset, and the downsampled ImageNet~\citep{van2016pixel, deng2009imagenet} dataset, where we focus on maximizing likelihood. For our result with data augmentation we used random flips, 90-degree rotations, and color channel swapping. More details of our model specifications are in Appendix~\ref{sec:hps}.

\subsection{Likelihood and samples}
Table~\ref{table:results} shows our results on modeling the CIFAR-10 dataset, and the downsampled ImageNet dataset. We establish a new state-of-the-art in terms of test set likelihood on all considered benchmarks, by a significant margin. Our model for CIFAR-10 without data augmentation surpasses the previous best result of $2.80$ about 10x faster than it takes the Sparse Transformer to reach this, in wall clock time on equivalent hardware.
Our CIFAR-10 model, whose hyper-parameters were tuned for likelihood, results in a FID (perceptual quality) score of 7.41. This would have been state-of-the-art until recently, but is worse than recent diffusion models that specifically target FID scores ~\citep{nichol2021improved,song2020score,ho2020denoising}. By instead using a weighted diffusion loss, with the weighting function $w(\text{SNR})$ used by \cite{ho2020denoising} and described in Appendix~\ref{sec:ncsn_ddpm_weights}, our FID score improves to 4.0. We did not pursue further tuning of the model to improve FID instead of likelihood.
A random sample of generated images from our model is provided in Figure~\ref{fig:imagenet64}. We provide additional samples from this model, as well as our other models for the other datasets, in Appendix~\ref{sec:app_more_samples}.
\begin{figure}[t]
    \centering
    \includegraphics[width=\textwidth]{media/images/imn64.png}
    \caption{Non cherry-picked unconditional samples from our Imagenet 64x64 model, trained in continuous time and generated using $T=1000$. The model's hyper-parameters and parameters are optimized w.r.t. the likelihood bound, so the model is not optimized for synthesis quality.}
    \label{fig:imagenet64}
\end{figure}

\subsection{Ablations}
Next, we investigate the relative importance of our contributions. In Table~\ref{table:t_results} we compare our discrete-time and continuous-time specifications of the diffusion model: When evaluating our model with a small number of steps, our discretely trained models perform better by learning the diffusion schedule to optimize the VLB. However, as argued theoretically in Section~\ref{sec:moresteps}, we find experimentally that more steps $T$ indeed gives better likelihood. When $T$ grows large, our continuously trained model performs best, helped by training its diffusion schedule to minimize variance instead. %

Minimizing the variance also helps the continuous time model to train faster, as shown in Figure~\ref{fig:fourier-ablation}. This effect is further examined in Table~\ref{table:variance}, where we find dramatic variance reductions compared to our baselines in continuous time. Figure~\ref{fig:snr_schedule} shows how this effect is achieved: Compared to the other schedules, our learned schedule spends much more time in the high $\text{SNR}(t)$ / low $\sigma^{2}_{t}$ range.

In Figure~\ref{fig:fourier-ablation} we further show training curves for our model including and excluding the Fourier features proposed in Appendix~\ref{sec:fourier}: with Fourier features enabled our model achieves much better likelihood. For comparison we also implemented Fourier features in a PixelCNN++ model~\citep{salimans2017pixelcnn++}, where we do not see a benefit. In addition, we find that learning the SNR is necessary to get the most out of including Fourier features: if we fix the SNR schedule to that used by \cite{ho2020denoising}, the maximum log-SNR is fixed to approximately 8 (see figure~\ref{fig:weighting_funs}), and test set negative likelihood stays above 4 bits per dim. When learning the SNR endpoints, our maximum log-SNR ends up at $13.3$, which, combined with the inclusion of Fourier features, leads to the SOTA test set likelihoods reported in Table~\ref{table:results}.

\subsection{Lossless compression}
For a fixed number of evaluation timesteps $T_{eval}$, our diffusion model in discrete time is a hierarchical latent variable model that can be turned into a lossless compression algorithm using bits-back coding~\citep{hinton1993keeping}.
As a proof of concept of practical lossless compression, Table~\ref{table:t_results} reports net codelengths on the CIFAR10 test set for various settings of $T_{eval}$ using BB-ANS~\citep{townsend2018practical}, an implementation of bits-back coding based on asymmetric numeral systems~\citep{duda2009asymmetric}. Details of our implementation are given in Appendix~\ref{sec:compression}. We achieve state-of-the-art net codelengths, proving our model can be used as the basis of a lossless compression algorithm. However, for large $T_{eval}$ a gap remains with the theoretically optimal codelength corresponding to the negative VLB, and compression becomes computationally expensive due to the large number of neural network forward passes required. Closing this gap with more efficient implementations of bits-back coding suitable for very deep models is an interesting avenue for future work.

\begin{figure}[t]
    \centering\footnotesize
    \begin{subfigure}[b]{0.35\textwidth}
        \centering
        \includegraphics[width=\textwidth]{media/images/cifar10/snr.png}
        \caption{$\log\text{SNR}$ vs time $t$ }
        \label{fig:snr_schedule}
    \end{subfigure}\hspace{1.5cm}
    \begin{subfigure}[b]{0.35\textwidth}
        \centering
        \begin{tabular}{rr}
        $\text{SNR}(t)$~schedule & $\mathrm{Var}(\text{BPD})$\\
        \midrule
        \textbf{Learned (ours)} & \textbf{0.53}\\
        $\log\text{SNR}$-linear & 6.35\\
        $\beta$-Linear [1] & 31.6\\
        $\alpha$-Cosine [2] & 31.1\\
        \bottomrule
        \end{tabular}
        \caption{Variance of VLB estimate}
        \label{table:variance}
    \end{subfigure}
    \caption{Our learned continuous-time variance-minimizing noise schedule $\text{SNR}(t)$ for CIFAR-10, compared to its log-linear initialization and to schedules from the literature: [1] The $\beta$-Linear schedule from \cite{ho2020denoising}, [2] The $\alpha$-Cosine schedule from \cite{nichol2021improved}. All schedules were scaled and shifted on the log scale such that the resulting $\text{SNR}_{\text{min}}, \text{SNR}_{\text{max}}$ were the equal to our learned endpoints, resulting in the same VLB estimate of 2.66. We report the variance of our VLB estimate per data point, computed on the test set, and conditional on the data: This does not include the noise due to sampling minibatches of data.}
\end{figure}

\begin{figure}[t]
\begin{floatrow}
\capbtabbox{%
\scriptsize
\begin{tabular}{rrrc}
$T_{train}$ & $T_{eval}$ & BPD & \makecell{Bits-Back \\ Net BPD} \\
\midrule
10 & 10 & 4.31 \\
100 & 100 & 2.84 \\
250 & 250 & 2.73 \\
500 & 500 & 2.68 \\
1000 & 1000 & 2.67 \\
10000 & 10000 & 2.66 \\
\midrule
$\infty$ & 10 & 7.54 & 7.54 \\
$\infty$ & 100 & 2.90 & 2.91 \\
$\infty$ & 250 & 2.74 & 2.76 \\
$\infty$ & 500 & 2.69 & 2.72 \\
$\infty$ & 1000 & 2.67 & 2.72 \\
$\infty$ & 10000 & \textbf{2.65} \\
\midrule
$\infty$ & $\infty$ & \textbf{2.65}\\
\bottomrule
\end{tabular}
}{%
\caption{Discrete versus continuous-time training and evaluation with CIFAR-10, in terms of bits per dimension (BPD).}%
\label{table:t_results}
}
\ffigbox{%
  \includegraphics[width=0.5\textwidth]{media/images/cifar10/fourier_ablation.pdf}%
}{%
  \caption{Test set likelihoods during training, with/without Fourier features, and with/without learning the noise schedule to minimize variance.}%
  \label{fig:fourier-ablation}
}
\end{floatrow}
\end{figure}

\section{Conclusion}
\label{sec:conclusion}
We presented state-of-the-art results on modeling the density of natural images using a new class of diffusion models that incorporates a learnable diffusion specification, Fourier features for fine-scale modeling, as well as other architectural innovations. In addition, we obtained new theoretical insight into likelihood-based generative modeling with diffusion models, showing a surprising invariance of the VLB to the forward time diffusion process in continuous time, as well as an equivalence between various diffusion processes from the literature previously thought to be different.

\section*{Acknowledgments}
We thank Yang Song, Kevin Murphy, Mohammad Norouzi and Chin-Yun Yu for helpful feedback on the paper, and Ruiqi Gao for helping with writing an open source version of the code.

\bibliographystyle{plainnat}