text:
\section{Standard Maximum Likelihood Estimation and Links to \textsc{I-MLE}\@\xspace}

\label{apx-mle}

In the standard MLE setting \citep[see, e.g., ][Ch. 9]{murphy2012machine} we are interested in learning the parameters of a probability distribution, here assumed to be from the (constrained) exponential family (see \cref{def-constrained-exp-family}), given a set of example states.
More specifically,
given training data  $\mathcal{D} = \{\hat{\mathbf{z}}_j\}_{j=1}^N$, with $\hat{\mathbf{z}}_j\in \mathcal{C}\subseteq \{0, 1\}^m$, maximum-likelihood learning aims to minimize the empirical risk
\begin{equation}
\label{eq-apx-std-MLE-loss}
\mathcal{L}(\bm{\theta}, q_\mathcal{D})
= \mathbb{E}_{\hat{\bz}\sim q_{\mathcal{D}}}[- \log p(\hat{\mathbf{z}}; \bm{\theta}) ]
=
\frac{1}{N} \sum_{j=1}^{N} -\log p(\hat{\mathbf{z}}_j; \bm{\theta}) =
\frac{1}{N} \sum_{j=1}^{N} \left( A(\bm{\theta})
- \langle\hat{\mathbf{z}}_j,\bm{\theta}\rangle
\right)
\end{equation}
with respect to $\bm{\theta}$, where $q_{\mathcal{D}}(\bm{z})=\sum_j \delta_{\hat{\bz}_j}(\bm{z}) /N $
is the empirical data distribution and $\delta_{\hat{\bz}}$ is the Dirac delta centered in $\hat{\bz}$.
In \cref{eq-apx-std-MLE-loss}, the point-wise loss $\ell$ is the negative log likelihood $- p(\hat{\bz}, \bm{\theta})$, and $q_{\mathcal{D}}$ may be seen as a (data/empirical) \emph{target distribution}.
Note that in the main paper, as we assumed $q$ to be from the exponential family with parameters $\bm{\theta}'$, we used the notation $\mathcal{L}(\bm{\theta}, \bm{\theta}')$ to indicate the MLE objective rather than $\mathcal{L}(\bm{\theta}, q)$.
These two definitions are, however,  essentially equivalent.

\cref{eq-apx-std-MLE-loss} is a smooth objective that can be optimized with a (stochastic) gradient descent procedure.
For a data point $\hat{\mathbf{z}}$, the gradient of the point-wise loss is given by $\grad{\bm{\theta}}{\ell} = \bm{\mu}(\bm{\theta}) - \hat{\mathbf{z}}$, since $\grad{\bm{\theta}}{A}=\bm{\mu}$.
For the entire dataset one has
\begin{equation}
    \label{eq-apx-grad_smle}
    \grad{\bm{\theta}}{\mathcal{L}(\bm{\theta}, q_{\mathcal{D}})} = \mu(\bm{\theta}) - \frac{1}{N} \sum_{j=1}^N \hat{\bz}_j =
    \mu(\bm{\theta}) - \mathbb{E}_{\hat{\bz}\sim q_{\mathcal{D}}}[ \hat{\bz}_j ]
\end{equation}
which is (cf.  \cref{eq-grad-mleobj}) the difference between the marginals of $p(\bz; \btheta)$ (the mean of $\mathbf{Z}$) and the empirical mean of $\mathcal{D}$,
$
\bm{\mu}(q_{\mathcal{D}}) =
\mathbb{E}_{\hat{\bz}\sim\mathcal{D}}[\hat{\bz}]
$.
As mentioned in the main paper, the main computational challenge when evaluating \cref{eq-apx-grad_smle} is to compute the marginals (of $p(\bz; \btheta)$).
There are many approximate schemes, one of which is the so-called \emph{perceptron rule}, which approximate \cref{eq-apx-grad_smle} as
\begin{equation*}
    \widehat{\nabla}_{\btheta}\mathcal{L}(\bm{\theta}, q_{\mathcal{D}} ) = \fMAP(\bm{\theta}) - \frac{1}{N} \sum_{j=1}^N \hat{\bz}_j
\end{equation*}
and it is frequently employed in a stochastic manner by sampling one or more points form $\mathcal{D}$, rather than computing the full dataset mean.

We may interpret the standard MLE setting described in this section from the perspective of the problem setting we presented in \cref{sec-ps}.
The first indeed
amounts to the special case of the latter where there are no inputs ($\mathcal{X}=\emptyset$), the  (target) output space coincides with the state space of the distribution ($\mathcal{Y}=\mathcal{Z}$), $f$ is the identity mapping, $\ell$ is the negative log-likelihood and the model's parameter coincide with the distribution parameters, that is $\bm{\omega}=\bm{\theta}$.

\section{Algorithm} \label{app:algorithm}

\cref{algo:main} shows the pseudo-code of the algorithm implementing \cref{eqn-perturb-map-approx} for $S=1$, using perturbation-based implicit differentiation target distributions (see \cref{subsection-implicit-differentiation}).
This demonstrates that \textsc{I-MLE}\@\xspace can be implemented and used as a layer in deep learning frameworks such as PyTorch and TensorFlow.

\input{passes.tex}

\section{Proofs of Section \ref{section-perturb-and-map} and Section~\ref{section-optimal-q-co}}

\label{proofs-imle}

This section contains the proofs of the results relative to the perturb and map section (\cref{section-perturb-and-map}) and the section on optimal target distributions for typical loss functions when backpropagating through combinatorial optimization problems (section~\ref{section-optimal-q-co}).
We repeat the statements here, for the convenience of the reader.

\begin{mirrorproposition}
Let $p(\bm{z}; \bm{\theta})$ be a discrete exponential family distribution with constraints $\mathcal{C}$ and temperature $\tau$, and let $\langle\bm{\theta},\bm{z}\rangle$ be the unnormalized weight of each $\bm{z}$ with $\bm{z} \in \mathcal{C}$. Moreover, let $\tilde{\bm{\theta}}$ be such that, for all $\bm{z} \in \mathcal{C}$,
$$\langle\bm{z},\tilde{\bm{\theta}}\rangle =  \langle\bm{z},\bm{\theta}\rangle + \epsilon(\bm{z})$$ with each
$\epsilon(\bm{z})$ i.i.d. samples from $\GumbelDist(0,\tau)$. Then, $$\Pr\left(\mathtt{MAP}(\tilde{\bm{\theta}})=\bm{z}\right) = p(\bm{z}; \bm{\theta}).$$
\end{mirrorproposition}

\begin{proof}
Let $\epsilon_i \sim \GumbelDist(0,\tau)$ i.i.d. and $\tilde{\theta}_i = \theta_i + \epsilon_i$. Following a derivation similar to one made in \citet{Papandreou:2011}, we have:
\begin{align*}
    & \Pr\{\argmax(\tilde{\theta}_1, \ldots, \tilde{\theta}_m) = n\}  = \\
    = & \Pr\{ \tilde{\theta}_n \geq \max_{j\neq n} \{\tilde{\theta}_j\} \}  \\
    = & \int_{-\infty}^{+\infty} g(t; \theta_n) \prod_{j\neq n} G(t; \theta_j) \,dt \\
    = & \int_{-\infty}^{+\infty} \frac{1}{\tau} \exp\left(\frac{\theta_n - t}{\tau} - e^{\frac{\theta_n - t}{\tau}}\right) \prod_{j\neq n} \exp\left( -e^{\frac{\theta_j - t}{\tau}}\right) \,dt \\
    = &  \int_{-\infty}^{+\infty} \frac{1}{\tau} e^{\frac{\theta_n - t}{\tau}}\exp\left(-e^{\frac{\theta_n-t}{\tau}}\right) \prod_{j\neq n} \exp\left(-e^{\frac{\theta_j - t}{\tau}}\right) \,dt \\
    = & \int_{0}^{1} \prod_{j \neq n} z^{\exp\left( \frac{\theta_j - \theta_n}{\tau}  \right)} \,dz \qquad \text{with } z \triangleq \exp\left( -e^{\frac{\theta_n - t}{\tau}} \right) \\
    = & \frac{1}{1 + \sum_{j \neq n} e^{\frac{\theta_j - \theta_n}{\tau}} } \\
    = & \frac{e^{\frac{\theta_n}{\tau}}}{\sum_{j=1}^{m} e^{\frac{\theta_j}{\tau}}, }
  \end{align*}
  where $g$ and $G$ are respectively the Gumbel probability density function and the Gumbel cumulative density function.
  The proposition now follows from arguments made in \citet{Papandreou:2011} using the maximal equivalent re-parameterization of $p(\bm{z}; \bm{\theta})$ where we specify a parameter $\langle \bm{\theta},\bm{z}\rangle$ for each $\bm{z}$ with $\mathcal{C}(\bm{z})$ and perturb these parameters.
\end{proof}

\begin{mirrorlemma}
Let $X \sim \mbox{Gumbel}(0 ,\tau)$ and let $\kappa \in \mathbb{N} \setminus \{0\}$.  Then we can write
\begin{equation*}
X \sim \sum_{j=1}^{\kappa} \frac{\tau}{\kappa} \left[ \lim_{s \rightarrow \infty}\sum_{i=1}^{s} \left\{\GammaDist(1/\kappa, \kappa/i)\right\} - \log(s)\right],
\end{equation*}
where $\GammaDist(\alpha, \beta)$ is the gamma distribution with shape $\alpha$ and scale $\beta$.
\end{mirrorlemma}

\begin{proof}
Let $\kappa \in \mathbb{N} \setminus \{0\}$ and let $X \sim \GumbelDist(0, \tau)$. Its moment generating function has the form
\begin{equation}
\label{eq:m-gen-f-gumb}
\mathbb{E}[\exp(tX)] = \Gamma(1 - \tau t).
\end{equation}
As mentioned in \citet[p. 443,][]{johnson1998advances} we know that we can write the Gamma function as
\begin{equation}
\label{eq-gamma-book}
\Gamma(1 - \tau t) = e^{\gamma \tau t} \prod_{i=1}^{\infty} \left(1 - \frac{\tau t}{i}\right)^{-1} e^{\frac{-\tau t}{i}}
\end{equation}
where $\gamma$ is the Euler-Mascheroni constant.
We have that
$$\left(1 - \frac{\tau t}{i}\right)^{-1} = \frac{i}{i - \tau t} = \frac{\frac{i}{\tau}}{\frac{i - \tau t}{\tau}} =  \frac{\frac{i}{\tau}}{\frac{i}{\tau} - \frac{\tau t}{\tau}} =
\frac{\frac{i}{\tau}}{\frac{i}{\tau} - t}.$$
The last term is the moment generating function of an exponential distribution with scale $\frac{\tau}{i}$.
We can now take the logarithm on both sides of \cref{eq-gamma-book} and obtain
\begin{equation*}
  tX = \gamma\tau t + \lim_{s\rightarrow \infty}\sum_{i=1}^{s} \left( t \ExpDist(\tau/i) - \frac{\tau t}{i} \right),
\end{equation*}
where $\ExpDist(\alpha)$ is the exponential distribution with scale $\alpha$.
Hence,
\begin{align*}
 X \sim & \lim_{s\rightarrow \infty}\sum_{i=1}^{s}  \left( \ExpDist(\tau/i) - \frac{\tau}{i} \right) +  \gamma\tau \\
  = & \lim_{s\rightarrow \infty}\sum_{i=1}^{s}  \left( \ExpDist(\tau/i) - \frac{\tau}{i} \right) + \tau \lim_{s\rightarrow \infty}\sum_{i=1}^{s} \frac{1}{i} - \log(s) \\
  = & \lim_{s\rightarrow \infty}\sum_{i=1}^{s}  \left( \ExpDist(\tau/i) - \frac{\tau}{i} + \frac{\tau}{i} \right) - \tau\log(s)   \\
    = & \lim_{s\rightarrow \infty}\sum_{i=1}^{s}  \ExpDist(\tau/i) - \tau\log(s)
\end{align*}

Since $\ExpDist(\alpha) \sim \GammaDist(1, \alpha)$, and due to the scaling and summation properties of the Gamma distribution (with shape-scale parameterization), we can write for all $r>1$:
$$\ExpDist(\alpha) \sim \sum_{j=1}^{r} \GammaDist(1/r, \alpha r)/r.$$
Hence, picking $r=\kappa$ from the hypothesis, we have

\begin{align*}
    X \sim & \lim_{s\rightarrow \infty} \left\lbrace
        \sum_{i=1}^{s} \sum_{j=1}^{\kappa} \GammaDist(1/\kappa, \tau \kappa/i)/\kappa
    \right\rbrace - \tau \log(s) \\
    = &
    \lim_{s\rightarrow \infty} \left\lbrace
        \sum_{j=1}^{\kappa} \frac{\tau}{k}\sum_{i=1}^{s} \GammaDist(1/\kappa, \kappa/i)
    \right\rbrace - \sum_{j=1}^\kappa \frac{\tau}{\kappa} \log(s) \\
    = &
    \lim_{s\rightarrow \infty}
    \sum_{j=1}^{\kappa} \frac{\tau}{\kappa}
    \left\lbrace
        \left[
            \sum_{i=1}^{s} \GammaDist(1/\kappa, \kappa/i)
        \right] - \log(s)
    \right\rbrace \\
    = &
    \sum_{j=1}^{\kappa}

    \frac{\tau}{\kappa}
    \left\lbrace
        \lim_{s\rightarrow \infty}
        \left[
            \sum_{i=1}^{s} \GammaDist(1/\kappa, \kappa/i)
        \right] - \log(s)
    \right\rbrace

\end{align*}
This concludes the proof.
Parts of the proof are inspired by a post on stackexchange~\cite{214875}.
\end{proof}

\begin{mirrortheorem}
Let $p(\bm{z}; \bm{\theta})$ be a discrete exponential family distribution with constraints $\mathcal{C}$ and temperature $\tau$, and let $k \in \mathbb{N} \setminus \{0\}$. Let us assume that if $\mathcal{C}(\bm{z})$ then $\langle\bm{z}, \mathbf{1}\rangle = k$.
Let $\tilde{\bm{\theta}}$ be the perturbation obtained by $\tilde{\bm{\theta}}_i = \bm{\theta}_i + \epsilon_i$
with
\begin{equation}
    \epsilon_i \sim \frac{\tau}{k} \left[\lim_{s\to\infty} \sum_{i=1}^s \left\{\GammaDist(1/k,k/i)\right\} - \log(s)\right],
\end{equation}
where $\GammaDist(\alpha, \beta)$ is the gamma distribution with shape $\alpha$ and scale $\beta$.
Then, for every $\bm{z}$ we have that $\langle \bm{z}, \tilde{\bm{\theta}}\rangle = \langle \bm{z}, \bm{\theta}\rangle + \epsilon(\bm{z})$ with  $\epsilon(\bm{z}) \sim \GumbelDist(0, \tau)$.
\end{mirrortheorem}

\begin{proof}

Since we perturb each $\theta_i$ by $\varepsilon_i$
 we have, by assumption, that $\langle\bm{\theta}, \mathbf{1}\rangle = k$, for every $\bm{z}$ with $\mathcal{C}(\bm{z})$, that
\begin{equation}
    \langle \bm{z}, \tilde{\bm{\theta}} \rangle = \langle  \bm{z}, \bm{\theta} \rangle + \sum_{j=1}^k \varepsilon_j.
\end{equation}
Since by \cref{lemma-gamma-gumbel} we know that $\sum^{k}_{j=1} \varepsilon_i \sim \GumbelDist(0, \tau)$, the statement of the theorem follows.
\end{proof}

The following  theorem shows that the infinite series from \cref{lemma-gamma-gumbel} can be well approximated by a finite sum using convergence results for the Euler-Mascheroni series.

\begin{theorem}
Let $X\sim \GumbelDist(0, \tau)$ and $\tilde{X}(m) \sim \sum_{j=1}^{\kappa}  \frac{\tau}{\kappa} \left[ \sum_{i=1}^{m} \{ \GammaDist(1/\kappa, \kappa/i)\} - \log(m) \right]$. Then
$$\frac{\tau}{2(m+1)} < \mathbb{E}[\tilde{X}(m)] - \mathbb{E}[X]  < \frac{\tau}{2m}.$$
\end{theorem}

\begin{proof}
We have that
\begin{align*}
& \mathbb{E}[\tilde{X}(m)] = \\
& \quad = \mathbb{E} \left[ \sum_{j=1}^{\kappa}  \frac{\tau}{\kappa} \left[ \sum_{i=1}^{m} \{ \GammaDist(1/\kappa, \kappa/i)\} - \log(m) \right]\right] \\
& \quad = \sum_{i=1}^{m}   \mathbb{E}\left[\GammaDist(1/\kappa, \tau \kappa/i)\right] - \tau \log(m) \\
& \quad = \left[ \sum_{i=1}^{m} \frac{1}{\kappa} \frac{\tau \kappa}{i} - \tau \log(m) \right] \\
& \quad = \tau \left[  \sum_{i=1}^{m} \frac{1}{i} - \log(m) \right].
\end{align*}
If $X\sim \GumbelDist(0, \tau)$, we know that $\mathbb{E}\left[X \right] = \tau \gamma$.
Hence,
\begin{align*}
  \mathbb{E}[\tilde{X}(m)] - \mathbb{E}[X] & = \tau \left[  \sum_{i=1}^{m} \frac{1}{i} - \log(m) \right] -   \tau \gamma\\
& = \tau \left[ \sum_{i=1}^{m} \frac{1}{i} - \log(m) - \gamma  \right].
\end{align*}
The theorem now follows from convergence results of the Euler-Mascheroni series~\cite{MORTICI2010}.
\end{proof}

\begin{mirrorfact}
If one uses
$\ell_H$, then \textsc{I-MLE}\@\xspace with the target distribution
of \cref{eqn-target-q-optimal}
and $\rho(\bepsilon)=\delta_0(\bm{\epsilon})$
is equivalent to
the perceptron-rule estimator of the MLE objective between $p(\bm{z}; h_{\bm{v}}(\hat{\bx}_j))$ and $\hat{\by}_j$.
\end{mirrorfact}

\begin{proof}
Rewriting the definition of the Hamming loss gives us
$$\ell_H(\bm{z}, \bm{y}) = \frac{1}{m}\sum_{i=1}^{m} \left( \bm{z}_i + \bm{y}_i - 2 \bm{z}_i\bm{y}_i \right).$$
Hence, we have that $$\grad{\bm{z}_i}{\ell_H} = \frac{1}{m}\left(1 - 2\bm{y}_i\right).$$
Therefore, $\grad{\bm{z}_i}{\ell_H} = -\frac{1}{m}$ if $\bm{y}_i = 1$ and $\grad{\bm{z}_i}{H} = \frac{1}{m}$ if $\bm{y}_i = 0$.
Since, by definition $\bm{y} \in \mathcal{C}$, we have that
$$\fMAP\left(-\grad{\bm{z}}{\ell_H}\right) = \bm{y}.$$
Now, when using \textsc{I-MLE}\@\xspace with $S=1$ and $\rho(\bepsilon)=\delta_0(\bm{\epsilon})$ we approximate the gradients as
\[\widehat{\nabla}_{\btheta}\mathcal{L}(\bm{\theta}, \bm{\theta}') = \fMAP(\bm{\theta}) - \fMAP(\bm{\theta}') = \fMAP(\bm{\theta}) - \fMAP\left(-\grad{\bm{z}}{\ell_H}\right) = \fMAP(\bm{\theta})  - \bm{y}.\]
This concludes the proof.
\end{proof}

\begin{mirrorfact}
If one uses $\ell_R$
then \textsc{I-MLE}\@\xspace with the target distribution
of \cref{eqn-target-q-optimal} is equivalent to the perturb-and-MAP estimator of the MLE objective between $p(\bm{z}; h_{\bm{v}}(\hat{\bx}_j))$ and $p(\bm{z}; - \hat{\bc}_j)$.
\end{mirrorfact}

\begin{proof}
We have that $\grad{\bm{z}_i}{\ell_R} = \bm{c}_i$ for all $i$.
Now, when using \textsc{I-MLE}\@\xspace with target distribution $\hat{q}(\bm{z}; \bm{\theta}')$ of \cref{eqn-target-q-optimal} (and without loss of generality, for $S=1$) we have that $\hat{q}(\bm{z}; \bm{\theta}') = p(\bm{z}; -\bm{c})$, and we approximate the gradients as
\[ \widehat{\nabla}_{\btheta}\mathcal{L}(\bm{\theta}, \bm{\theta}') =  \fMAP(\bm{\theta} + \bm{\epsilon}_i) - \fMAP(\bm{\theta}' + \bm{\epsilon}_i) = \fMAP(\bm{\theta} + \bm{\epsilon}_i) - \fMAP(-\bm{c} + \bm{\epsilon}_i),  \;  \text{where} \;
    \bm{\epsilon}_i \sim \rho(\bepsilon). \]
Hence, \textsc{I-MLE}\@\xspace approximates the gradients of the maximum likelihood estimation problem between $p(\bm{z}; h_{\bm{v}}(\hat{\bx}_j))$ and $p(\bm{z}; - \hat{\bc}_j)$ using perturb-and-MAP. This concludes the proof.
\end{proof}

\clearpage

\section{Experiments: Details and Additional Results}

\subsection{Synthetic Experiments}
\label{ax-toy-topk}

\begin{figure}[t]
    \centering
    \includegraphics[width=0.32\textwidth]{toy_topk/cm_I-MLE-PaM-SoG.pdf}

    \includegraphics[width=0.32\textwidth]{toy_topk/cm_I-MLE-PaM-Gumbel.pdf}

    \includegraphics[width=0.32\textwidth]{toy_topk/cm_I-MLE-SoG-Gum-means.pdf}
    \caption{Average (over 100 runs) values of $L(\bm{\theta})$ after 50 steps of stochastic gradient descent (with momentum) using single-sample \textsc{I-MLE}\@\xspace with $\SoGDist(1, 5, 10)$ noise (left) and $\GumbelDist(0, 1)$ noise (center) varying the perturbation intensity $\lambda$ (see \cref{eq-pid-q}) and learning rate. The rightmost heat-map depicts the (point-wise) difference between the two methods (blue = better SoG).}
    \label{fig:toy_exps_apx_means}
\end{figure}
\begin{figure}[t]
    \centering
    \includegraphics[width=0.32\textwidth]{toy_topk/cm_I-MLE-PaM-SoG-std.pdf}

    \includegraphics[width=0.32\textwidth]{toy_topk/cm_I-MLE-PaM-Gumbel-std.pdf}

    \includegraphics[width=0.32\textwidth]{toy_topk/cm_I-MLE-SoG-Gum-std.pdf}
    \caption{Same as above, but reporting standard deviations.}
    \label{fig:toy_exps_apx_1}
\end{figure}

In this series of experiments  we analyzed the behaviour of various discrete gradient estimators, comparing our proposed \textsc{I-MLE}\@\xspace with standard straight-trhough (STE) and score-function (SFE) estimators.
We also study of the effect of using Sum-of-Gamma perturbations rather than standard Gumbel noise.
In order to be able to compute exactly (up to numerical precision) all the quantities involved, we chose a  tractable $5$-subset distribution (see Example 2) of size $m=10$.
\begin{wrapfigure}[13]{r}{0.38\textwidth}
\centering
\vspace{-4mm}
\includegraphics[width=0.37\textwidth]{toy_topk/times.pdf}%
\vspace{-4mm}
\caption{\small \label{fig:time} Runtime (mean and standard deviation) for computing $L$ and samples of it, as the dimensionality of the $k$-subset distribution increases (with $k=m/2$).
}
\end{wrapfigure}
We set the loss to $L(\bm{\theta}) = \mathbb{E}_{\hat{\bz} \sim p(\mathbf{z}; \bm{\theta})} [ \norm{\hat{\bz} - \mathbf{b}}^2 ]$, where $\mathbf{b}$ is a fixed vector sampled (only once) from $\mathcal{N}(0,\mathbf{I})$.
This amounts to an unconditional setup where there are no input features (as in the standard MLE setting of \cref{apx-mle}), but where the point-wise loss $\ell(\bm{z})$ is the (Euclidian) distance between the distribution output and a fixed vector $\mathbf{b}$.
In \cref{fig:time} we plot the runtime (mean and standard deviation over 10 evaluations) of the full objective $L$ (Expect., in the plot), of a (faithful) sample of $\ell$ and of a perturb-and-MAP sample with Sum-of-Gamma noise distribution (P\&M) for increasing size $m$, with $k=m/2$.
As it is evident from the plot, the runtime for both expectation and faithful samples, which require computing all the states in $\mathcal{C}$, increases exponentially, while perturb and MAP remains almost constant.

Within this setting, the one-sample \textsc{I-MLE}\@\xspace  estimator is
\begin{equation*}
    \widehat{\nabla}_{\mathrm{I\mbox{-}MLE}} L(\bm{\theta}) = \fMAP(\bm{\theta} + \epsilon) - \fMAP(\bm{\theta}' + \epsilon), \; \text{with} \; \epsilon\sim \SoGDist(1, 5, 10)
\end{equation*}
where $\bm{\theta}'=\bm{\theta} - \lambda [2 ( \hat{\bz}  - \mathbf{b})]$, where   $\hat{\bz}=\fMAP(\bm{\theta} + \epsilon)$ is a (perturb-and-MAP) sample, while the one-sample straight through estimator is
\begin{equation*}
    \widehat{\nabla}_{\mathrm{STE}}{L(\bm{\theta})} =  2 ( \hat{\bz}  - \mathbf{b}), \; \text{with} \; \hat{\bz} = \fMAP(\bm{\theta} + \epsilon), \; \epsilon \sim \GumbelDist(0, 1).
\end{equation*}
For the score function estimator, we have used an expansive faithful sample/full marginal implementation  given by
\begin{equation*}
     \widehat{\nabla}_{\mathrm{SFE}}{L(\bm{\theta})} =   \norm{\hat{\bz} - \mathbf{b}}^2  \grad{\bm{\theta}}{\log p(\hat{\bz}; \bm{\theta})}
     = \norm{\hat{\bz} - \mathbf{b}}^2 [\hat{\bz} - \bm{\mu}(\bm{\theta})], \; \text{with} \; \hat{\bz} \sim p(\bz; \btheta)
\end{equation*}
since, in preliminary experiments, we did not manage to obtain meaningful results with SFE using perturb-and-MAP for sampling and/or marginals approximation.
These equations give the formulae for the estimators which we used for the results plotted in \cref{fig:toy_experiments} (top) in the main paper.

In \cref{fig:toy_exps_apx_means} we plot the heat-maps for the sensitivity results comparing between \textsc{I-MLE}\@\xspace with SoG and \textsc{I-MLE}\@\xspace with Gumbel perturbations.
The two leftmost heat-maps depict the average value (over 100 runs) of $L(\bm{\theta})$ after 50 steps of stochastic gradient descent, for various choices of $\lambda$ and learning rates (momentum factor was fixed at 0.9 for all experiments).
The rightmost plot of \cref{fig:toy_exps_apx_means} is the same as the one in the main paper, and represents the difference between the first and the second heat-maps.
\cref{fig:toy_exps_apx_1} refers to the same setting, but this time showing standard deviations.
The rightmost plot of \cref{fig:toy_exps_apx_1} suggests that using SoG perturbations results also in reduced variance (of the final loss) for most of the tried hyperparameter combinations.
Finally, in \cref{fig:toy_exps_apx_ste_sfe} we show sensitivity plots for STE (both with SoG and Gumbel perturbations) and SFE, where we vary the learning rate.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.19\textwidth]{toy_topk/cm_STE-PaM-SoG.pdf}
    \includegraphics[width=0.19\textwidth]{toy_topk/cm_STE-PaM-Gumbel.pdf}

    \includegraphics[width=0.19\textwidth]{toy_topk/cm_STE-SoG-Gum-means.pdf}

    \includegraphics[width=0.19\textwidth]{toy_topk/cm_SFE-SM-mean.pdf}
    \includegraphics[width=0.19\textwidth]{toy_topk/cm_SFE-SM-std.pdf}
    \caption{First three plots, from  left to right: average (over 100 runs) final values of $L(\bm{\theta})$ after 50 steps of optimization using the straight-through estimator with SoG noise, varying the learning rate; same but using Gumbel noise; difference of averages between the first and the second heat-maps. Last two plots, from left to right: average (over 20 runs) final values of $L(\bm{\theta})$ after 500 steps of optimization using the score function estimator (with faithful samples and exact marginals), varying the learning rate; standard deviation for the same setting.}
    \label{fig:toy_exps_apx_ste_sfe}
\end{figure}

\subsection{Learning to Explain}

\begin{table*}[t]
\small
    \centering
    \resizebox{\textwidth}{!}{
    \begin{tabular}{ccccccc}
    \toprule
    \multirow{2}{*}{\bf Method} & \multicolumn{2}{c}{\bf Appearance} & \multicolumn{2}{c}{\bf Palate} & \multicolumn{2}{c}{\bf Taste} \\

    \cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7}
    & \multicolumn{1}{c}{\bf Test MSE}  & \multicolumn{1}{c}{\bf Subset precision} & \multicolumn{1}{c}{\bf Test MSE}  & \multicolumn{1}{c}{\bf Subset precision} & \multicolumn{1}{c}{\bf Test MSE}  & \multicolumn{1}{c}{\bf Subset precision}  \\

    \midrule
    L2X ($t=0.1$) & 10.70 $\pm$ 4.82 & 30.02 $\pm$ 15.82 & 6.70 $\pm$ 0.63 & 50.39 $\pm$ 13.58 & 6.92 $\pm$ 1.61 & 32.23 $\pm$ 4.92 \\
    SoftSub ($t=0.5$) & \ \textbf{2.48} $\pm$ 0.10 & 52.86 $\pm$ \ \ 7.08 &  \textbf{2.94} $\pm$ 0.08 &  39.17 $\pm$ \ \ 3.17 & \textbf{2.18} $\pm$ 0.10 & \textbf{41.98} $\pm$ 1.42 \\
    \textsc{I-Mle} ($\tau=30$) & \ \textbf{2.51} $\pm$ 0.05 & \textbf{65.47} $\pm$ \ \ 4.95 &  \textbf{2.96} $\pm$ 0.04 & 40.73 $\pm$ \ \ 3.15 & 2.38 $\pm$ 0.04 & \textbf{41.38} $\pm$ 1.55 \\
    \bottomrule
\end{tabular}
}
\caption{\label{table-all-aspects-k-10} Experimental results (mean $\pm$ std. dev.) for the learning to explain experiments for $k=10$ and various aspects.}
\end{table*}

Experiments were run on a server with Intel(R) Xeon(R) CPU E5-2637 v4 @ 3.50GHz, 4 GeForce GTX 1080 Ti, and 128 GB RAM.

The pre-trained word embeddings and data set can be found here: \url{http://people.csail.mit.edu/taolei/beer/}.
Figure~\ref{fig-appendix-l2x} depicts the neural network architecture used for the experiments. As in prior work, we use a batch size of 40. The maximum review length is 350 tokens. We use the standard neural network architecture from prior work~\cite{chen2018learning,paulus2020gradient}. The dimensions of the token embeddings (of the embedding layers) are 200. All 1D convolutional layers have 250 filters with a kernel size of $3$. All dense layers have a dimension of 100. The dropout layer has a dropout rate of 0.2. The layer Multiply perform the multiplication between the token mask (output of I-MLE) and the embedding matrix. The Lambda layer computes the mean of the selected embedding vectors. The last dense layer has a sigmoid activation. IMLESubsetkLayer is the layer implementing \textsc{I-Mle}. We train for 20 epochs using the standard Adam settings in Tensorflow 2.4.1 (learning rate=0.001, beta1=0.9, beta2=0.999, epsilon=1e-07, amsgrad=False), and no learning rate schedule. The training time (for the 20 epochs) for \textsc{I-MLE},  with sum-of-Gamma perturbations, is $380$ seconds, for SoftSub $360$ seconds, and for L2X $340$ seconds. We always evaluate the model with the best validation MSE among the $20$ epochs.

Implementations of \textsc{I-MLE}\@\xspace and all experiments will soon be made available.
\cref{table-all-aspects-k-10} lists the results for L2X, SoftSub, and \textsc{I-Mle} for three additional aromas and $k=10$.

\begin{table}[t!]
\centering
\resizebox{\columnwidth}{!}{
\setlength{\tabcolsep}{2pt}
\begin{tabular}{cccc}
\includegraphics[width=.25\linewidth]{figures/maps/0.pdf} & \includegraphics[width=.25\linewidth]{figures/maps/12.pdf} & \includegraphics[width=.25\linewidth]{figures/maps/2.pdf} & \includegraphics[width=.25\linewidth]{figures/maps/13.pdf}\\
\includegraphics[width=.25\linewidth]{figures/maps/3.pdf} & \includegraphics[width=.25\linewidth]{figures/maps/4.pdf} & \includegraphics[width=.25\linewidth]{figures/maps/18.pdf} & \includegraphics[width=.25\linewidth]{figures/maps/10.pdf}\\
\includegraphics[width=.25\linewidth]{figures/maps/6.pdf} & \includegraphics[width=.25\linewidth]{figures/maps/7.pdf} & \includegraphics[width=.25\linewidth]{figures/maps/8.pdf} & \includegraphics[width=.25\linewidth]{figures/maps/11.pdf}\\
\includegraphics[width=.25\linewidth]{figures/maps/20.pdf} & \includegraphics[width=.25\linewidth]{figures/maps/21.pdf} & \includegraphics[width=.25\linewidth]{figures/maps/22.pdf} & \includegraphics[width=.25\linewidth]{figures/maps/23.pdf}\\
\end{tabular}
}
\caption{\label{app:warcraft} Sample of Warcraft maps, and corresponding shortest paths from the upper left to the lower right corner of the map.}
\end{table}

\subsection{Discrete Variational Auto-Encoder}

\begin{figure}[t]%
\centering
\subfigure{%
\label{fig:original-images}%
\includegraphics[width=0.3\textwidth]{figures/original_numbers.pdf}}%
\hspace{2.0mm}
\subfigure{%
\label{fig:reconstructed-images-gamma}%
\includegraphics[width=0.3\textwidth]{figures/recons_numbers_lambda-1.pdf}}%
\hspace{2.0mm}
\subfigure{%
\label{fig:reconstructed-images-lambda-1}%
\includegraphics[width=0.3\textwidth]{figures/recons_numbers_lambda-10.pdf}}%
\caption{\label{figure-dist-plots-mnist}Original MNIST digits from the test set and their reconstructions using the discrete $10$-subset VAE trained with Sum-of-Gamma perturbations for $\lambda=1$ (center) and $\lambda=10$ (right).}
\end{figure}

Experiments were run on a server with Intel(R) Xeon(R) CPU E5-2637 v4 @ 3.50GHz, 4 GeForce GTX 1080 Ti, and 128 GB RAM.

The data set can be loaded in Tensorflow 2.x with \textit{tf.keras.datasets.mnist.load\_data()}. As in prior work, we use a  batch size of $100$ and train for $100$ epochs, plotting the test loss after each epoch. We use the standard Adam settings in Tensorflow 2.4.1 (learning rate=0.001, beta1=0.9, beta2=0.999, epsilon=1e-07, amsgrad=False), and no learning rate schedule. The MNIST dataset consists in black-and-white $28\times28$ pixels images of hand-written digits. The encoder network consists of an input layer with dimension $784$ (we flatten the images), a dense layer with dimension $512$ and ReLu activation, a dense layer with dimension $256$ and ReLu activation, and a dense layer with dimension $400$ ($20\times20$) which outputs the $\bm{\theta}$ and no non-linearity.
The IMLESubsetkLayer takes $\bm{\theta}$ as input and outputs a discrete latent code of size $20\times20$. The decoder network, which takes this discrete latent code as input, consists of a dense layer with dimension $256$ and ReLu activation, a dense layer with dimension $512$ and ReLu activation, and finally a dense layer with dimension $784$ returning the logits for the output pixels. Sigmoids are applied to these logits and the binary cross-entropy loss is computed.  The training time (for the 100 epochs) was 21 minutes with the sum-of-Gamma perturbations and 18 minutes for the standard Gumbel perturbations.

The code is available in the submission system as a Jupyter notebook.
\begin{figure}[t!]
\centering
\includegraphics[width=1.0\textwidth]{figures/model_plot.pdf}%
\caption{\label{fig-appendix-l2x}The neural network architecture for the learning to explain experiments. (Please zoom into the vector graphic for more details.) We use the standard architecture and settings from prior work~\cite{chen2018learning}. The maximum review length is 350 tokens. The dimensions of the token embeddings (of the embedding layers) are 200. All 1D convolutional layers have 250 filters with a kernel size of $3$. All dense layers have a dimension of 100. The dropout layer has a dropout rate of 0.2. The layer Multiply perform the multiplication between the token mask (output of I-MLE) and the embedding matrix. The Lambda layer computes the mean of the selected embedding vectors The last dense layer has a sigmoid activation. IMLESubsetkLayer is the layer implementing \textsc{I-Mle}. Code is available in the submission system.}
\end{figure}
\begin{figure}[t]
\centering
\includegraphics[width=0.9\columnwidth]{figures/warcraft/warcraft_neurips_12-crop.pdf}
\includegraphics[width=0.9\columnwidth]{figures/warcraft/warcraft_neurips_18-crop.pdf}
\includegraphics[width=0.9\columnwidth]{figures/warcraft/warcraft_neurips_24-crop.pdf}
\includegraphics[width=0.9\columnwidth]{figures/warcraft/warcraft_neurips_30-crop.pdf}
\caption{Training dynamics for different models on $K \times K$ shortest path tasks on Warcraft maps, with $K \in \{ 12, 18, 24, 30 \}$.} \label{app:dynamics}
\end{figure}

\subsection{Differentiating through Combinatorial Solvers} \label{sec-axp:warcraft}

The experiments were run on a server with  Intel(R) Xeon(R) Silver 4208 CPU @ 2.10GHz CPUs, 4 NVIDIA Titan RTX GPUs, and 256 GB main memory.

\cref{app:warcraft} shows a set of $30 \times 30$ Warcraft maps, and the corresponding shortest paths from the upper left to the lower right corner of the map.
In these experiments, we follow the same experimental protocol of \citet{poganvcic2019differentiation}: optimisation was carried out via the Adam optimiser, with scheduled learning rate drops dividing the learning rate by 10 at epochs 30 and 40.
The initial learning rate was $5 \times 10^{-4}$, and the models were trained for 50 epochs using 70 as the batch size.
As in \citep{poganvcic2019differentiation}, the $K \times K$ weights matrix is produced by a subset of the ResNet18~\citep{DBLP:conf/cvpr/HeZRS16} architecture
whose weights are trained on the task.
For training BB, in all experimental results in \cref{section-experiment-blackbox} and \cref{app:warcraft}, the hyperparameter $\lambda$ was set to $\lambda = 20$.

\cref{app:dynamics} shows the training dynamics of different models, including the method proposed by \citet{poganvcic2019differentiation} (BB) with different choices of the $\lambda$ hyperparameter, the ResNet18 baseline proposed by \citet{poganvcic2019differentiation}, and \textsc{I-Mle}.
Code and data for the experiments will soon be made available.

\section{Introduction}

While deep neural networks excel at perceptual tasks, they tend to generalize poorly whenever the problem at hand requires
some level of symbolic manipulation, reasoning, or known algorithmic structure.
Logic, relations, and explanations, as well as decision processes, frequently find natural abstractions in discrete structures, ill-captured by the continuous mappings of standard neural nets.
Several application domains, ranging from relational and explainable ML to (discrete) decision-making~\citep{mivsic2020data}, could benefit from general-purpose learning algorithms whose inductive biases are more amenable to integrating symbolic and neural computation.
Motivated by these considerations, there is a growing interest in end-to-end learnable models incorporating discrete components that allow e.g. to sample from discrete latent distributions~\citep{jang2016categorical,paulus2020gradient} or solve combinatorial optimization problems~\citep{poganvcic2019differentiation,Mandi_Guns:2020}.

For complex discrete distributions, it is intractable to compute the exact gradients of the expected loss.
For combinatorial optimization problems, there are discontinuities, and the gradients are zero
almost everywhere.
The standard approach revolves around problem-specific smooth relaxations, which allow one to fall back to (stochastic) backpropagation.
These strategies, however, require tailor-made relaxations, presuppose access to the constraints and are, therefore, not always feasible nor tractable for large state spaces.
Moreover, reverting to discrete outputs at test time may cause unexpected behavior. In other situations, discrete outputs are required at training time because one has to make one of a number of discrete choices, such as accessing discrete memory or deciding on an action in a game.

With this paper, we take a step towards the vision of general-purpose algorithms for hybrid learning systems.
Specifically, we consider settings where the discrete component(s), embedded in a larger computational graph, are discrete random variables from the (constrained) exponential family\footnote{This includes integer linear programs via a natural link that we outline in \cref{ex:ILP}.}.
Grounded in concepts from Maximum Likelihood Estimation (MLE) and perturbation-based implicit differentiation, we propose Implicit Maximum Likelihood Estimation (\textsc{I-MLE}\@\xspace).
To approximate the gradients of the discrete distributions' parameters, \textsc{I-MLE}\@\xspace computes, at each update step, a \emph{target distribution} $q$ that depends on the loss incurred from the discrete output in the forward pass.
In the backward pass, we approximate maximum likelihood gradients by treating $q$ as the empirical distribution.
We propose ways to derive target distributions and introduce a novel family of noise perturbations well-suited for approximating marginals via perturb-and-MAP.
\textsc{I-MLE}\@\xspace is general-purpose as it only requires the ability to compute most probable states and not faithful samples or probabilistic inference.
We show that \textsc{I-MLE}\@\xspace simplifies to \emph{explicit} maximum-likelihood learning when used in some recently studied learning settings involving combinatorial optimization solvers.  Experimental results suggest that \textsc{I-MLE}\@\xspace is flexible and competitive compared to the straight-through and relaxation-based estimators, while being simple to implement.

\section{Problem Statement and Motivation}
\label{sec-ps}

\begin{wrapfigure}[8]{r}{0.48\textwidth}
\centering
\vspace{-8mm}
\includegraphics[clip,width=0.47\textwidth]{figures/dtexp-family-illustration}%
\caption{\small \label{setting-ilustration-main} Illustration of the addressed learning problem. $\bm{z}$ is the discrete (latent) structure.
}
\end{wrapfigure}

We consider models described by the equations
\begin{equation} \label{eq:hybridm}
 \bm{\theta} = h_{\bm{v}}(\bm{x}), \quad \bm{z}\sim p(\bz; \btheta), \quad  \bm{y}=f_{\bm{u}}(\bm{z}),
\end{equation}
where $\bm{x}\in\mathcal{X}$ and $\bm{y}\in\mathcal{Y}$ denote feature inputs and target outputs,
$h_{\bm{v}}:\mathcal{X}\to\Theta$ and $f_{\bm{u}}:\mathcal{Z}\to\mathcal{Y}$ are smooth parameterized maps,
and $p(\bz; \btheta)$ is a discrete probability distribution.

Given a set of examples $\mathcal{D}=\{(\hat{\bx}_j, \hat{\by}_j)\}_{j=1}^N$, we are concerned with learning the parameters $\bm{\omega} = (\bm{v}, \bm{u})$ of \eqref{eq:hybridm} by finding approximate solutions of
$\min_{\bm{\omega}}  \sum_j L(\hat{\bx}_j, \hat{\by}_j; \bm{\omega})/N$.
The training error $L$ is typically
defined as:
\begin{equation} \label{eq:loss_cs}
L(\hat{\bx}, \hat{\by}; \bm{\omega}) = \mathbb{E}_{\hat{\bz}\sim p(\bm{z}; \hat{\btheta})}\left[\ell(f_{\bm{u}}(\hat{\bz}), \hat{\by})\right] \quad \text{with} \quad \hat{\btheta} = h_{\bm{v}}(\hat{\bx}),
\end{equation}
where $\ell:\mathcal{Y}\times\mathcal{Y}\to\mathbb{R}^+$ is a point-wise loss function. \cref{setting-ilustration-main} illustrates the setting.
For example, an interesting instance
of \eqref{eq:hybridm} and \eqref{eq:loss_cs} arises in
\emph{learning to explain} user reviews~\citep{chen2018learning} where the task is to infer a target sentiment score (e.g.\@\xspace w.r.t. the quality of a product) from a review while also providing a concise explanation of the predicted score by selecting a subset of exactly $k$ words (cf. \cref{ex-topk}). In \cref{sec-exps}, we present experiments precisely in this setting.
As anticipated in the introduction, we
restrict the discussion to instances in which
$p(\bz; \btheta)$ belongs
to the (constrained) \emph{discrete exponential family}, which we now formally introduce.

Let $\bm{Z}$ be a vector of discrete random variables over a state space $\mathcal{Z}$ and let $\mathcal{C} \subseteq \mathcal{Z}$ be the set of states  that satisfy a given set of \emph{linear constraints}.\footnote{For the sake of simplicity we assume $\mathcal{Z} \subseteq \{0, 1\}^{m}$. The set $\mathcal{C}$ is the integral polytope spanned by the problem-specific and given  linear constraints.}
Let $\bm{\theta} \in \Theta \subseteq \mathbb{R}^m$ be a real-valued parameter vector.
The probability mass function (PMF) of a discrete constrained exponential family r.v. is:
\begin{equation} \label{def-constrained-exp-family}
p(\bm{z}; \bm{\theta}) = \left\lbrace
\begin{array}{ll}
     \exp\left(\langle\bm{z},\bm{\theta}\rangle/\tau - A(\bm{\theta}) \right) & \text{if } \bm{z} \in \mathcal{C}, \\
     0 & \text{otherwise.}

\end{array}
\right.
\end{equation}
Here, $\langle\cdot, \cdot\rangle$ is the inner product and $\tau$ the temperature, which, if not mentioned otherwise, is assumed to be $1$.
$A(\bm{\theta})$ is the log-partition function defined as $A(\bm{\theta}) = \log\left(\sum_{\bm{z} \in \mathcal{C}} \exp \left(\langle\bm{z},\bm{\theta}\rangle/\tau  \right)\right)$.
We call $\langle \bm{z}, \bm{\theta}\rangle$ the \emph{weight} of the state $\bm{z}$.
The \emph{marginals} (expected value, mean) of
the
r.v.s
$\mathbf{Z}$ are defined as
$\bm{\mu}(\bm{\theta}) \coloneqq \mathop{\mathbb{E}}_{\hat{\bz} \sim p(\bz; \btheta)}[\hat{\bz}]$.
Finally, the most probable or Maximum A-Posteriori (MAP) states are
defined as $\fMAP(\bm{\theta}) \coloneqq \argmax_{\bm{z} \in \mathcal{C}} \;\langle\bm{z},\bm{\theta}\rangle$.
The family of probability distributions we define here captures a broad range of settings and subsumes probability distributions such as positive Markov random fields and statistical relational formalisms~\citep{wainwright2008graphical,raedt2016statistical}.
We now discuss some examples which we will use in the experiments.
Crucially, in \cref{ex:ILP} we establish the link between the constrained exponential family and
integer linear programming (ILP)
identifying the ILP cost coefficients with the
distribution's parameters $\bm{\theta}$.

\begin{example}[Categorical Variables] \label{ex-catvars}
An $m$-way (one-hot) categorical variable
corresponds to $p(\bm{z}; \bm{\theta}) = \exp\left( \langle \bm{z}, \bm{\theta}\rangle - A(\bm{\theta}) \right)$, subject to the constraint $\langle \bm{z}, \mathbf{1} \rangle = 1$, where $\mathbf{1}$ is a vector of ones.
\end{example}

As
$\mathcal{C} = \{\mathbf{e}_i\}_{i=1}^m$, where $\mathbf{e}_i$ is the $i$-th vector of the canonical base, the parameters of the above distribution coincide with the weights, which are often called \emph{logits} in this context.
The marginals $\bm{\mu}$ coincide with the PMF and can be expressed through a closed-form smooth function of $\bm{\theta}$: the softmax.
This facilitates a natural relaxation that involves using $\bm{\mu}(\bm{\theta})$ in place of $\bm{z}$~\citep{jang2016categorical}.
The
properties
of the categorical distribution, however, quickly disappear
even for slightly more complex distributions, as the following example shows.
\begin{example}[$k$-subset Selection]  \label{ex-topk}
Assume we want to sample binary $m$-dimensional vectors with $k$ ones.
This amounts to replacing the constraint in \cref{ex-catvars} by the constraint
$\langle \bm{z}, \mathbf{1}\rangle = k$.

\end{example}
Here, a closed-form expression for the marginals does not exist:
sampling from this distribution requires computing
the $\binom{m}{k}=O(m^k)$ weights (if $k\leq m/2$). Computing
MAP states instead
takes time linear in $m$.

\begin{example}[%
Integer Linear Programs%
]
\label{ex:ILP}
Consider the combinatorial optimization problem given by the integer linear program $\argmin_{\bm{z} \in \mathcal{C}} \langle \bm{z}, \bm{c} \rangle$, where $\mathcal{C}$ is an integral polytope and $\bm{c}\in\mathbb{R}^m$ is a vector of  cost coefficients, and let $\bm{z}^*(\bm{c})$ be the set of its solutions.
We can associate to the ILP the family (indexed by $\tau > 0$) of probability distributions  $p(\bz; \btheta)$ from  \eqref{def-constrained-exp-family},  with $\mathcal{C}$ the ILP polytope and $\bm{\theta}=-\bm{c}$.
Then, for every $\tau>0$, the solutions of the ILP correspond to the MAP states: $\fMAP(\bm{\theta}) =  \argmax_{\bm{z} \in \mathcal{C}} \langle\bm{z}, \bm{\theta} \rangle = \bm{z}^*(\bm{c})$ and for $\tau\to 0$ one has that $\Pr(\mathbf{Z}\in \bm{z}^*(\bm{c}))\to 1$.
\end{example}

Many problems of practical interest can be expressed as ILPs, such as finding shortest paths,
planning and scheduling problems, and inference in propositional logic.

\section{The Implicit Maximum Likelihood Estimator} \label{section-imle-framework}

In this section, we develop and motivate a family of general-purpose gradient estimators for \cref{eq:loss_cs}
that respect the structure of $\mathcal{C}$
.~\footnote{
The derivations are adaptable to other types of losses defined over the outputs of \cref{eq:hybridm}.
}
Let $(\hat{\bx}, \hat{\by})\in\mathcal{D}$ be a training example and $\hat{\bz} \sim p(\bm{z}; h_{\bm{v}}(\hat{\bx}))$.
The gradient
of $L$ w.r.t. $\bm{u}$ is given by
$
\grad{\bm{u}}{L
(\hat{\bx}, \hat{\by}; \bm{\omega})}
= \mathbb{E}_{\hat{\bz}} [
    \diff{\bm{u}}{f_{\bm{u}}(\hat{\bz})^{\intercal}
    \grad{\bm{y}}{\ell(\bm{y}, \hat{\by})}
    }
]
$ with $\bm{y} = f_{\bm{u}}(\hat{\bz})$,
which may be estimated by drawing one or more samples from $p$.
Regarding $\grad{\bm{v}}{L}$, one has
\begin{equation}
\label{eq-chain-rule-for-w}
    \grad{\bm{v}}{L
(\hat{\bx}, \hat{\by}; \bm{\omega})} =
    \diff{\bm{v}}{h_{\bm{v}}(\hat{\bx})^\intercal}
    \nabla_{\bm{\theta}}L(\hat{\bx}, \hat{\by}; \bm{\omega}),

\end{equation}
where the major challenge is to compute $\nabla_{\bm{\theta}}L$.
A
standard approach is to employ the score function estimator (SFE) which typically suffers from high variance. Whenever a pathwise derivative estimator (PDE) is available it is usually the preferred choice \citep{schulman2015gradient}.
In our setting, however, the PDE is not readily applicable
since $\bm{z}$ is discrete and, therefore, every
(exact) reparameterization path would be discontinuous.
Various authors developed (biased) adaptations of the PDE for discrete r.v.s (see \cref{sec:rw}).
These involve either smooth approximations of $p(\bz; \btheta)$
or approximations of the derivative of the reparameterization map.
Our proposal departs from these two routes and instead involves the formulation
of an implicit maximum likelihood estimation problem.
In a nutshell, I-MLE is a (biased) estimator that replaces
$\nabla_{\bm{\theta}}L$ in \cref{eq-chain-rule-for-w}
with $\widehat{\nabla}_{\btheta}\mathcal{L}$, where $\mathcal{L}$ is an implicitly defined MLE objective and $\widehat{\nabla}$ is an estimator of the gradient.

We now focus on deriving the (implicit) MLE objective $\mathcal{L}$. Let us assume we can, for a given $\hat{\by}$, construct an exponential family distribution $q(\bm{z}; \bm{\theta}')$,
which we call the \emph{target distribution}, such that:
\begin{equation*}
    \mathbb{E}_{\hat{\bz}\sim q(\bm{z}; \bm{\theta}')}\left[\ell(f_{\bm{u}}(\hat{\bz}), \hat{\by})\right] \leq \mathbb{E}_{\hat{\bz}\sim p(\bm{z}; \bm{\theta})}\left[\ell(f_{\bm{u}}(\hat{\bz}), \hat{\by})\right].
\end{equation*}
We then define
$\mathcal{L}(\bm{\theta}, \bm{\theta}') \coloneqq - \mathbb{E}_{\hat{\bz}\sim q(\bm{z}; \bm{\theta}')} [\log p(\hat{\bz}; \bm{\theta})]
=
\mathbb{E}_{\hat{\bz}\sim q(\bm{z}; \bm{\theta}')} [
    A(\bm{\theta}) - \langle \hat{\bz}, \bm{\theta} \rangle
]
$.
$\mathcal{L}$ can be interpreted\footnote{
We expand on this in \ifbool{printApdx}{\cref{apx-mle}}{the appendix} where we also review the classic MLE setup \cite[][Ch. 9]{murphy2012machine}.
} as an MLE objective between the (current) model distribution $p$ and the target distribution
$q$.
The intuition is that making $p$
more similar to $q$
will also reduce the loss $L(\hat{\bx}, \hat{\by}; \bm{\omega})=\mathbb{E}_{\hat{\bz}\sim p(\bm{z}; \bm{\theta})}\left[\ell(f_{\bm{u}}(\hat{\bz}), \hat{\by})\right]$.
Exploiting the fact that $\nabla_{\bm{\theta}} A(\bm{\theta}) = \bm{\mu}(\bm{\theta})$, the gradient of $\mathcal{L}$ is computed as:
\begin{equation} \label{eq-grad-mleobj}
    \grad{\bm{\theta}}{\mathcal{L}(\bm{\theta}, \bm{\theta}')} =
    \bm{\mu}(\bm{\theta}) - \mathbb{E}_{\hat{\bz}\sim q(\bm{z}; \bm{\theta}')} [\hat{\bz}]
    =
    \bm{\mu}(\bm{\theta}) - \bm{\mu}(\bm{\theta}'),
\end{equation}
that is%
, the difference between the marginals of the current distribution $p$ and the marginals of the  target distribution $q$,
also equivalent to the gradient of the KL divergence between $p$ and $q$.

We will not use \cref{eq-grad-mleobj} directly, as computing the marginals is, in general, a \#P-hard problem and scales poorly with the dimensionality $m$.
MAP states are typically less expensive to compute (e.g.\@\xspace see \cref{ex-topk}) and are often used to approximate $\bm{\mu}(\bm{\theta})$\footnote{This is known as the \emph{perceptron learning rule} in standard MLE.}
or to compute
perturb-and-MAP approximations,
where $\bm{\mu}(\bm{\theta})\approx \mathbb{E}_{\bm{\epsilon}\sim \rho(\bepsilon)}\fMAP(\bm{\theta} + \bm{\epsilon})$ where $\bm{\epsilon}\sim\rho(\bepsilon)$ is an appropriate \emph{noise distribution} with domain $\mathbb{R}^{m}$.
In this work we follow -- and explore in more detail in \cref{section-perturb-and-map} -- the latter approach (also referred to as the Gumbel-max trick \citep[cf.][]{Papandreou:2011}),
which
retains most of the computational  advantages of the MAP approximation but may be far less crude.
Henceforth, we only assume access to an algorithm to compute MAP states (this may be a standard ILP solver in the case of \cref{ex:ILP}) and rephrase, with a slight abuse of notation, \cref{eq:hybridm} as
\begin{equation} \label{eq:hybridm2}
 \bm{\theta} = h_{\bm{v}}(\bm{x}), \quad
 \bm{z}
 =
 \fMAP(\bm{\theta} + \bm{\epsilon})
 \; \text{with} \; \bm{\epsilon} \sim p(\bm{\epsilon}), \quad  \bm{y}=f_{\bm{u}}(\bm{z}).
\end{equation}
With \cref{eq:hybridm2} in place,
the general expression for the I-MLE estimator is $\widehat{\nabla}_{\bm{v}} L(\bm{x}, \bm{y}; \bm{\omega})=    \diff{\bm{v}}{h_{\bm{v}}(\hat{\bx})^\intercal}
    \widehat{\nabla}_{\btheta}\mathcal{L}(\bm{\theta}, \bm{\theta}')$
with $\bm{\theta}=h_{\bm{v}}(\hat{\bx})$ where, for $S\in\mathbb{N}^+$:
\begin{equation}
\label{eqn-perturb-map-approx}
    \widehat{\nabla}_{\btheta}\mathcal{L}(\bm{\theta}, \bm{\theta}') =
    \frac{1}{S}\sum_{i=1}^S
    [
        \fMAP(\bm{\theta} + \bm{\epsilon}_i) -
        \fMAP(\bm{\theta}' + \bm{\epsilon}_i)
    ],  \;  \text{with} \;
    \bm{\epsilon}_i \sim \rho(\bepsilon) \; \text{for} \; i\in\{1, \dots, S\}.
\end{equation}
Note that, if the states of both $p$ and $q$ are binary vectors,
$\widehat{\nabla}_{\btheta}\mathcal{L}(\bm{\theta}, \bm{\theta}') \in [-1, 1]^{m}$ and when $S=1$ $\widehat{\nabla}_{\btheta}\mathcal{L}(\bm{\theta}, \bm{\theta}') \in \{-1, 0, 1\}^{m}$.
A pseudo-code implementation of I-MLE is available in \ifbool{printApdx}{\cref{app:algorithm}}{the appendix}.
In the following, we discuss the problem of constructing the target distributions $q$.

\subsection{Target Distributions via Perturbation-based Implicit Differentiation}
\label{subsection-implicit-differentiation}

The efficacy of the I-MLE estimator hinges on a
proper
choice of $q$, a hyperparameter of our framework.
In this section we derive and motivate a class of general-purpose target distributions, rooted in perturbation-based implicit differentiation (PID):
\begin{equation}
    \label{eq-pid-q}
    q(\bm{z}; \bm{\theta}') = p(\bm{z}; \bm{\theta} - \lambda \grad{\bm{z}}{\ell}(f_{\bm{u}}(\bm{z}), \hat{\by}))
    \; \text{with} \;
    \bm{z} = \fMAP(\bm{\theta} + \bm{\epsilon}) \; \text{and} \;
    \bm{\epsilon}\sim \rho(\bepsilon),
\end{equation}
where $(\hat{\bx}, \hat{\by})\in\mathcal{D}$ is a data point,  $\bm{\theta}=h_{\bm{v}}(\hat{\bx})$, and $\lambda>0$ is a hyperparameter that controls the perturbation intensity.

To motivate \cref{eq-pid-q}, consider the setting where the inputs to $f$ are the marginals of $p(\bz; \btheta)$ (rather than discrete perturb-and-MAP samples as in \cref{eq:hybridm2}), that is, $\bm{y} = f_{\bm{u}}(\bm{\mu}(\bm{\theta}))$ with $\bm{\theta}=h_{\bm{v}}(\hat{\bx})$, and redefine the training error $L$ of \cref{eq:loss_cs} accordingly.
A seminal result by \cite{domke:2010} shows that, in this case, we can obtain  $\grad{\bm{\theta}}{L}$ by perturbation-based differentiation as:
\begin{equation}
\label{eqn-domke}
\grad{\bm{\theta}}{L
(\hat{\bx}, \hat{\by}; \bm{\omega})} = \lim_{\lambda \rightarrow
0}
\left\lbrace\frac{1}{\lambda}
 \left[ \bm{\mu}(\bm{\theta}) - \bm{\mu}\left(\bm{\theta} - \lambda  \grad{\bm{\mu}}{
L(\hat{\bx}, \hat{\by}; \bm{\omega})
} \right)\right]\right\rbrace,
\end{equation}
where
$\grad{\bm{\mu}}{L}
=
\diff{\bm{\mu}}{f_{\bm{u}}(\bm{\mu})}^\intercal \grad{\bm{y}}{\ell(\bm{y}, \hat{\by})}
$.
The expression inside the limit may be interpreted as the gradient of an implicit MLE objective
(see \cref{eq-grad-mleobj})  between the distribution $p$ with (current) parameters $\bm{\theta}$  and
$p$
with parameters perturbed in the negative direction of the downstream gradient $\grad{\bm{\mu}}{L}$.
Now, we can adapt \eqref{eqn-domke} to our setting of \cref{eq:hybridm2} by resorting to the straight-through estimator (STE) assumption \citep{bengio2013estimating}.
Here, the STE assumption translates into reparameterizing $\bm{z}$ as a function of $\bm{\mu}$ and approximating $\diff{\bm{\mu}}{\bm{z}}\approx \bm{I}$.
Then, $\grad{\bm{\mu}}{L}=\diff{\bm{\mu}}{\bm{z}}^\intercal \grad{\bm{z}}{L}\approx \grad{\bm{z}}{L}$ and we approximate \cref{eqn-domke} as
\begin{equation}
\label{eqn-approx-domke-stright-through}
\grad{\bm{\theta}}{L
(\hat{\bx}, \hat{\by}; \bm{\omega})} \approx
\frac{1}{\lambda}
\left[ \bm{\mu}(\bm{\theta}) - \bm{\mu}\left(\bm{\theta} - \lambda \grad{\bm{z}}{L}\right)\right]
 = \frac{1}{\lambda}\nabla \mathcal{L}_{\bm{\theta}}(\bm{\theta}, \bm{\theta} - \lambda \grad{\bm{z}}{L}),
\end{equation}
for some $\lambda>0$.
From \cref{eqn-approx-domke-stright-through} we derive \eqref{eq-pid-q} by taking a single sample estimator of $\grad{\bm{z}}{L}$
(with perturb-and-MAP sampling)
and by incorporating the constant $1/\lambda$ into a global learning rate.
I-MLE with PID target distributions may be seen as a way to generalize the STE to more complex distributions.
Instead of using the gradients $\grad{\bm{z}}{L}$ to backpropagate directly,
\textsc{I-MLE}\@\xspace uses them to construct a target distribution $q$.
With that, it defines an implicit maximum likelihood objective, whose gradient (estimator) propagates the supervisory signal upstream,
critically, taking the constraints into account.
When using \cref{eq-pid-q} with $\rho(\bepsilon)=\delta_{0}(\bm{\epsilon})$\footnote{$\delta_{0}$ is the Dirac delta centered around $0$ -- this is equivalent to approximating the marginals with $\fMAP$.},
the I-MLE estimator also recovers a recently proposed gradient estimation rule to differentiate through black-box combinatorial optimization problems~\citep{poganvcic2019differentiation}. \textsc{I-MLE}\@\xspace unifies existing gradient estimation rules in one framework.

We will resume the discussion about target distributions in \cref{section-optimal-q-co}, where we analyze more closely the setup in \cref{ex:ILP}. Next, we focus on the perturb-and-MAP strategies and derive a class of noise distributions that are particularly apt to the settings we consider in this work.

\subsection{A Novel Family of Perturb-and-MAP Noise Distributions}
\label{section-perturb-and-map}

We propose a novel and general way to design tailored perturb and MAP sampling strategies. While the proposed family of noise distributions works with \textsc{I-MLE}, the results of this Section are of independent interest and can also be used in other (relaxed) perturb-and-MAP based gradient estimators~\citep{paulus2020gradient}.
The following proposition is a generalization of the classic perturb and MAP result by \citet{Papandreou:2011} (also known as the Gumbel-max trick) in that it considers the temperature $\tau$.

\begin{proposition}
\label{prob-full-potential}
Let $p(\bm{z}; \bm{\theta})$ be a discrete exponential family distribution with integer polytope $\mathcal{C}$ and temperature $\tau$, and let $\langle\bm{z},\bm{\theta}\rangle$ be the unnormalized weight of each $\bm{z} \in \mathcal{C}$. Moreover, let $\tilde{\bm{\theta}}$ be such that, for all $\bm{z} \in \mathcal{C}$, $\langle\bm{z},\tilde{\bm{\theta}}\rangle =  \langle\bm{z},\bm{\theta}\rangle + \epsilon(\bm{z})$ with each
$\epsilon(\bm{z})$ sampled i.i.d. from $\GumbelDist(0,\tau)$. Then we have that  $\Pr(\fMAP(\tilde{\bm{\theta}})=\bm{z}) = p(\bm{z}; \bm{\theta}).$
\end{proposition}

\ifbool{proofsIMLE}{
    \begin{proof}
    Let $\epsilon_i \sim \GumbelDist(0,\tau)$ i.i.d. and $\tilde{\theta}_i = \theta_i + \epsilon_i$. Following a derivation similar to one made in \citet{Papandreou:2011}, we have:
    %
    \begin{align*}
        & \Pr\{\argmax(\tilde{\theta}_1, \ldots, \tilde{\theta}_m) = n\}  = \\
        = & \Pr\{ \tilde{\theta}_n \geq \max_{j\neq n} \{\tilde{\theta}_j\} \}  \\
        = & \int_{-\infty}^{+\infty} g(t; \theta_n) \prod_{j\neq n} G(t; \theta_j) \,dt \\
        = & \int_{-\infty}^{+\infty} \frac{1}{\tau} \exp\left(\frac{\theta_n - t}{\tau} - e^{\frac{\theta_n - t}{\tau}}\right) \prod_{j\neq n} \exp\left( -e^{\frac{\theta_j - t}{\tau}}\right) \,dt \\
        = &  \int_{-\infty}^{+\infty} \frac{1}{\tau} e^{\frac{\theta_n - t}{\tau}}\exp\left(-e^{\frac{\theta_n-t}{\tau}}\right) \prod_{j\neq n} \exp\left(-e^{\frac{\theta_j - t}{\tau}}\right) \,dt \\
        = & \int_{0}^{1} \prod_{j \neq n} z^{\exp\left( \frac{\theta_j - \theta_n}{\tau}  \right)} \,dz \qquad \text{with } z := \exp\left( -e^{\frac{\theta_n - t}{\tau}} \right) \\
        = & \frac{1}{1 + \sum_{j \neq n} e^{\frac{\theta_j - \theta_n}{\tau}} } \\
        = & \frac{e^{\frac{\theta_n}{\tau}}}{\sum_{j=1}^{m} e^{\frac{\theta_j}{\tau}}, }
      \end{align*}
      where $g$ and $G$ are respectively the Gumbel probability density function and the Gumbel cumulative density function.
      The proposition now follows from arguments made in \citet{Papandreou:2011} using the maximal equivalent re-parameterization of $p(\bm{z}; \bm{\theta})$ where we specify a parameter $\langle \bm{\theta},\bm{z}\rangle$ for each $\bm{z} \in \mathcal{C}$ and perturb these parameters.
    \end{proof}
}{
All proofs can be found in
\ifbool{printApdx}{\cref{proofs-imle}.}{the appendix.}
}
The proposition states that if we can perturb the weights $\langle \bm{z},\bm{\theta} \rangle$ of each $\bm{z} \in \mathcal{C}$ with independent $\GumbelDist(0, \tau)$ noise, then obtaining MAP states from the perturbed model is equivalent to sampling from $p(\bm{z}; \bm{\theta})$ at temperature\footnote{
Note that the temperature here is different to the temperature of the Gumbel softmax trick~\citep{jang2016categorical}
which scales \emph{both} the sum of the logits \emph{and} the samples from $\GumbelDist(0, 1)$.}
 $\tau$.
For general constrained exponential distributions,  perturbing the weights $\langle\bm{z},\bm{\theta}\rangle$ for each state $\bm{z} \in \mathcal{C}$ is at least as expensive as computing the marginals exactly. Hence, one usually \emph{locally} perturbs each $[\bm{\theta}]_i$ (the $i$-th component of the vector $\bm{\theta}$) with Gumbel noise.
Fortunately, we can prove that,  for a large class of problems, it is possible to design more suitable \emph{local} perturbations.
First, we show that, for any $\kappa \in \mathbb{N}^+$, a Gumbel distribution can be written as a finite sum of $\kappa$
i.i.d. (implicitly defined) random variables.
\begin{lemma}
\label{lemma-gamma-gumbel}
Let $X \sim \GumbelDist(0 ,\tau)$ and let $\kappa \in \mathbb{N}^+$.
Define the Sum-of-Gamma distribution as
\begin{equation}
    \label{eqn-sum-of-gammas}
    \SoGDist(\kappa, \tau, s) \coloneqq \frac{\tau}{\kappa}

    \left\lbrace\sum_{i=1}^{s} \left\{\GammaDist(1/\kappa, \kappa/i)\right\} - \log(s)\right\rbrace,
\end{equation}
where $s\in\mathbb{N}^+$ and $\GammaDist(\alpha, \beta)$ is the Gamma distribution with shape $\alpha$ and scale $\beta$, and let $\SoGDist(\kappa, \tau)\coloneqq\lim_{s \rightarrow \infty} \SoGDist(\kappa, \tau, s)$.
Then we have that $X\sim \sum_{j=1}^{\kappa} \bm{\epsilon}_j$, with $\bm{\epsilon}_j\sim \SoGDist(\kappa, \tau).$
\end{lemma}

\ifbool{proofsIMLE}{
    \begin{proof}
    Let $\kappa \in \mathbb{N}^+$ and let $X \sim \GumbelDist(0, \tau)$. Its moment generating function has the form
    \begin{equation}
    \label{eq:m-gen-f-gumb}
    \mathbb{E}[\exp(tX)] = \Gamma(1 - \tau t).
    \end{equation}
    As mentioned in \citet[p. 443,][]{johnson1998advances} we know that we can write the Gamma function as
    \begin{equation}
    \label{eq-gamma-book}
    \Gamma(1 - \tau t) = e^{\gamma \tau t} \prod_{i=1}^{\infty} \left(1 - \frac{\tau t}{i}\right)^{-1} e^{\frac{-\tau t}{i}}
    \end{equation}
    where $\gamma$ is the Euler-Mascheroni constant.
    We have that
    $$\left(1 - \frac{\tau t}{i}\right)^{-1} = \frac{i}{i - \tau t} = \frac{\frac{i}{\tau}}{\frac{i - \tau t}{\tau}} =  \frac{\frac{i}{\tau}}{\frac{i}{\tau} - \frac{\tau t}{\tau}} =
    \frac{\frac{i}{\tau}}{\frac{i}{\tau} - t}.$$
    The last term is the moment generating function of an exponential distribution with scale $\frac{\tau}{i}$.

    We can now take the logarithm on both sides of \cref{eq-gamma-book} and obtain
    \begin{equation*}
      tX = \gamma\tau t + \lim_{s\rightarrow \infty}\sum_{i=1}^{s} \left( t \ExpDist(\tau/i) - \frac{\tau t}{i} \right),
    \end{equation*}
    where $\ExpDist(\alpha)$ is the exponential distribution with scale $\alpha$.
    Hence,
    \begin{align*}
     X \sim & \lim_{s\rightarrow \infty}\sum_{i=1}^{s}  \left( \ExpDist(\tau/i) - \frac{\tau}{i} \right) +  \gamma\tau \\
      = & \lim_{s\rightarrow \infty}\sum_{i=1}^{s}  \left( \ExpDist(\tau/i) - \frac{\tau}{i} \right) + \tau \lim_{s\rightarrow \infty}\sum_{i=1}^{s} \frac{1}{i} - \log(s) \\
      = & \lim_{s\rightarrow \infty}\sum_{i=1}^{s}  \left( \ExpDist(\tau/i) - \frac{\tau}{i} + \frac{\tau}{i} \right) - \tau\log(s)   \\
        = & \lim_{s\rightarrow \infty}\sum_{i=1}^{s}  \ExpDist(\tau/i) - \tau\log(s)
    \end{align*}

    Since $\ExpDist(\alpha) \sim \GammaDist(1, \alpha)$, and due to the scaling and summation properties of the Gamma distribution (with shape-scale parameterization), we can write for all $r>1$:
    %
    $$\ExpDist(\alpha) \sim \sum_{j=1}^{r} \GammaDist(1/r, \alpha r)/r.$$
    %

    Hence, picking $r=\kappa$ from the hypothesis, we have

    \begin{align*}
        X \sim & \lim_{s\rightarrow \infty} \left\lbrace
            \sum_{i=1}^{s} \sum_{j=1}^{\kappa} \GammaDist(1/\kappa, \tau \kappa/i)/\kappa
        \right\rbrace - \tau \log(s) \\
        = &
        \lim_{s\rightarrow \infty} \left\lbrace
            \sum_{j=1}^{\kappa} \frac{\tau}{k}\sum_{i=1}^{s} \GammaDist(1/\kappa, \kappa/i)
        \right\rbrace - \sum_{j=1}^\kappa \frac{\tau}{\kappa} \log(s) \\
        = &
        \lim_{s\rightarrow \infty}
        \sum_{j=1}^{\kappa} \frac{\tau}{\kappa}
        \left\lbrace
            \left[
                \sum_{i=1}^{s} \GammaDist(1/\kappa, \kappa/i)
            \right] - \log(s)
        \right\rbrace \\
        = &
        \sum_{j=1}^{\kappa}

        \frac{\tau}{\kappa}
        \left\lbrace
            \lim_{s\rightarrow \infty}
            \left[
                \sum_{i=1}^{s} \GammaDist(1/\kappa, \kappa/i)
            \right] - \log(s)
        \right\rbrace

    \end{align*}
    This concludes the proof.
    Parts of the proof are inspired by a post on stackexchange~\cite{214875}.
    \end{proof}
}

Based on \cref{lemma-gamma-gumbel}, we can show that for exponential family distributions where every $\bm{z} \in \mathcal{C}$ has exactly $k$ non-zero entries
we can design perturbations of $\langle\bm{z},\bm{\theta}\rangle$ following a  Gumbel distribution.

\begin{theorem}
\label{thm-perturb-distribution}
Let $p(\bm{z}; \bm{\theta})$ be a discrete exponential family distribution with integer polytope $\mathcal{C}$ and temperature $\tau$.
Assume that if $\bm{z} \in \mathcal{C}$ then $\langle\bm{z}, \mathbf{1}\rangle = k$ for some constant $k \in \mathbb{N}^+$.
Let $\tilde{\bm{\theta}}$ be the perturbation obtained by $[\tilde{\bm{\theta}}]_j = [\bm{\theta}]_j + \bm{\epsilon}_j$ with $\bm{\epsilon}_j \sim \SoGDist(k, \tau)$ from \cref{eqn-sum-of-gammas}.
Then,
$\forall\bm{z}\in\mathcal{C}$ we have that $\langle \bm{z}, \tilde{\bm{\theta}}\rangle = \langle \bm{z}, \bm{\theta}\rangle + \epsilon(\bm{z})$, with $\epsilon(\bm{z}) \sim \GumbelDist(0, \tau)$.
\end{theorem}

\ifbool{proofsIMLE}{
    \begin{proof}

    Since we perturb each $\theta_i$ by $\varepsilon_i$

     we have, by assumption, that $\langle\bm{\theta}, \mathbf{1}\rangle = k$, for every $\bm{z} \in \mathcal{C}$, that
    \begin{equation}
        \langle \bm{z}, \tilde{\bm{\theta}} \rangle = \langle  \bm{z}, \bm{\theta} \rangle + \sum_{j=1}^k \varepsilon_j.
    \end{equation}
    Since by \cref{lemma-gamma-gumbel} we know that $\sum^{k}_{j=1} \varepsilon_i \sim \GumbelDist(0, \tau)$, the statement of the theorem follows.
    \end{proof}
}

\begin{wrapfigure}[10]{r}{0.489\textwidth}
\centering
\vspace{-5.5mm}
\subfigure{\includegraphics[scale=0.24]{figures/gumbel-k.pdf}}
\subfigure{\includegraphics[scale=0.24]{figures/gamma-k.pdf}}
\vspace{-4mm}
\caption{\small \label{figure-dist-plots} Histograms for $10$k samples where each sample is (left) the sum of $5$
$\epsilon_j\sim\GumbelDist(0, 1)$ or (right) the sum of $5$ $\epsilon_j\sim\SoGDist(5, 1, 10)$.
}
\end{wrapfigure}
Many problems such as $k$-subset selection, traveling salesman, spanning tree, and graph matching strictly satisfy the assumption of \cref{thm-perturb-distribution}.
We can, however, also apply the strategy in cases where the variance of
$\langle \mathbf{Z}, \mathbf{1} \rangle$
is small (e.g.\@\xspace shortest weighted path).
The Sum-of-Gamma perturbations provide a more fine-grained approach to  noise perturbations. For $\tau=\kappa=1$, we obtain the standard Gumbel perturbations.
In contrast to the standard $\GumbelDist(0, 1)$ noise, the proposed local Sum-of-Gamma
perturbations result in
weights' perturbations
that follow the Gumbel distribution. \cref{figure-dist-plots} shows histograms of $10$k samples, where each sample is either the sum of $5$ samples from $\GumbelDist(0, 1)$ (the standard approach) or the sum of $k=5$ samples from $\SoGDist(5, 1, 10) = \frac{1}{5}\sum_{i=1}^{10} \{\GammaDist(1/5, 5/i) - \log(10)\}$.
While we still cannot sample faithfully from $p(\bm{z}; \bm{\theta})$ as the perturbations are \emph{not} independent, we can counteract the problem of partially dependent perturbations by increasing the temperature $\tau$ and, therefore, the variance of
the noise distribution.
We explore and verify the importance of tuning $\tau$ empirically.
In the appendix, we also show that the infinite series from \cref{lemma-gamma-gumbel} can be well approximated by a finite sum using convergence results for the Euler-Mascheroni series~\citep{MORTICI2010}.

\section{Target Distributions for Combinatorial Optimization Problems}
\label{section-optimal-q-co}

In this section, we explore the setting where the discrete computational component arises from a combinatorial optimization (CO) problem, specifically an integer linear program (ILP).
Many authors have recently  considered
the setup where the CO component occupies the last layer of the model defined by \cref{eq:hybridm} (where $f_{\bm{u}}$ is the identity) and the supervision is available  in terms of examples of either optimal solutions~\citep[e.g.][]{poganvcic2019differentiation} or optimal cost coefficients (conditioned on the inputs)~\citep[e.g.][]{elmachtoub2020smart}.
We have seen in  \cref{ex:ILP} that we can naturally associate
to each ILP a probability distribution  \eqref{def-constrained-exp-family} with $\bm{\theta}$ given by the negative  cost coefficients $\bm{c}$ of the ILP and $\mathcal{C}$ the integral polytope. Letting $\tau \rightarrow 0$ is equivalent to taking the MAP in the forward pass. Furthermore, in \cref{subsection-implicit-differentiation} we showed that the I-MLE framework subsumes a recently propose method by  \citet{poganvcic2019differentiation}.
Here, instead, we show that, for a certain choice of the target distribution, \textsc{I-MLE}\@\xspace estimates the gradient of an explicit maximum likelihood learning
loss $\mathcal{L}$ where the data distribution is ascribed to either (examples of) optimal solutions or optimal cost coefficients.

Let $q(\bm{z}; \bm{\theta}')$ be the distribution $p(\bm{z}; \bm{\theta}')$, with parameters
\begin{equation}
\label{eqn-target-q-optimal}
[{\bm{\theta}'}]_i :=\left\lbrace
\begin{array}{ll}
     [\bm{\theta}]_i  & \mbox{ if } [\grad{\bm{z}}{L}]_i = 0  \\
     - [\grad{\bm{z}}{L}]_i & \mbox{ otherwise}.
\end{array}
\right.
\end{equation}
In the first CO setting, we observe training data $\mathcal{D}=\{(\hat{\bx}_j, \hat{\by}_j)\}_{j=1}^N$ where $\hat{\by}_j\in\mathcal{C}$ and the loss $\ell$ measures a distance between a discrete  $\hat{\bz}_j \sim p(\bm{z}; \bm{\theta}_j)$ with $\bm{\theta}_j = h_{v}(\hat{\bx})$ and
a given optimal solution of the ILP
$\hat{\by}_j$.
An example is the Hamming loss $\ell_H$~\citep{poganvcic2019differentiation} defined as $\ell_H(\bm{z}, \bm{y}) = \bm{z} \circ (\mathbf{1} - \bm{y}) + \bm{y} \circ (\mathbf{1} - \bm{z})$, where $\circ$ denotes the Hadamard (or entry-wise) product.

\begin{fact}
\label{thm-optimal-target-distribution-y}
If one uses
$\ell_H$, then \textsc{I-MLE}\@\xspace with the target distribution
of \cref{eqn-target-q-optimal}
and $\rho(\bepsilon)=\delta_0(\bm{\epsilon})$
is equivalent to
the perceptron-rule estimator of the MLE objective between $p(\bm{z}; h_{\bm{v}}(\hat{\bx}_j))$ and $\hat{\by}_j$.
\end{fact}

\ifbool{proofsIMLE}{
\begin{proof}
Rewriting the definition of the Hamming loss gives us
$$\ell_H(\bm{z}, \bm{y}) = \frac{1}{m}\sum_{i=1}^{m} \left( \bm{z}_i + \bm{y}_i - 2 \bm{z}_i\bm{y}_i \right).$$
Hence, we have that $$\grad{\bm{z}_i}{\ell_H} = \frac{1}{m}\left(1 - 2\bm{y}_i\right).$$
Therefore, $\grad{\bm{z}_i}{\ell_H} = -\frac{1}{m}$ if $\bm{y}_i = 1$ and $\grad{\bm{z}_i}{H} = \frac{1}{m}$ if $\bm{y}_i = 0$.
Since, by definition $\bm{y} \in \mathcal{C}$, we have that
$$\fMAP\left(-\grad{\bm{z}}{\ell_H}\right) = \bm{y}.$$
Now, when using \textsc{I-MLE}\@\xspace with $S=1$ and $\rho(\bepsilon)=\delta_0(\bm{\epsilon})$ we approximate the gradients as
\[\widehat{\nabla}_{\btheta}\mathcal{L}(\bm{\theta}, \bm{\theta}') = \fMAP(\bm{\theta}) - \fMAP(\bm{\theta}') = \fMAP(\bm{\theta}) - \fMAP\left(-\grad{\bm{z}}{\ell_H}\right) = \fMAP(\bm{\theta})  - \bm{y}.\]
This concludes the proof.
\end{proof}
}
It follows that the method by \citet{poganvcic2019differentiation} returns, for a large enough $\lambda$, the maximum-likelihood gradients  (scaled by $1/\lambda$) approximated by the perceptron rule.
Proofs in \ifbool{printApdx}{\cref{proofs-imle}.}{the Appendix.}

In the second CO setting, we observe training data $\mathcal{D}=\{(\hat{\bx}_j, \hat{\bm{c}}_j)\}_{j=1}^N$, where $\hat{\bm{c}}_j$ is the optimal cost conditioned on input $\hat{\bx}_j$. Here,
various authors \citep[e.g.][]{elmachtoub2020smart,Mandi_Guns:2020,DBLP:conf/nips/MandiG20} use as point-wise loss
the regret $\ell_R(\bm{\theta}, \bm{c}) = \bm{c}^{\top} \left(\bm{z}(\bm{\theta}) - \hat{\bz}^{*}(\bm{c})\right)$ where $\bm{z}(\bm{\theta})$ is a state  sampled from $p(\bm{z}; \bm{\theta})$ (possibly with temperature $\tau\to0$, that is, a MAP state) and $\hat{\bz}^{*}(\bm{c})\in \bm{z}^*(\bm{c})$ is an optimal state for $\bm{c}$.

\begin{fact}
\label{thm-optimal-target-distribution-c}
If one uses $\ell_R$
then \textsc{I-MLE}\@\xspace with the target distribution
of \cref{eqn-target-q-optimal} is equivalent to the perturb-and-MAP estimator of the MLE objective between $p(\bm{z}; h_{\bm{v}}(\hat{\bx}_j))$ and $p(\bm{z}; - \hat{\bc}_j)$.
\end{fact}

\ifbool{proofsIMLE}{
\begin{proof}
We have that $\grad{\bm{z}_i}{\ell_R} = \bm{c}_i$ for all $i$.
Now, when using \textsc{I-MLE}\@\xspace with target distribution $\hat{q}(\bm{z}; \bm{\theta}')$ of \cref{eqn-target-q-optimal} (and without loss of generality, for $S=1$) we have that $\hat{q}(\bm{z}; \bm{\theta}') = p(\bm{z}; -\bm{c})$, and we approximate the gradients as
\[ \widehat{\nabla}_{\btheta}\mathcal{L}(\bm{\theta}, \bm{\theta}') =  \fMAP(\bm{\theta} + \bm{\epsilon}_i) - \fMAP(\bm{\theta}' + \bm{\epsilon}_i) = \fMAP(\bm{\theta} + \bm{\epsilon}_i) - \fMAP(-\bm{c} + \bm{\epsilon}_i),  \;  \text{where} \;
    \bm{\epsilon}_i \sim \rho(\bepsilon). \]
Hence, \textsc{I-MLE}\@\xspace approximates the gradients of the maximum likelihood estimation problem between $p(\bm{z}; h_{\bm{v}}(\hat{\bx}_j))$ and $p(\bm{z}; - \hat{\bc}_j)$ using perturb-and-MAP. This concludes the proof.
\end{proof}
}

This last result also implies that when using the target distribution $q$ from \eqref{eqn-target-q-optimal} in conjunction with the regret, \textsc{I-MLE}\@\xspace performs maximum-likelihood learning minimizing the KL divergence between the current distribution and the distribution whose parameters are the optimal cost.

\section{Related Work}
\label{sec:rw}

Several papers address the gradient estimation problem for discrete r.v.s, many resorting to relaxations.
\citet{maddison2016concrete, jang2016categorical} propose the Gumbel-softmax distribution to relax categorical r.v.s;
\citet{paulus2020gradient} study extensions to more complex probability distributions.
\citet{tucker2017rebar, grathwohl2017backpropagation} develop parameterized control variates based on continuous relaxations for the score-function estimator.
In contrast to the above, we focus explicitly on problems where \emph{only} discrete samples are used during training.
Approaches that do not rely on relaxations are specific to certain distributions~\citep{bengio2013estimating, franceschi2019learning, liu2019rao} or assume knowledge of $\mathcal{C}$~\citep{kool2020estimating}.
We provide a general-purpose framework that does not require access to the linear constraints and the corresponding integer polytope $\mathcal{C}$.
Experiments in the next section show that while \textsc{I-MLE}\@\xspace only requires a MAP solver, it is competitive and sometimes outperforms tailor-made relaxations.
SparseMAP~\citep{Niculae2018SparseMAPDS} is an approach to structured prediction and latent variables, replacing the exponential distribution (specifically, the softmax) with a sparser distribution. Similar to our work, it only presupposes the availability of a MAP oracle.
LP-SparseMAP~\citep{Niculae2020LPSparseMAPDR} is an extension that uses a relaxation of the optimization problem rather than a MAP solve. Sparsity can also be exploited for efficient marginal inference in latent variable models~\citep{correia2020efficient}.

A series of works
about differentiating through CO problems
\citep{wilder2019melding, elmachtoub2020smart, ferber2020mipaal, DBLP:conf/nips/MandiG20} relax ILPs by adding $L^1$, $L^2$ or log-barrier regularization terms and differentiate through the KKT conditions deriving from the application of the cutting plane or the interior-point methods.
These approaches are conceptually linked to techniques for differentiating through smooth programs \citep{amos2017optnet, donti2017task, agrawal2019differentiable, chen2020understanding, domke2012generic, franceschi2018bilevel} that arise not only in modelling but also in hyperparameter optimization and meta-learning.  \citet{poganvcic2019differentiation, rolinek2020deep, berthet2020learning} propose methods that are not tied to a specific ILP solver.
As we saw
above,
the former two, originally derived from a continuous interpolation argument, may be interpreted as special instantiations of \textsc{I-MLE}\@\xspace. The latter addresses the theory of perturbed optimizers and discusses perturb and MAP in the context of the Fenchel-Young loss.
All the CO-related works assume that either optimal costs or solutions are given as training data, while \textsc{I-MLE}\@\xspace may be also applied in the absence of such supervision by making use of implicitly generated target distributions.
Other authors focus on devising differentiable relaxations for specific CO problems such as SAT~\citep{evans2018learning} or MaxSAT~\citep{wang2019satnet}.
Machine learning intersects with CO also in other contexts, e.g.\@\xspace in learning heuristics to improve the performances of CO solvers or  differentiable models such as GNNs to \textquote{replace} them; see \citet{bengio2020machine} and references therein.

Direct loss minimization (DLM)~\citep{mcallester2010direct, song2016training} is also related to our work, but the assumption there is that examples of optimal states $\hat{\bz}$ are given. \citet{lorberbom:2019} extend the DLM framework to discrete VAEs using coupled perturbations.
Their approach is tailored to VAEs and not general-purpose.
Under a methodological viewpoint,
\textsc{I-MLE}\@\xspace inherits from classical MLE~\citep{wainwright2008graphical} and perturb-and-MAP~\citep{Papandreou:2011}.
The theory of perturb-and-MAP was used to derive general-purpose upper bounds for log-partition functions~\citep{hazan2012partition,shpakova:2016}.

\section{Experiments}
\label{sec-exps}

The set of experiments can be divided into three parts. First, we analyze and compare the behavior of \textsc{I-MLE}\@\xspace with (i) the score function and (ii) the straight-through estimator using a toy problem. Second, we explore the latent variable setting where both $h_{\bm{v}}$ and $f_{\bm{u}}$ in \cref{eq:hybridm} are neural networks and the optimal structure is \emph{not} available during training.
Finally, we address the problem of differentiating through black-box combinatorial optimization problems. Here, we use the target distribution derived in \cref{section-optimal-q-co}. More experimental details for available in the  appendix.

\paragraph{Synthetic Experiments.}

We conducted a series of experiments with a tractable $5$-subset distribution (see Example 2) where $\mathbf{z}\in\{0,1\}^{10}$.
We set the loss to $L(\bm{\theta}) = \mathbb{E}_{\hat{\bz} \sim p(\mathbf{z}; \bm{\theta})} [ \norm{\hat{\bz} - \mathbf{b}}^2 ]$, where $\mathbf{b}$ is a fixed vector sampled from $\mathcal{N}(0,\mathbf{I})$.
In \cref{fig:toy_experiments} (Top), we plot optimization curves with means and standard deviations, comparing the proposed estimator with the straight-through (STE) and the score function (SFE) estimators.~\footnote{Hyperparameters are optimized against $L$ for all methods independently. Statistics are over 100 runs. We found STE slightly better with Gumbel rather than SoG noise. SFE failed with all tested PaM strategies.}
For STE and I-MLE, we use perturb and MAP (PaM) with Gumbel and $\SoGDist(1, 5, 10)$ noise, respectively.
The SFE uses faithful samples and exact marginals (which is feasible only when $m$ is very small) and converges much more slowly than the other methods, while the STE converges to worse solutions than those found using \textsc{I-MLE}\@\xspace.
Figure~\ref{fig:toy_experiments} (Bottom) shows the benefits of using SoG rather than Gumbel perturbations with \textsc{I-MLE}\@\xspace.
While the best configurations for both are comparable, SoG noise achieves in average (over 100 runs) strictly better final values of $L$ for more than $50\%$ of the tested  configurations (varying $\lambda$ from \cref{eq-pid-q} and the learning rate) and exhibit smaller variance (see \ifbool{printApdx}{\cref{fig:toy_exps_apx_1}}{Appendix}). Additional details and  results in
\ifbool{printApdx}{\cref{ax-toy-topk}.}{the Appendix.}

\begin{table}[t!]
  \begin{minipage}[b]{0.37\linewidth}
  \centering
     \includegraphics[width=0.78\columnwidth]{toy_topk/ToyTopK_wSFE_final_2.pdf}
    \\
    \includegraphics[width=0.9\columnwidth]{toy_topk/cm_I-MLE-SoG-Gum-means.pdf}
    \captionof{figure}{\small \label{fig:toy_experiments}Top:
    Gradient-based optimization of $L$ with various
    estimators.

    Bottom: Mean difference of the final value of $L$ between I-MLE with SoG or Gumbel $\rho(\epsilon)$, varying $\lambda$ and the learning rate (blue $=$ better SoG).

    }
  \end{minipage}
  \hspace{2mm}
  \begin{varwidth}[b]{0.57\linewidth}
  \small
    \begin{tabular}{rrrrr}
    \toprule
    \multicolumn{1}{c}{\multirow{3}{*}{\bf Method}} & \multicolumn{2}{c}{\bf Test MSE}  & \multicolumn{2}{c}{\bf Subset Precision}  \\
    \cmidrule(lr){2-3} \cmidrule(lr){4-5}
    & \multicolumn{1}{c}{\bf Mean} & \multicolumn{1}{c}{\bf Std. Dev.} &  {\bf Mean} &  {\bf Std. Dev.} \\

    \midrule
    \multicolumn{5}{c}{$k=10$} \\
    \midrule
    L2X ($t=0.1$) & 6.68 & 1.08 & 26.65 & 9.39 \\

    SoftSub ($t=0.5$) & \textbf{2.67} & 0.14 & 44.44 & 2.27 \\

    STE  ($\tau=30$) & 4.44 & 0.09 & 38.93 & 0.14 \\
    \textsc{I-MLE}\@\xspace $\fMAP$ & 4.08 & 0.91 & 14.55 & 0.04 \\
    \textsc{I-MLE}\@\xspace $\GumbelDist$ & \textbf{2.68} & 0.10 & 39.28 & 2.62 \\

    \textsc{I-MLE}\@\xspace ($\tau=30$) & \textbf{2.71} & 0.10 & \textbf{47.98} & 2.26 \\
    \midrule
    \multicolumn{5}{c}{$k=5$} \\

    \midrule
    L2X ($t=0.1$) & 5.75 & 0.30 & 33.63 & 6.91 \\
    SoftSub ($t=0.5$) & \textbf{2.57} & 0.12 & \textbf{54.06} & 6.29 \\
    \textsc{I-MLE}\@\xspace ($\tau=5$) & \textbf{2.62} & 0.05 & \textbf{54.76} & 2.50 \\
    \midrule
    \multicolumn{5}{c}{$k=15$} \\
    \midrule
    L2X ($t=0.1$) & 7.71 & 0.64 & 23.49 & 10.93 \\
    SoftSub ($t=0.5$) & \textbf{2.52} & 0.07 & 37.78 & 1.71 \\
    \textsc{I-MLE}\@\xspace ($\tau=30$) & 2.91 & 0.18 & \textbf{39.56} & 2.07 \\
    \bottomrule
    \end{tabular}
    \vspace{4mm}
    \caption{\small \label{tab:l2x-aroma}Detailed results for the aspect \textsc{Aroma}. Test MSE and subset precision, both $\times 100$, for $k \in \{5, 10, 15\}$.}
  \end{varwidth}%
\vspace{-7mm}
\end{table}

\paragraph{Learning to Explain.}
\label{section-experiments-l2x}

The \textsc{BeerAdvocate} dataset \citep{McAuley2012LearningAA} consists of free-text reviews and ratings for $4$ different aspects of beer: appearance, aroma, palate, and taste. Each sentence in the test set has annotations providing the words that best describe the various aspects. Following the experimental setup of recent work~\citep{paulus2020gradient}, we address the problem introduced by the L2X paper~\citep{chen2018learning} of learning a distribution over $k$-subsets of words that best explain a given aspect rating.
The training set has 80k reviews for the aspect \textsc{Appearance} and 70k reviews for all other aspects.  Since the original dataset~\citep{McAuley2012LearningAA} did not provide separate validation and test sets, we compute $10$ different evenly sized validation/test splits of the 10k held out set and compute mean and standard deviation over $10$ models, each trained on one split.
Subset precision was computed using a subset of 993 annotated reviews.
We use pre-trained word embeddings from~\cite{Lei2016RationalizingNP}.
Prior work used non-standard neural networks for which an implementation is not available~\citep{paulus2020gradient}.
Instead, we used the neural network from the L2X paper with $4$ convolutional and one dense layer.
This neural network outputs the parameters $\bm{\theta}$ of the distribution $p(\bm{z}; \bm{\theta})$ over $k$-hot binary latent masks with $k \in \{5, 10, 15\}$.
We compare to relaxation-based baselines L2X~\citep{chen2018learning} and SoftSub~\citep{Xie2019ReparameterizableSS}.
We also compare the straight-through estimator (STE) with Sum-of-Gamma (SoG) perturbations.
We used the standard hyperparameter settings of \citet{chen2018learning} and choose the temperature parameter $t \in \{0.1, 0.5, 1.0, 2.0\}$. For \textsc{I-MLE}\@\xspace we choose $\lambda \in \{10^1, 10^2, 10^3\}$, while for both \textsc{I-MLE}\@\xspace and STE we choose $\tau \in \{k, 2k, 3k\}$ based on the validation MSE.
We used the standard Adam settings.
We trained separate models for each aspect using MSE as point-wise loss $\ell$.

\cref{tab:l2x-aroma} lists detailed results for the aspect \textsc{Aroma}.
\textsc{I-MLE}\@\xspace's MSE values are competitive with those of the best baseline, and its subset precision is significantly higher than all other methods (for $\tau=30$).
Using only MAP as the approximation of the marginals leads to poor results. This shows that using the tailored perturbations with tuned temperature is crucial to achieve state of the art results. The Sum-of-Gamma perturbation introduced in this paper outperforms the standard local Gumbel perturbations. %
More details and results can be found in the appendix.

\paragraph{Discrete Variational Auto-Encoder.}
\label{section-experiment-vae}

We evaluate various perturbation strategies for a discrete $k$-subset Variational Auto-Encoder (VAE) and compare them to the straight-through estimator (STE) and the Gumbel-softmax trick. The latent variables model a probability distribution over $k$-subsets of (or top-$k$ assignments too) binary vectors of length $20$.  The special case of $k=1$ is equivalent to a categorical variable with $20$  categories. For $k>1$, we use \textsc{I-MLE}\@\xspace using the class of PID target distributions of \cref{eq-pid-q} and compare various perturb and MAP noise sampling strategies.
The experimental setup is similar to those used in prior work on the Gumbel softmax tricks~\citep{jang2016categorical}. The loss is the sum of the reconstruction loss (binary cross-entropy loss on output pixels) and the KL divergence between the marginals of the variables and the uniform distribution. The encoding and decoding functions of the VAE consist of three dense layers (encoding: 512-256-20x20; decoding: 256-512-784). We do not use temperature annealing. Using \cref{eqn-perturb-map-approx} with $S=1$, we use either $\GumbelDist(0, 1)$ perturbations (the standard approach)\footnote{Increasing the temperature $\tau$ of $\GumbelDist(0, \tau)$ samples increased the test loss.} or Sum-of-Gamma (SoG) perturbations at a temperature of $\tau=10$. We run $100$ epochs and record the loss on the test data. The difference in training time is negligible.
\cref{fig-autoencoder-plots} shows that using the SoG noise distribution is beneficial.
The test loss using the SoG perturbations is lower despite the perturbations having higher variance and, therefore, samples of the model being more diverse.
This shows that using weights' perturbations that follow a
proper Gumbel distribution is indeed beneficial.
\textsc{I-MLE}\@\xspace significantly outperforms the STE, which does not work in this setting and is competitive with the Gumbel-Softmax trick for the $1$-subset (categorical) distribution where marginals can be computed in closed form.

\begin{figure*}[t!]%
\centering
\subfigure{%
\label{fig:r1}%
\includegraphics[width=0.332\textwidth]{figures/plot-subset-k-autoencoder_r10.pdf}}%
\subfigure{%
\label{fig:r10}%
\includegraphics[width=0.332\textwidth]{figures/plot-subset-k-autoencoder-ste.pdf}}%
\subfigure{%
\label{fig:softmax}%
\includegraphics[width=0.332\textwidth]{figures/plot-subset-1-autoencoder.pdf}}%
\caption{\label{fig-autoencoder-plots} \small Plots of the sum of the binary reconstruction loss and the KL divergence as a function of the number of epochs (lower is better). (Left) Discrete $10$-subset VAE trained with \textsc{I-MLE}\@\xspace with $\lambda=10$ (\textsc{I-MLE}\@\xspace). (Center) Discrete $10$-subset VAE trained with the straight-through estimator (STE). (Right) Discrete $1$-subset VAE using the Gumbel softmax trick (GSMT). The down-up-down artifact is due to temperature annealing. Sum-of-Gamma (SoG) perturbations have the lowest test loss for the $10$-subset VAEs. For $\lambda=10$ and SoG perturbations, the test loss is similar to that of the categorical ($1$-subset) VAE trained with the Gumbel softmax trick.
}
\vspace{-5mm}
\end{figure*}

\paragraph{Differentiating through Combinatorial Solvers.} \label{section-experiment-blackbox}

\begin{wraptable}[10]{r}{0.64\textwidth}
\vspace{-4.4mm}
\small
	\caption{\small Results for the Warcraft shortest path task. Reported is the accuracy, i.e.\@\xspace percentage of paths with the optimal costs. Standard deviations are over five runs.} \label{tab:sp-results}
	\centering
	\begin{tabular}{cccccc}
	\toprule
	$K$ & {\bf \textsc{I-Mle} ($\bm{\mu}$-$\bm{\mu}$)} & {\bf \textsc{I-Mle} ($\mathtt{M}$-$\mathtt{M}$)} & {\bf BB} & {\bf DPO} \\
	\midrule
	12 & ${\bf 97.2}\std{0.5}$ & $95.2\std{0.3}$ & $95.2\std{0.7}$ & $94.8\std{0.3}$ \\
	18 & ${\bf 95.8}\std{0.7}$ & $94.4\std{0.5}$ & $94.7\std{0.4}$ & $92.3\std{0.8}$ \\
	24 & ${\bf 94.3}\std{1.0}$ & $93.2\std{0.2}$ & $93.8\std{0.3}$ & $91.5\std{0.4}$ \\
	30 & ${93.6}\std{0.4}$ & ${\bf 93.7}\std{0.6}$ & $93.6\std{0.5}$ & $91.5\std{0.8}$ \\
	\bottomrule
	\end{tabular}
\end{wraptable}
In these experiments, proposed by \citet{poganvcic2019differentiation}, the training datasets consists of 10,000 examples of randomly generated images of terrain maps from the Warcraft II tile set~\citep{warcraft_map_editor}.
Each example has an underlying $K \times K$ grid whose cells represent terrains with a fixed cost.
The shortest (minimum cost) path between the top-left and bottom-right cell in the grid is encoded as an indicator matrix and serves as the target output.
An image of the terrain map is presented to a CNN, which produces a $K \times K$ matrix of vertex costs.
These costs are then given to \emph{Dijkstra's algorithm} (the MAP solver) to compute the  shortest path.
We closely follow the evaluation protocol of \citet{poganvcic2019differentiation}.
We considered two instantiations of \textsc{I-MLE}\@\xspace: one  derived from \cref{thm-optimal-target-distribution-y} ($\mathtt{M}$-$\mathtt{M}$ in \cref{tab:sp-results}) using $\ell_H$ and one derived from \cref{thm-optimal-target-distribution-c}  ($\bm{\mu}$-$\bm{\mu}$) using $\ell_R$, with $\rho(\epsilon)=\SoGDist(k, 1, 10)$ where $k$ is the empirical mean of the path lengths (different for each grid size $K$).
We compare with the method proposed by \citet{poganvcic2019differentiation} (BB\footnote{
Note that this is the same as using \textsc{I-MLE}\@\xspace with PID target distribution form \cref{eq-pid-q} and $\rho(\epsilon)=\delta_0$.}) and \citet{berthet2020learning} (DPO).
The results are listed in \cref{tab:sp-results}.
\textsc{I-MLE}\@\xspace obtains results comparable to (BB) with  $\mathtt{M}$-$\mathtt{M}$
and is more accurate with  $\bm{\mu}$-$\bm{\mu}$.
We believe that the $\bm{\mu}$-$\bm{\mu}$ advantage may be partially due to an implicit form of data augmentation since we know from \cref{thm-optimal-target-distribution-c} that, by using \textsc{I-MLE}\@\xspace, we obtain samples from the distribution whose parameters are the optimal cost.
Training dynamics, showing faster convergence of \textsc{I-MLE}\@\xspace ($\bm{\mu}$-$\bm{\mu}$), and additional details are available in \ifbool{printApdx}{\cref{app:warcraft}}{the Appendix}.

\begin{comment}
\input{tables/energy_v2}

\subsubsection{Energy-Cost Aware Scheduling}
This is a resource-constrained day-ahead job scheduling problem to minimise total energy cost~\citep{csplib:prob059}.
Tasks must be assigned to a given number of machines, where each task has a duration, an earliest start, a latest end, a resource requirement, and a power usage, and each machine has a resource capacity constraint.
Also, tasks cannot be interrupted once started, nor migrated to another machine, and must be completed before midnight.
Time is discretised in 48 time-slots.
As day-ahead energy prices are not known, they must be predicted for each time-slot first.
The energy price data is taken from the Irish Single Electricity Market Operator~\citep[SEMO,][]{DBLP:conf/cp/IfrimOS12} and it consists of historical energy price data at 30-minute intervals from 2011-2013.
Out of the available 789 days, 552 are used for training, 60 for validation, and 177 for testing.
Each time-slot instance has calendar attributes; day-ahead estimates of weather characteristics; SEMO day-ahead forecasted energy-load, wind-energy production and prices; and actual wind-speed, temperature, CO$_{2}$ intensity and price.
Of the actual attributes, we keep only the actual price and use it as a target for prediction.

Following \citet{DBLP:conf/nips/MandiG20}, we use two kinds of neural networks, one without a hidden layer (0-layer) and one with a single hidden layer (1-layer), and compare different approaches in \cref{tab:exp-energy}.
\end{comment}

\section{Conclusions}

\textsc{I-MLE}\@\xspace is an efficient, simple-to-implement, and general-purpose framework for learning
hybrid
models.
\textsc{I-MLE}\@\xspace is competitive with relaxation-based approaches for discrete latent-variable models and with approaches to backpropagate through CO solvers. Moreover, we showed empirically that \textsc{I-MLE}\@\xspace outperforms the straight-through estimator. A limitation of the work is its dependency on computing MAP states which is, in general, an NP-hard problem (although for many interesting cases there are efficient algorithms).
Future work includes devising target distributions when $\grad{\bm{z}}{L}$ is not available, studying
the properties (including the bias) of the proposed estimator, developing adaptive strategies for $\tau$ and $\lambda$, and integrating and testing \textsc{I-MLE}\@\xspace in several challenging application domains.

redpajama_set_name: RedPajamaArXiv


around 1.25km to the south of the village of West Firle.
has a mound 8m in diameter and 0.3m high, also with a slight, central hollow.
but survive as buried features around 1.5m wide.
practices over of a period of around 3000 years.


text:
Domestic Bonds
Commodity Costs
Is the pandemic remodeling boom looming? Refinancing of withdrawals reaches highest level in 13 years - Forbes Advisor
A year later, Nigeria still awaits a Chinese loan of 5.3 billion dollars for the railway to Kano
70-year-old woman with $ 134,000 in student loan debt trapped in student loan scam
Federal Reserve Expands Main Street Loan Program | Cole schotz
Jota's return could lead to a 4-2-3-1 – Predicting Liverpool's formation against Fulham – Liverpool FC
PPP loan terms pose a challenge for businesses
Would you rather be rich than poor?
Water footprint and socio-economic implications
Highlights of the UNGA Event Potential Sources of Funding for Road Safety | News | SDG Knowledge Center
How Chairs can support positive changes in teaching and learning (opinion)
'Impressive recovery' for Israel's diamond industry
Media: Sanya Science City speeds up creation of Hainan Free Trade Port in 2021 – Business & Economy
HAIKOU / China /, December 1. / TASS /. Yazhouwan Innovation Development Zone spurred the establishment of the Hainan Free Trade Port this year, becoming an important factor in increasing the gross regional p...
Reimagining the future of Sri Lanka's garment industry – adaderana biz english
By Shirendra Lawrence The emergence of clothing and textiles as a major contributor to Sri Lanka's economy began after independence in 1948, when a few pioneering i...
"The Higher Council of Free Zones is seeking the creation of an international stock exchange in free zones"
TEHRAN - The secretary of the High Council of Free Zones of Iran said: "We are seeking the establishment of an international stock exchange in the free zones (FTAs) of the country to help the stock market." ...
Special independent economic zones to steer the battery and electric vehicle value chain in the DRC
Kinshasa, November 25, 2021 (CEA) - Members of the government of the Democratic Republic of Congo (DRC) and leaders of the private sector agreed that for an efficient ecosystem for ...
Another look at the prospects for a Eurasian digital platform
The prophecy should come true, but I will not explain it. Then, in due course, it will be called back and picked up. He who has ears to hear, let him hear. -FM DostoyevskyIn November 2021,...
Nigeria: NEPZA employee convicted of naira 342 million fraud
The Independent Commission on Corrupt Practices and Other Related Offenses (ICPC) obtained the conviction of Mr. Dasel Nanjwan, an employee of the Nigeria Export Processing Zones Authority (NEPZA). Nanjwan, a...
"A favorable business environment will help attract FDI"
Speakers say at an event of the Chamber of Commerce and Industry of foreign investorsFeatured activity report Sun Nov 21, 2021 6:22 PM Last update: Sun Nov 21, 2021 7:17 PMThe Fo...
Bepza seeks more South Korean investments
Featured activity report Thu Nov 18, 2021 6:49 PM Last update: Thu Nov 18, 2021 7:04 PM ...
Government wants to ensure tax breaks for free ports are not wasted
EENTER THE Tilbury Harbor and you'll see building materials piled up, imported, processed and waiting to be transported to construction sites. With London nearby and an increase in demand for home extensions l...
Will Bangladesh fall victim to the Chinese debt trap? – Cashmere reader
Bangladesh, the second largest recipient of Chinese investment in South Asia after Pakistan, imports the largest volume of goods from China, making it Bangladesh's largest trading partner. This is the beginn...

redpajama_set_name: RedPajamaCommonCrawl


text:
September Songs – The Music of Kurt Weill is a music video of 94 minutes recorded on VHS in 1994 for Rhombus Media, ZDF (Germany), CBC (Canada) and RTP (Portugal). It was produced and directed by Larry Weinstein, and written by Weinstein and David Mortin. The film was conceived as a follow-up to the album Lost in the Stars: The Music of Kurt Weill whose producer Hal Willner served as the music supervisor in this project. The film was nominated for the 1995 Emmy Awards for Outstanding Cultural Program; it won five Gemini Awards in 1997. An album was released in 1997.

Background
Larry Weinstein was inspired to film September Songs after he heard Lost in the Stars: The Music of Kurt Weill. The producer of Lost in the Stars Hal Willner was engaged to supervise the recording.
The backdrop of the performance consists of an abandoned warehouse in Toronto where musicians perform a series of songs by Kurt Weill. Between songs, parts of Weill's biography are narrated.

The cast of consists of the Brodsky Quartet, William S. Burroughs, Betty Carter, Nick Cave and Spanish Fly, Elvis Costello, Kathy Dalton, Bob Dorough, Charlie Haden, PJ Harvey, David Johansen, Lou Reed, Mary Margaret O'Hara, The Persuasions, Stan Ridgway, Ralph Schuckett, Ellen Shipley, Teresa Stratas, the Brodsky Quartet, the Y Chamber Symphony. It also features the voices of Lotte Lenya, Bertolt Brecht, and Weill.

Two artists reprised their performances from Lost in the Stars: The Music of Kurt Weill: Lou Reed recorded a second version of "September Song", and Charlie Haden re-recorded "Speak Low" with new arrangement by Fred Hersch and archived voice of Kurt Weill singing the song added. Betty Carter was originally slated to appear in the first album, but scheduling issues precluded her inclusion. Here she performed "Lonely House", a song which she also added to her album I'm Yours, You're Mine. Elvis Costello and Brodsky Quartet were asked to record "Lost in the Stars" as Willner had seen them performing the song as an encore while on tour.

The original video recording was released in 1994 and was shown as part of Great Performances on PBS. An album was released in 1997; many of the performances in the album were recorded live in Toronto, but some were recorded elsewhere. Although Teresa Stratas performed the two songs in the film, "Youkali Tango" and "Surabaya Johnny" were licensed for inclusion in this album from her own albums but with a little accordion added to "Youkali Tango". Also included in the album are two historical performances: a 1955 recording of Weill's wife Lotte Lenya singing "Pirate Jenny", and a 1930 recording of Bertolt Brecht reciting "Mack the Knife". However, no single could be released from the album due to licensing issues.

Track list
The order of songs is taken from the Sony Classical CD. The order in the film is different and there are more songs in the film: Kathy Dalton sings "Aggie's Sewing Machine" ("Aggie's Song") from Johnny Johnson, Stan Ridgway sings the "Cannon Song" from The Three Penny Opera with Tom and Bruce Fowler and a choreography by Janice Hladki, Ghettoriginal Dance Company performs "Mandalay-B-Boy-Parlay" inspired by the "Mandalay Song" from Happy End. The film also features the Esprit Orchestra playing the "Nocturne" and "Weill Variations" (by Alex Pauk) from Weill's incidental music for the play Konjunktur by Léo Lania and the instrumental "Song of the Goddess" from Johnny Johnson played by the MGM Studio Orchestra.

References

External links

 September Songs – The Music of Kurt Weill, cast and crew, fandango.com
Review by Dennis Harvey, Variety, 7 August 1995
Review by Stephanie Zacharek, salon.com, 15 September 1997
Review by Parry Gettelman, Orlando Sentinel, 19 September 1997
Review, Time Out

Music video compilation albums
Kurt Weill tribute albums
1994 musical films
1994 films
Canadian musical films
German musical films
Portuguese musical films
ZDF original programming
CBC Television original films
Rádio e Televisão de Portugal original programming
1990s Canadian films
1990s German films

redpajama_set_name: RedPajamaWikipedia


data on the St. Lawrence global system.
Changing the way scientists work and think in a complex world.
Unlock the Power of Data Integration.
Presence, abundance and distribution of living species in the St. Lawrence River.
Collaborative data acquisition tools for the general public.