equ[][]

# Direct measurement of the Husimi-Q function of the electric-field in the
time-domain

Sho Onoe<EMAIL_ADDRESS>_femtoQ_ Lab, Department of Engineering
Physics, Polytechnique Montréal, Montréal, QC H3T 1J4, Canada  Stéphane
Virally  _femtoQ_ Lab, Department of Engineering Physics, Polytechnique
Montréal, Montréal, QC H3T 1J4, Canada  Denis V. Seletskiy  _femtoQ_ Lab,
Department of Engineering Physics, Polytechnique Montréal, Montréal, QC H3T
1J4, Canada

###### Abstract

We develop the theoretical tools necessary to promote electro-optic sampling
to a time-domain quantum tomography technique. Our proposed framework
implements detection of the time evolution of both the electric-field of a
propagating electromagnetic wave and its Hilbert transform (quadrature).
Direct detection of either quadrature is not strictly possible in the time-
domain, detection efficiency approaching zero when an exact mode-matching to
either quadrature is reached. As all real signals have a limited bandwidth, we
can trace out the irrelevant sampling bandwidth to optimize the detection
efficiency while preserving quantum information of the relevant signal.
Through the developed understanding of the mode structure of the amplitude and
Hilbert transform quadratures, we propose multiplexing and mode-matching
operations on the gating function to extract full quantum information on both
quantities, simultaneously. The proposed methology is poised to open a novel
path toward quantum state tomography and quantum spectroscopy directly in the
time domain.

##  I Introduction

Light can be used to detect and transport quantum information. In order to
reveal this information, various quantum metrological techniques must be
applied. All information can in fact be recovered via quantum-state tomography
of the incoming light. The photocounting theory of Glauber, Kelley, and
Kleiner [1, 2], was one of the first attempts at a quantum mechanical
interpretation of the detection of the radiation field. However, even with
unit quantum efficiency, this technique cannot extract the full quantum
information of mixed or even pure states due to its insensitivity to
phase/quadratures. The nominal work by Yuen and Shapiro [3] introduced the
quantum mechanical interpretation of homodyne detection, which can detect both
quadrature fields and their correlations at a desired frequency $\omega_{0}$.
This interpretation has promoted our understanding of quantum field theory
from a particle- to a field-based point-of-view incorporating the important
role of the phase/quadratures.

Although its statistics can be utilized to reproduce the Wigner function [4,
5, 6], allowing full quantum-state reconstruction of the desired frequency,
experimental implementations and extraction of quantum information have been
met with difficulty due to the presence of noise. A leap in experimental
progress was made with the introduction of the balanced-homodyne detection
scheme [7], capable of suppressing the technical noise in the self-homodyne
method. Rigorous analysis of the (quantum) noise associated with the balanced
homodyne detection [8, 9, 10, 11] together with the introduction of the
spectral decomposition method [12] paved the way toward quantum sensing in the
frequency domain. These research led to the realization that efficient
frequency domain homodyne detection of wavepackets require an understanding of
the space-time character of the signal, to which the detection must be mode-
matched. This insight led to the use of a quasi-monochromatic probe with a
duration that is matched to the temporal duration of the signal, allowing
efficient extraction of quantum information [13, 14] from the state-under-
study. Detection of non-classical states, starting with squeezed vacuum states
[10, 15, 16, 17], and the reproduction of their Wigner function [18]
established balanced homodyne detection as one of the most reliable quantum
sensing techniques in modern physics [19, 20, 21].

Despite its success in the (quasi-)frequency domain, homodyne detection has
not been successfully implemented for quasi-time-domain measurements [22]. For
a monochromatic light, the quadratures are well-defined, and the transition
between the two (e.g. $sin$ and $cos$) can be achieved via swap of parity or a
simple time-delay, typically implemented in experiments via a change of group-
velocity phase. In the time-domain, we can no longer swap between the two
quadratures via a time-delay [23]: this simply changes the detection time and
does not affect the parity of the measurement. As a result, the technical
hurdle towards time-domain measurement is not only a broadening of the probe
spectrum, but we must also correctly identify the quadrature fields in the
quasi-time-domain, and seek the optimal spectral phases and amplitudes to
mode-match the probe to those quadratures.

Electro-optic sampling (EOS)[24, 25, 26, 27, 28] has the ability, via
nonlinearity, to isolate the electric-field from the Hilbert field (the
quadrature orthogonal to the electric-field) at sub-cycle resolution. These
features established EOS as one of the most reliable classical time-domain
sensing techniques in the mid-infrared (MIR) frequency range, relevant for
spectroscopic studies of molecular fingerprints [29], semiconductors [30] and
various light-matter interactions [31]. Its success in the classical regime
has motivated its promotion to the quantum domain, where a team lead by
Leitenstorfer achieved the first milestone towards quantum sensing: the direct
detection of the MIR quantum vacuum in the time-domain with sub-cycle
resolution [32]. Since then, direct detection of terahertz quantum vacuum as
well as its two-point correlations have been reported [33] based on the EOS
variants. However, there remain several milestones this technique must reach
in order to establish itself as a reliable MIR quantum sensing technique in
the time domain. Prominently among those are the identification of the quasi-
time-domain quadratures and the ability to detect both field quadratures [23]
and their correlations, and the experimental demonstration of extracting of
quantum information from the field.

Our Letter develops the theoretical framework for the former by accurately
identifying the quadrature fields in the time domain and recasting the
measurement protocols of EOS to be mode-matched to the detection of both
simultaneously, hence achieving the direct measurement of the Husumi-Q
function. Gained insights of mode-matching are exemplified through a predicted
improvement upon a recent experimental proposal for the detection of a quantum
squeezed state[34], which addresses the latter.

The remainder of this Letter is organized as follows: We first identify the
quasi-time-domain quadratures as frequency-filtered time-domain electric-field
and its Hilbert transform in Sec. II. We then optimize the mode-matching and
coupling efficiency in Sec. III, to mitigate the introduction of undesirable
vacuum noise and increase sensitivity to squeezing. In Sec. IV we promote EOS
to efficiently couple to both the frequency-filtered electric field and its
Hilbert transform simultaneously, giving the statistics of the Husimi-Q
function which can utilized for full quantum-state reconstruction. We
numerically demonstrate its effectiveness in quantum-state reconstruction of a
squeezed state.

##  II Electric Field and its conjugate in the time-domain

Fields in (3+1) dimensions can be simplified to travelling-wave modes in (1+1)
dimensions in the paraxial approximation, where they are considered to
propagate along a single direction [35, 19, 36, 23]. In such scenario we can
decompose the electric field operator, $\hat{E}$, as

$\begin{gathered}\hat{E}_{\sigma}(t,x)=\int_{-\infty}^{\infty}\mathrm{d}\omega\,E_{\omega}(t,x)(\,\hat{a}_{\omega,\sigma})\,\text{,
with}\\\
E_{\omega}(t,x)=-i\,\text{sign}(\omega)\sqrt{\frac{\hbar\lvert\omega\rvert
c}{4\pi n_{\omega}A}}e^{-i\omega(t-\frac{n_{\omega}x}{c})}\,,\end{gathered}$
(1)

where $A$ accounts for the transversal extension of the field, $n_{\omega}$ is
the frequency-dependent index of refraction of the propagation medium, and the
ladder operators $\hat{a}_{\omega,I}$ satisfy the canonical commutation
relations
$[\hat{a}_{\omega,\sigma},\hat{a}_{\omega^{\prime},\sigma^{\prime}}]=\text{sign}(\omega)\delta(\omega+\omega^{\prime})\delta^{\sigma}_{\sigma^{\prime}}$.
Here, we have adopted the notation where $\sigma$ refers to the polarisation
degree of freedom, and
$\hat{a}_{\omega,\sigma}=\hat{a}_{-\omega,\sigma}^{{\dagger}}$. The Hilbert
quadrature ($H(\hat{E})(t)=\hat{H}(t)$) of the electric field is then
decomposed as

$\begin{gathered}{H}(\hat{E}_{\sigma})(t,x)=\int_{-\infty}^{\infty}\mathrm{d}\omega\,i\,\text{sign}(\omega)E_{\omega}(t,x)\,(\hat{a}_{\omega,\sigma})\,.\end{gathered}$
(2)

Figure 1: The waveform of bandwidth-limited electric/Hilbert field are plotted
in red/blue. These can be compared to the standard electric and Hilbert field
which is a $\delta$-function and $1/t$ respectively. The red/blue dashed lines
represent the coupled-mode via EOS.

The fundamental reason for the impossibility of measuring either quadratures
strictly in the time-domain is the requirement for infinite detection
bandwidth. This mirrors the fundamental reason for the impossibility of
measuring both quadratures strictly in the frequency-domain, namely the
requirement for infinite detection time. Any detection process will feature
both finite bandwidth and a finite amount of interaction time with the field
under probe.

One may naively believe that approaching infinite bandwidth/time is desirable
for the direct detection of quadratures in time/frequency-domain. However, in
a real set-up this leads to a coupling with fewer photons (e.g. there are
fewer photons in a shorter time-window), reducing the coupling efficiency to
the signal [34]: the preciseness in frequency/time comes at the cost of
detection efficiency. The detection efficiency approaches zero for the direct
detection of quadrature in the time/frequency-domain.

In the frequency-domain, the conventional approach is to impose a time-window
based off information of the signal that is known _a priori_. As the signal
does not exist outside this time-window, we do not lose any information by
imposing this time-window to the probe. In the same light, we can impose an
appropriate bandwidth when extracting the information of the light signal. The
bandwidth-limited electric field and its Hilbert transform are defined as:

$\begin{gathered}\hat{E}_{BL,I}(t,x)=\int_{-\omega_{m}}^{\omega_{m}}\mathrm{d}\omega\,E_{\omega}(t,x)\hat{a}_{\omega,I}\,,\text{,
with}\\\
\hat{H}_{BL,I}(t,x)=i\int_{-\omega_{m}}^{\omega_{m}}\mathrm{d}\omega\,\text{sign}(\omega)E_{\omega}(t,x)\hat{a}_{\omega,I}\,.\end{gathered}$
(3)

This simplification renders both $\hat{E}(t)$ and $\hat{H}(t)$ to be
detectable (refer to Fig. 1), without losing any relevant information of the
signal under study. Furthermore, the finite bandwidth introduces the notion of
mode-matching, which allows us to quantify how accurately we couple to these
quadrature modes.

##  III Mode-matching and detection efficiency of EOS

Electro-optic sampling (EOS) is one of the most promising candidates in ultra-
fast photonics for the detection of MIR electric field in the time-domain [24,
26, 27, 28]. It utilizes a $\chi^{(2)}$ based interaction between near-
infrared (NIR) and mid-infrared (MIR) via sum- and difference- frequency
generation. By driving this interaction with a coherent probe-pulse with an
envelope smaller than a single wavelength of the signal (i.e. sub-cycle with
respect to the MIR signal), the sub-cycle quantum information is up-converted
to the NIR. [25]. This information is then analyzed via the ellipsometry
scheme. For a long time, experimental efforts have been made towards the
direct detection of MIR electric field in the time-domain by utilizing an
ultra-short pulse, reaching sub-cycle resolution under 6 fs [32]. This
experimental effort for shorter pulses was motivated by the idea that the
ideal probe is a Dirac-delta distribution. The introduction of $\hat{E}_{BL}$
allows us to take a different narrative approach.

Figure 2: We demonstrate the effect of bandwidth on coupling intensity and
mode-matching efficiency. The black (dashed) line correspond to the constant
photon number (bandwidth-limited) coupling intensity plots. The grey (dashed)
line correspond to the constant intensity (bandwidth-limited), with coupling
intensity plots. The red line correspond to the bandwidth-limited mode-
matching efficiency as a function of bandwidth.

Let us consider a NIR signal with a central frequency of $\Omega_{0}=20$ THz.
The information of such pulse is generally limited to be under $2\Omega_{0}$
[37, 22]. When this assumption is satisfied, we can neglect the quantum
information of frequencies above $2\Omega_{0}$ when analyzing this signal, and
as a result we can utilize the bandwidth-limited electric field (refer to Eq.
3) with $\omega_{m}=2\Omega_{0}$. Our goal is now tailored towards the
detection of $\hat{E}_{BL}$. We can detect $\hat{E}_{BL}$ extremely
accurately, with a mode-matching efficiency of $99\%$ by utilizing the EOS
scheme with a sinc probe pulse. In Fig. 2 we investigate the optimal bandwidth
of the probe pulse to detect $\hat{E}_{BL}$ by analysing the effect of probe-
pulse bandwidth on bandwidth-limited coupling intensity, $\theta_{E}$, and
mode-matching efficiency, $\gamma_{E}$, which are defined as:

$\displaystyle\gamma_{E}=:[\hat{E}_{BL},Tr_{BL}(\hat{\textbf{E}})]:/\sqrt{\theta_{e}\theta_{E}}\,,$
(4)
$\displaystyle\theta_{E}=:[Tr_{\omega_{m}}({\hat{\textbf{E}}}),Tr_{BL}({\hat{\textbf{E}})}]:\,,$
(5) $\displaystyle\theta_{e}=:[\hat{E}_{BL},{\hat{E}_{BL}}]:\,.$ (6)

Here $\hat{\textbf{E}}$ is the detected electric-field via EOS (refer to Eq.
(10)),
$Tr_{BL}[\hat{a}_{\omega}]=(1-\Pi((\omega)/(2\omega_{m})))\hat{a}_{\omega}$,
$\Pi(x)=1\,\forall-0.5\leq x\leq 0.5$, otherwise $\Pi(x)=0$, and we have
defined $:[,]:$ as the normal-ordered commutator, with
$:[\hat{a}_{\pm|\omega|},\hat{a}_{\mp|\omega^{\prime}|}]:=[\hat{a}_{|\omega|},\hat{a}_{|\omega^{\prime}|}^{{\dagger}}]$.

We find in Fig. 2, that when we keep the probe photon-number constant (at
$5\times 10^{9}$) the bandwidth-limited coupling intensity plateaus at about
$4\%$ for large bandwidth. In comparison, the bandwidth-unlimited coupling
intensity keeps increasing, implying that a larger bandwidth THz introduces
additional noise from undesired frequencies. The amount of unwanted noise can
be explicitly taken as the difference between to the bandwidth -limited and
-unlimited detection efficiency. This implies we are using our resources in an
inefficient manner, and we follow our analysis by considering the same plot,
but with a constant intensity of the probe pulse motivated by the damage
threshold of EOX crystal. In this scenario, a smaller bandwidth (roughly $70$
THz) probe pulse shows a significant increase in coupling intensity, reaching
roughly $9\%$ in coupling intensity. It is noted that increasing the photon
number does not change the ratio between the wanted and unwanted noise.

We also plot the bandwidth-limited mode-matching efficiency, which shows that
a mode-matching of $99\%$ is achieved when the full-bandwidth is utilized. We
constrain ourselves to a mode-matching efficiency that is at least of that
achieved by the full-bandwidth. Under this condition, the optimal coupling
intensity of $8\%$ is achieved at around $120$ THz. By considering the
appropriate spectral support we suppress the admixture of unwanted vacuum
fluctuation and by efficiently mode-matching our probe specifically to the
signal in question, we improved the detection efficiency [34] by roughly two
folds without the cost of accuracy.

##  IV Full Tomography

(a) Spectral decomposition of probe-pulse electric-field

(b) Experimental set-up

Figure 3: a) We show the probe pulse electric-field that is utilised to drive
the EOX interaction. The signal below/above $275$ THz is dedicated towards the
detection of the electric-field/Hilbert-field quadrature. b) We present a
novel schematic set-up for the simultaneous detection of the electric- and
Hilbert- quadrature via EOS. The probe pulse goes into the EOX crystal in the
Z-polarisation , while the signal in the S-polarisation . This is followed by
a spectral filter beam-splitter , splitting the three frequency regimes (red
for below $275$ THz, purple for above $340$ THz and blue for frequencies in-
between). The first port has a half-waveplate in the 22.5 degrees from the
z-polarisation. The second and third port has a quarter-waveplate 45 degrees
from the z-polarisation. This is followed by a Wollaston prism which separates
the S- and Z- polarisation fields, where the intensity difference between the
polarisations are recorded.

One of EOS’s strength over homodyne detection is its ability to isolate the
electric-field from the Hilbert quadrature regardless of the probe pulse’s
phase. This introduces the notion of an absolute phase [23, 38] of the signal
and we lose the meaning of relative phase.

In contrast, for homodyne detection we are agnostic to whether we are coupling
directly to the electric-/Hilbert- field, but we obtain information on the
relative phase between the signal and probe. In homodyne detection, the
relative phase is altered by a time-delay to the signal. For EOS, this method
cannot be used to alter the absolute phase, as it simply shifts the time that
is sampled without affecting the quadrature (e.g.
$\hat{E}_{BL}(t)\rightarrow\hat{E}_{BL}(t^{\prime})$). In the following sub-
section we identify how EOS can be altered to detect the Hilbert transform.

###  IV.1 Detection of the Hilbert transform

When we are analysing the electric field, we can obtain its localized property
in the time-domain by a mode with an even parity. Its Hilbert transform is
represented by a mode with odd parity which is a delocalized property of the
electric-field (refer to Fig. 1). From the motivational stand point of view,
EOS aimed at extracting the local property of the electric-field and therefore
fails to extract the Hilbert transform. There are two fundamental pillars of
EOS that needs to be removed for the detection of the Hilbert transform: its
ability to detect phase-locked (in frequency) electric-field mode and
localized modes. The prior can be resolved by the introduction of spectral
filters, and adding a different phase contribution to different frequencies.
Changing from a flat phase with respect to the electric field to an odd-phase
contribution (positive above the central frequency, negative below the central
frequency) is desirable for the detection of a odd-function (e.g. the Hilbert
transform). As the Hilbert quadrature is a delocalized quantity of the
electric-field, we benefit from utilizing a delocalized probe-pulse with a
strong central frequency contribution (refer to Fig. 3(a) for the exact
waveform). In the App. A, we analyse the effect of bandwidth on detection
efficiency and coupling strength, and find that a peak mode-matching of
$99.9\%$ is achieved at $80$ THz bandwidth with a coupling intensity of
roughly $2\%$.

###  IV.2 Simultaneous detection of both field quadratures

So far, we have established a mean to detect the field and its Hilbert
transform independently. We can conduct both simultaneously via multiplexing
and assigning different frequency regime to the detection of either. We
demonstrate this experimental setup in Fig. 3(b). In this paper we evaluate
the effect of the optical elements via the Heisenberg picture. Driving the EOX
crystal utilizing a strong coherent signal $\alpha_{z}(t,x)$ leads to a
Hamiltonian of the following form [36]:

$\hat{H}_{\chi}(t)=\int_{-\infty}^{\infty}\mathrm{d}x\lambda{\alpha}_{z}(t,x)\hat{E}_{S}(t,x)\hat{E}_{S}(t,x)\Pi\left(\frac{x}{L}\right),$
(7)

where $\lambda=\frac{A\epsilon_{0}d}{2}$, with the coupling constant
$d=-n^{4}r_{41}$, expressed in terms of the electro-optic susceptibility
coefficient, $r_{41}$, and the refractive index of the crystal, $n$ [39, 36].
We consider, without the loss of generality, a zincblende-type crystal with a
length of $L=7\,\mathrm{\mu m}$. Its electro-optic coefficient is taken as
$r_{41}=4~{}\mathrm{pm/V}$ (for the particular case of ZnTe). The refractive
index, $n_{\Omega}$, varies only slightly (from 2.55 to 2.59) in the MIR [39].
We utilize a fit for the refractive index in the NIR frequency range
$n_{{\omega}}$ [40, 36]. We compute the effect of this Hamiltonian by
utilizing a unitary correction [36] to second order perturbation [41, 42].
This is followed by a spectral filter, splitting the output into three
frequency components:

$\begin{gathered}\hat{N}_{i,\sigma}=\int^{\omega_{i,max}}_{\omega_{i,min}}\mathrm{d}\omega\;\hat{a}_{\omega,\sigma}^{{\dagger}}\hat{a}_{\omega,\sigma}\,,\end{gathered}$
(8)

Each port goes through either a half- or quarter- waveplate:

$\hat{U}(\theta_{i},\phi_{i})=\hat{e}^{i\phi_{i}(\cos(\theta_{i})\hat{N}_{i,S}+\sin(\theta_{i})(\hat{N}_{i,Z})}\,.$
(9)

Details of $\theta_{i}$ and $\phi_{i}$ can be found in Fig. 3(b). In each
branch, the waveplate is followed by a Wollaston prism and a balanced pair of
photo-detectors. The E-field is given by

$\hat{\textbf{E}}(t)=\frac{\hat{N}_{z,1}^{\prime}-\hat{N}_{s,1}^{\prime}}{\sqrt{\braket{\hat{N}_{z,1}^{\prime}+\hat{N}_{s,1}^{\prime}}}}\,,$
(10)

and its Hilbert transform,

$\hat{\textbf{H}}(t)=\frac{\hat{N}_{z,2}^{\prime}-\hat{N}_{s,2}^{\prime}}{\sqrt{\braket{\hat{N}_{z,2}^{\prime}+\hat{N}_{s,2}^{\prime}}}}+\frac{\hat{N}_{s,3}^{\prime}-\hat{N}_{z,3}^{\prime}}{\sqrt{\braket{\hat{N}_{z,3}^{\prime}+\hat{N}_{s,3}^{\prime}}}}\,.$
(11)

Simultaneous detection of both allows the detection of the Husumi-Q function.
In App. B, we discuss two other promising schemes that can be implemented for
signals with larger bandwidth. One of these scheme is analogous to the method
discussed in this section, utilizing a beam-splitter instead of a spectral
filter (refer to App. B.1) , while the other can be utilized for full
tomography with an arbitrary phase (refer to App. B.2), which can be utilized
to extract the Wigner function.

##  V Numerical Results

The sensitivity of the $n$th moment utilizing the standard EOS setup scales as
$\epsilon^{n}$, where $\epsilon=\gamma_{E}$ is the coupling efficiency. This
makes it extremely challenging for EOS to be implemented for quantum sensing
beyond the second moment. However, a recent novel approaches to EOS utilizing
a quantum probe has shown an enhancement to sensitivity to the higher moments
[43]. For the purposes of extracting the information of a Gaussian state, it
is sufficient to compute the second moment [44, 45], extracted through:

$\begin{gathered}\Delta
V_{E}=(\braket{\hat{\textbf{E}}(t)^{2}}-\braket{\hat{\textbf{E}}(t)}^{2})-(\braket{\hat{\textbf{E}}_{0}(t)^{2}}-\braket{\hat{\textbf{E}}_{0}(t)}^{2})\,,\\\
\Delta
V_{H}=(\braket{\hat{\textbf{H}}(t)^{2}}-\braket{\hat{\textbf{H}}(t)}^{2})-(\braket{\hat{\textbf{H}}_{0}(t)^{2}}-\braket{\hat{\textbf{H}}_{0}(t)}^{2})\,,\end{gathered}$
(12)

where we have defined $\hat{\textbf{E}}_{0}$/ $\hat{\textbf{H}}_{0}$ as their
respective operators $\hat{\textbf{E}}$/ $\hat{\textbf{H}}$ in the absence of
the signal (i.e. in vacuum). We utilize this set-up to analyze the statistics
of a squeezed Gaussian signal:

$\begin{gathered}\hat{a}_{G}=\int_{0}^{\infty}\mathrm{d}\Omega\;G(\Omega)\hat{a}_{\Omega,Z}\\\
G(\Omega)=B\sqrt{\omega}(\frac{1}{2\pi\sigma_{G}^{2}})^{1/4}e^{(\Omega_{0}-\Omega)^{2}/(4\sigma_{G}^{2})}\end{gathered}$
(13)

with $\Omega_{0}=20$ THz, $\sigma_{G}=4$ THz and a squeezing strength of
$r=0.5$. Our approach utilizing simultaneous measurement of both
$\hat{E}_{BL}$ and $\hat{B}_{BL}$ allows the reconstruction of the Husimi
Q-function (with added noise) in the time-domain, showing the correlation
between the quadratures. The detected second moment and their correlation are
plotted in the Fig. 4, and compared to the re-scaled second moment of the
field. The results demonstrate EOS’s ability to conduct full tomography of a
squeezed state at a relatively high accuracy. The difference between the two
results is associated with the imperfect mode-matching, which can be traced
back to the use of the same probe-pulse for both the $\chi^{(2)}$ interaction
and ellipsometry, and other technical effects such as phase-matching.

Figure 4: The detected second moment (solid lines) are compared to the re-
scaled second moment of the field (dashed lines). The red correspond to the
electric-field quadrature, the blue correspond to the Hilbert quadrature and
the purple correspond to their positive correlation.

##  VI Conclusion

Our work establishes a method to conduct full tomography of a MIR squeezed
state in the time-domain via EOS and introduces a quantitative method to
analyse the quality of mode-matching, all of which are important milestones
towards promoting EOS to a time-domain quantum sensing technique. Although the
direct detection of $\hat{E}(t)$ and $\hat{H}(t)$ cannot be conducted in
theory [19, 20, 21], the quantum information of interest will exist within a
certain bandwidth. By introducing the same bandwidth to the electric field and
its Hilbert transform quadrature, they are rendered detectable:
$\hat{E}_{BL}(t)$ and $\hat{H}_{BL}(t)$. Contrary to the standard approach to
EOS which favors the use of a shorter probe pulse [34], we have demonstrated
that we can improve detection efficiency with the use of a longer probe-pulse,
without the cost of accuracy.

The standard EOS set-up was shown to be very inefficient in coupling with the
$\hat{H}_{BL}(t)$ field, and undetectable without an effective spectral filter
[46]. We have created an efficient and accurate scheme to detect the Hilbert
transform by introducing a phase-jump in frequency for the detection, and
utilizing a delocalized probe field. Furthermore, we achieved simultaneous
detection of both $\hat{E}_{BL}(t)$ and $\hat{H}_{BL}(t)$ via multiplexing. An
experimental realization of this set-up would be the first example of full-
tomography in the time-domain, and would be one of the most crucial advances
in quantum sensing for MIR and ultra-fast photonics.

We found that a perfect phase-matching to $\hat{E}_{BL}(t)$ and
$\hat{H}_{BL}(t)$ was not possible when the same probe pulse was utilized for
driving the $\chi^{(2)}$ interaction and ellipsometry, especially in the
presence of phase-matching. For future research, implementing the use of
template search [47] and a different probe-pulse for the $\chi^{(2)}$
interaction and ellipsometry as viable options to achieve perfect replication
of the Husumi-Q function.

## References

  * Glauber [1963] R. J. Glauber, The quantum theory of optical coherence, Phys. Rev. 130, 2529 (1963).
  * Kelley and Kleiner [1964] P. L. Kelley and W. H. Kleiner, Theory of electromagnetic field measurement and photoelectron counting, Phys. Rev. 136, A316 (1964).
  * Yuen and Shapiro [1978] H. P. Yuen and J. H. Shapiro, Quantum statistics of homodyne and heterodyne detection, in  _Coherence and Quantum Optics IV_ , edited by L. Mandel and E. Wolf (Springer US, Boston, MA, 1978) pp. 719–727.
  * Collett  _et al._ [1987] M. Collett, R. Loudon, and C. Gardiner, Quantum theory of optical homodyne and heterodyne detection, Journal of Modern Optics 34, 881 (1987).
  * Vogel and Risken [1989] K. Vogel and H. Risken, Determination of quasiprobability distributions in terms of probability distributions for the rotated quadrature phase, Phys. Rev. A 40, 2847 (1989).
  * Welsch  _et al._ [1999] D.-G. Welsch, W. Vogel, and T. OpatrnÃœ, Ii homodyne detection and quantum-state reconstruction (Elsevier, 1999) pp. 63–211.
  * Yuen and Chan [1983] H. P. Yuen and V. W. S. Chan, Noise in homodyne and heterodyne detection, Opt. Lett. 8, 177 (1983).
  * Abbas  _et al._ [1983] G. L. Abbas, V. W. S. Chan, and T. K. Yee, Local-oscillator excess-noise suppression for homodyne and heterodyne detection, Opt. Lett. 8, 419 (1983).
  * Schumaker [1984] B. L. Schumaker, Noise in homodyne detection, Opt. Lett. 9, 189 (1984).
  * Yurke [1985a] B. Yurke, Wideband photon counting and homodyne detection, Phys. Rev. A 32, 311 (1985a).
  * Braunstein [1990] S. L. Braunstein, Homodyne statistics, Phys. Rev. A 42, 474 (1990).
  * Yurke [1985b] B. Yurke, Squeezed-coherent-state generation via four-wave mixers and detection via homodyne detectors, Phys. Rev. A 32, 300 (1985b).
  * Aytür and Kumar [1992] O. Aytür and P. Kumar, Squeezed-light generation with a mode-locked q-switched laser and detection by using a matched local oscillator, Opt. Lett. 17, 529 (1992).
  * Shapiro and Shakeel [1997] J. H. Shapiro and A. Shakeel, Optimizing homodyne detection of quadrature-noise squeezing by local-oscillator selection, J. Opt. Soc. Am. B 14, 232 (1997).
  * Loudon and Knight [1987] R. Loudon and P. Knight, Squeezed light, Journal of Modern Optics 34, 709 (1987).
  * Slusher  _et al._ [1985] R. E. Slusher, L. W. Hollberg, B. Yurke, J. C. Mertz, and J. F. Valley, Observation of squeezed states generated by four-wave mixing in an optical cavity, Phys. Rev. Lett. 55, 2409 (1985).
  * Wu  _et al._ [1986] L.-A. Wu, H. J. Kimble, J. L. Hall, and H. Wu, Generation of squeezed states by parametric down conversion, Phys. Rev. Lett. 57, 2520 (1986).
  * Smithey  _et al._ [1993] D. T. Smithey, M. Beck, M. G. Raymer, and A. Faridani, Measurement of the wigner distribution and the density matrix of a light mode using optical homodyne tomography: Application to squeezed states and the vacuum, Phys. Rev. Lett. 70, 1244 (1993).
  * Scully and Zubairy [1997] M. O. Scully and M. S. Zubairy,  _Quantum Optics_ (Cambridge University Press, 1997).
  * [20] .
  * Milburn and Walls [2008] G. J. Milburn and D. Walls,  _Quantum Optics_ (Springer, 2008).
  * Gulla and Skaar [2021] J. Gulla and J. Skaar, Approaching single-photon pulses, Phys. Rev. Lett. 126, 073601 (2021).
  * Virally and Reulet [2019] S. Virally and B. Reulet, Unidimensional time-domain quantum optics, Phys. Rev. A 100, 023833 (2019).
  * Leitenstorfer  _et al._ [1999] A. Leitenstorfer, S. Hunsche, J. Shah, M. C. Nuss, and W. H. Knox, Detectors and sources for ultrabroadband electro-optic sampling: Experiment and theory, Applied Physics Letters 74, 1516 (1999).
  * Gallot and Grischkowsky [1999] G. Gallot and D. Grischkowsky, Electro-optic detection of terahertz radiation, J. Opt. Soc. Am. B 16, 1204 (1999).
  * Huber  _et al._ [2000] R. Huber, A. Brodschelm, F. Tauser, and A. Leitenstorfer, Generation and field-resolved detection of femtosecond electromagnetic pulses tunable up to 41 thz, Applied Physics Letters 76, 3191 (2000).
  * KÃŒbler  _et al._ [2004] C. KÃŒbler, R. Huber, S. TÃŒbel, and A. Leitenstorfer, Ultrabroadband detection of multi-terahertz field transients with gase electro-optic sensors: Approaching the near infrared, Applied Physics Letters 85, 3360 (2004).
  * Sell  _et al._ [2008] A. Sell, R. Scheu, A. Leitenstorfer, and R. Huber, Field-resolved detection of phase-locked infrared transients from a compact er:fiber system tunable between 55 and 107 thz, Applied Physics Letters 93, 251107 (2008).
  * Schweinberger  _et al._ [2016] W. Schweinberger, L. Vamos, M. Huber, J. Zhang, M. Trubetskov, N. Lilienfein, A. Apolonski, V. Pervak, N. Karpowicz, O. Pronin, and I. Pupeza, Broadband mid-infrared time-domain spectrometer for the molecular fingerprint region, in  _High-Brightness Sources and Light-Driven Interactions_ (Optica Publishing Group, 2016) p. MT1C.5.
  * Kira and Koch [2006] M. Kira and S. W. Koch, Quantum-optical spectroscopy of semiconductors, Phys. Rev. A 73, 013813 (2006).
  * Dorfman  _et al._ [2016] K. E. Dorfman, F. Schlawin, and S. Mukamel, Nonlinear optical signals and spectroscopy with quantum light, Rev. Mod. Phys. 88, 045008 (2016).
  * Riek  _et al._ [2015] C. Riek, D. Seletskiy, A. S. Moskalenko, J. Schmidt, P. Krauspe, S. Eckart, S. Eggert, G. Burkard, and A. Leitenstorfer, Direct sampling of electric-field vacuum fluctuations, Science 350, 420 (2015).
  * Settembrini  _et al._ [2022] F. F. Settembrini, F. Lindel, A. M. Herter, S. Y. Buhman, and J. Fais, Detection of quantum-vacuum field correlations outside the light cone, Nat Commun 13, 3383 (2022).
  * Kizmann  _et al._ [2019] M. Kizmann, T. L. d. M. Guedes, D. V. Seletskiy, A. S. Moskalenko, A. Leitenstorfer, and G. Burkard, Subcycle squeezing of light from a time flow perspective, Nature Physics 15, 960 (2019).
  * Blow  _et al._ [1990] K. J. Blow, R. Loudon, S. J. D. Phoenix, and T. J. Shepherd, Continuum fields in quantum optics, Physical Review A 42, 4102 (1990).
  * Onoe  _et al._ [2022] S. Onoe, T. L. M. Guedes, A. S. Moskalenko, A. Leitenstorfer, G. Burkard, and T. C. Ralph, Realizing a rapidly switched unruh-dewitt detector through electro-optic sampling of the electromagnetic vacuum, Phys. Rev. D 105, 056023 (2022).
  * Bialynicki-Birula [1998] I. Bialynicki-Birula, Exponential localization of photons, Phys. Rev. Lett. 80, 5247 (1998).
  * Onoe and Ralph [2019] S. Onoe and T. C. Ralph, Universal transformation of displacement operators and its application to homodyne tomography in differing relativistic reference frames, Phys. Rev. D 99, 116001 (2019).
  * Moskalenko  _et al._ [2015] A. S. Moskalenko, C. Riek, D. V. Seletskiy, G. Burkard, and A. Leitenstorfer, Paraxial theory of direct electro-optic sampling of the quantum vacuum, Physical Review Letters 115, 263601 (2015).
  * Marple [1964] D. T. F. Marple, Refractive index of znse, znte, and cdte, Journal of Applied Physics 35, 539 (1964).
  * Guedes  _et al._ [2022] T. L. M. Guedes, I. Vakulchyk, D. V. Seletskiy, A. Leitenstorfer, A. S. Moskalenko, and G. Burkard, Back action in quantum electro-optic sampling of electromagnetic vacuum fluctuations (2022).
  * Boyd [2019] R. W. Boyd,  _Nonlinear optics_ (Academic press, 2019).
  * Virally  _et al._ [2021] S. Virally, P. Cusson, and D. V. Seletskiy, Enhanced electro-optic sampling with quantum probes, Phys. Rev. Lett. 127, 270504 (2021).
  * Weedbrook  _et al._ [2012] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Gaussian quantum information, Reviews of Modern Physics 84, 621 (2012).
  * Adesso  _et al._ [2014] G. Adesso, S. Ragy, and A. R. Lee, Continuous variable quantum information: Gaussian states and beyond, Open Systems & Information Dynamics 21, 1440001 (2014).
  * Sulzer  _et al._ [2020] P. Sulzer, K. Oguchi, J. Huster, M. Kizmann, T. L. M. Guedes, A. Liehl, C. Beckh, A. S. Moskalenko, G. Burkard, D. V. Seletskiy, and A. Leitenstorfer, Determination of the electric field and its hilbert transform in femtosecond electro-optic sampling, Phys. Rev. A 101, 033821 (2020).
  * Chatwin-Davies  _et al._ [2022] A. Chatwin-Davies, A. Kempf, and P. Simidzija, Predictions for quantum gravitational signatures from inflation (2022).

##  Appendix A Hilbert transform

In this section we present the effect of bandwidth on the coupling intensity,
$\theta_{H}$, and mode-matching efficiency, $\gamma_{H}$, defined as:.

Figure 5: We demonstrate the effect of bandwidth on coupling intensity and
mode-matching efficiency. The black (dashed) line correspond to the constant
photon number (bandwidth-limited) coupling intensity plots. The grey (dashed)
line correspond to the constant intensity (bandwidth-limited) coupling
intensity plots. The red line correspond to the bandwidth-limited mode-
matching efficiency as a function of bandwidth.

$\displaystyle\theta_{H}=Tr_{\omega_{m}}(:[\hat{\mathbf{H}},\hat{\mathbf{H}}:)\,,$
(14)
$\displaystyle\gamma_{H}=Tr_{\omega_{m}}(:[\hat{a}_{H,S}^{\prime},\hat{\mathbf{H}}]:)/\sqrt{\theta\theta_{H}}\,.$
(15)

Fig. 5 shows a similar trend to that of Fig. 2. For constant probe photon
number, the bandwidth limited coupling intensity plateaus from a bandwidth of
$80$ ThZ, giving a coupling intensity of $0.005$. The bandwidth-unlimited
coupling intensity increases to $0.016$ when the full bandwidth is utilized,
implying the extra bandwidth introduces unwanted noise. In comparison, for
constant probe intensity, the bandwidth limited coupling intensity optimizes
at around $80$ THz, with a coupling intensity of roughly $0.015$. This is a
three fold increase in bandwidth-limited coupling intensity compared to when
the full bandwidth was utilized. The mode-matching efficiency optimizes at
about $99.9\%$ at around $80$ THz, then decreases for larger bandwidth due to
phase-matching considerations.

##  Appendix B Experimental set-ups

The experimental method discussed in the main text allocated different
bandwidths for the detection of the electric field and its Hilbert transform.
Splitting the bandwidth does not come with any significant costs when
considering a signal with $\Omega_{0}=20$THz, however for signals with large
bandwidth, there is not enough bandwidth of the probe-pulse to allow multi-
plexing. In this section, we discuss two alternate methods that can be
implemented for detection of fields with larger bandwidths (signal with a
central frequency of $25$ THz and above).

###  B.1 Beam-splitter method

This method works by utilizing a probe pulse that is sufficiently efficient
and accurate for both the detection of both the electric field and the Hilbert
transform. In this paper, we suggest the use of a coherent superposition of
the probe pulse to detect the electric field signal and the Hilbert transform
signal. In general, the coupling to the Hilbert transform is stronger than the
coupling to the electric-field. To account for this effect, we utilize a beam-
splitter with a suitable transmission coefficient so that more photons are
transmitted to the port extracting the Hilbert transform than the port
extracting the electric field (e.g. if the coupling intensity for the electric
field is $\theta_{E}=0.1$ and its Hilbert transform is $\theta_{H}=0.01$, then
we utilize a beam-splitter which sends 10 times more photons in the Hilbert
transform port). The electric field and Hilbert transform are simultaneously
analysed in each port. The experimental set-up is shown in Fig. 6.

Figure 6: Experimental set-up for simultaneous detection of both quadratures
for signals with a larger bandwidth.

###  B.2 Tomography method

There are several ways to conduct full tomography of the field. One is to
extract the electric field and its Hilbert transform simultaneously. Another
is to conduct an effective homodyne detection with various phases: an
independent measurement of
$\hat{X}(\phi)=\hat{a}e^{-i\phi}+\hat{a}^{{\dagger}}e^{i\phi}$ with
$\phi=\\{0,\pi/2,\pi,3\pi/2\\}$ is equivalent to a simultaneous measurement of
$\hat{X}(\phi)$ with $\phi=\\{0,\pi/2\\}$. The prior can be achieved by
utilizing a spectral filter to split the probe frequency above/below the
central frequency into two ports. The outputs are followed by both quarter-
and half-waveplate before the Wollaston prism. Careful tuning of both allows
the detection of $\hat{X}(\phi)$ with an arbitrary phase.

Figure 7: Experimental set-up for sub-cycle weak homodyne detection, where the
phase can be tuned via the angle of the half- and quarter-waveplates.


# Uncertainty in latent representations of variational autoencoders optimized
for visual tasks

Josefina Catoni
Research Institute for Signals, Systems
and Computational Intelligence sinc(i)
FICH-UNL/CONICET
Santa Fe, Argentina
<EMAIL_ADDRESS>
&Enzo Ferrante
Research Institute for Signals, Systems
and Computational Intelligence sinc(i)
FICH-UNL/CONICET
Santa Fe, Argentina
&Diego H. Milone
Research Institute for Signals, Systems
and Computational Intelligence sinc(i)
FICH-UNL/CONICET
Santa Fe, Argentina
&Rodrigo Evcheveste
Research Institute for Signals, Systems
and Computational Intelligence sinc(i)
FICH-UNL/CONICET
Santa Fe, Argentina

###### Abstract

Deep learning methods are increasingly becoming instrumental as modeling tools
in computational neuroscience, employing optimality principles to build
bridges between neural responses and perception or behavior. Developing models
that adequately represent uncertainty is however challenging for deep learning
methods, which often suffer from calibration problems. This constitutes a
difficulty in particular when modeling cortical circuits in terms of Bayesian
inference, beyond single point estimates such as the posterior mean or the
maximum a posteriori. In this work we systematically studied uncertainty
representations in latent representations of variational auto-encoders (VAEs),
both in a perceptual task from natural images and in two other canonical tasks
of computer vision, finding a poor alignment between uncertainty and
informativeness or ambiguities in the images. We next showed how a novel
approach which we call explaining-away variational auto-encoders (EA-VAEs),
fixes these issues, producing meaningful reports of uncertainty in a variety
of scenarios, including interpolation, image corruption, and even out-of-
distribution detection. We show EA-VAEs may prove useful both as models of
perception in computational neuroscience and as inference tools in computer
vision.

_K_ eywords Uncertainty $\cdot$ Perception $\cdot$ Variational Inference
$\cdot$ NeuroAI $\cdot$ VAEs

## Introduction

The interplay between neuroscience and artificial intelligence (AI) has been
steadily growing and has proven mutually beneficial for both fields. Indeed, a
large number of experts posit that research in ‘NeuroAI’ is key to help
develop the next generation of artificial intelligence methods [1]. In turn,
the use of modeling tools from machine learning (ML), has proven highly
effective in shedding light on how the brain acquires, processes and stores
information, [2, 3, 4, 5, 6, 7, 8]. In this vein, a fruitful line of modeling
work has been that of ‘goal-driven’ approaches, which are framed in terms of a
computational task that the system needs to solve. Methods employing
artificial neural networks (ANNs) have been able to successfully capture
features of cortical representations in multiple areas of the brain such as
the visual [2] and auditory systems [3]. The main hypothesis is that
optimizing an ANN for a task that is relevant for an area of the brain may
reveal properties of neural coding and perception in that area [9]. ML hence
brings powerful novel computational tools to a longstanding tradition of
explaining neural properties in terms of optimal information processing [10].

Perception involves making inference about non-observable variables of
interest ($\mathbf{y}$) from information coming through our senses
($\mathbf{x}$), which is always partial and incomplete. Bayes rule provides an
optimal way to perform computations in probabilistic scenarios [11], and hence
it is a natural tool to understand perception in the context of uncertainty.
Indeed, Bayesian theories of perception postulate that the brain performs
probabilistic inference to estimate posterior probability distributions
$\mathcal{P}(\mathbf{y}|\mathbf{x})$ for those non-observable (also called
latent) variables $\mathbf{y}$[12]. This computation allows to ‘fill-in the
gaps’ left by the uncertainty of the observation –captured by the _likelihood
function_ $\mathcal{P}(\mathbf{x}|\mathbf{y})$– with prior knowledge about
these variables, expressed in terms of a _prior distribution_
$\mathcal{P}(\mathbf{y})$ [13].

Processing of inputs by ANNs, in turn, is also often referred to as inference,
and indeed point estimates of the posterior, such as the posterior mean,
correspond to optimal estimates for typical loss functions such as the mean
squared error [14]. For decades, a popular and highly effective architecture
for performing inference has been that of deep convolutional neural networks
(CNNs), which originally drew inspiration from sensory processing in
biological systems[15]. It is then perhaps not surprising that for low level
perceptual tasks that have been preserved over long evolutionary time-scales,
optimized CNNs provide the currently best available predictors of mean
cortical neural responses [16].

There are however scenarios where point estimates are not enough. In
particular, optimal information fusion relies on knowledge about the
uncertainty of the sources [13]. A classical example of this in perceptual
tasks is multi-sensory integration, but its implications are not limited to
neuroscience. In the healthcare domain, where ML methods are increasingly
employed, diagnostics usually involve integrating results from multiple tests
and imaging modalities. While deep CNNs usually achieve high levels of
precision, they often suffer from calibration problems, resulting in
inaccurate uncertainty reports [17]. This problem becomes even more evident
for out of distribution examples where models may report a high degree of
certainty while systematically failing in their predictions. It is then
crucial to develop models that can represent uncertainty appropriately, and
the research community has been working on this problem for a number of years.
The dropout technique, which was initially designed as a regularization
method, for instance, has been proposed as a way to estimate model
uncertainty, which can be later used to compute output uncertainty in
regression and classification problems [18]. Biological neural networks are
also a source of inspiration for finding adequate solutions to this problem.
In this sense, a recurrent neural network mimicking the structure of the
primary visual cortex has been shown to represent full posterior distributions
via sampling, while at the same time displaying characteristic features of
cortical dynamics, such as stimulus dependent transient responses and
oscillations [8]. This approach, however required knowledge about the
generative model describing the likelihood function which may not always be
available.

(a)

(b)

Fig. 1: Expected behaviour of the posteriors over latent variables for low
level features in natural images.
a Three images (red, purple and blue frames) are presented at varying
contrast. As contrast increases, the inferred posteriors for each image
(colored filled lines) progressively deviate from the prior and separate from
one another. The prior is shown in dotted line. b Corresponding signal mean,
signal variance and noise variance as a function of contrast in this scenario.

In this context, Variational Autoencoders (VAEs) [19] become a valuable tool
to learn useful latent representations in an self-supervised fashion. These
models consist of an encoder network which maps inputs $\mathbf{x}$ to the
parameters of latent distributions $\mathcal{P}(\mathbf{y}|\mathbf{x})$ (the
_inference model_), and a decoder network which maps samples from the latent
representations $\mathbf{y}$ back to the original input space
$\mathcal{P}(\mathbf{x}|\mathbf{y})$ (the _generative model_). While the
original focus of these models has been to find well behaved and usually low-
dimensional latent representations, improving on those from deterministic
autoencoders, they also provide a means to simultaneously learn a generative
and an inference model from the statistics of the data.

Here we focus on visual inference tasks and explore settings relevant for
perception in biological and artificial systems. Classical work of Ref. [10]
has shown how sparse coding can explain natural receptive fields and
orientation selectivity in the cortex. Sparsity, in turn can be induced in
VAEs by employing Laplacian (rather than more standard Gaussian) priors,
resulting in sparse VAEs [20]. Indeed, when trained on the statistics of
natural images, VAEs with Laplacian (but not Gaussian) priors have been shown
to evolve latent representations closely mimicking natural receptive fields
and orientation selectivity in cortical visual processing[21].

As it happens in other domains of machine and deep learning, the optimization
of a VAE does not guarantee that uncertainty, in this case of the inferred
latent variables, will be represented in a meaningful way[17], and that
expected basic properties of the posterior distributions will be preserved.
From a Bayesian perspective it is expected that in a well behaved inference
model[22]:

  1. 1.

When no information is given by an observation, the inferred posterior is
equal to the prior.

  2. 2.

As observations become progressively more informative:

    1. 2.1

The posteriors corresponding to different images become increasingly distinct
from the prior, resulting in an increasing average distance between the
prior’s mean and the posteriors’ means.

    2. 2.2

The posteriors corresponding to different images also become increasingly
distinct from each other, resulting in an increasing variance of the
posteriors’ means.

    3. 2.3

The variance of individual posteriors (capturing uncertainty) decreases.

We note that there are two different sources of variance in latent
representations. The first refers to how the posterior mean changes with the
input, which we will henceforth refer to as the _signal variance_. The second
refers to the remaining uncertainty after having observed a given stimulus,
which we will call _noise variance_. Similarly, we call _signal mean_ the
distance between the prior’s mean and the posteriors’ mean. The _noise mean_
will be zero in our models, as habitually assumed for VAEs.

If we take for instance the case of natural images, contrast is one of the
features that regulates how informative an image is. In this context, the
posterior inferred from images with zero contrast should match the prior
distribution (Fig 1(a)). Moreover, one would expect an increase in signal mean
and variance for increasing contrast (Fig 1(b) top and middle). Finally,
uncertainty expressed in terms of noise variance should diminish for
increasing contrast (Fig 1(b) bottom). Indeed, cortical recordings in animals
presented with natural visual stimuli, have shown at the neural level how
internal models progressively adapt during development to the statistics of
natural stimuli, approximating statistically optimal representations as the
visual system matures[22]. Here we seek to test if standard VAE models, or
variations thereof, can capture this behavior and their inference modules be
used as models of cortical processing beyond mean responses. We also want to
assess and improve the quality of latent uncertainty representations in
variational models trained for traditional computer vision tasks.

To that end we audited multiple VAE models with different priors and
architectures in a variety of settings in terms of their capacity to represent
uncertainty. First, we show that a sparse VAE trained on natural images does
not necessarily present all of properties enumerated above. Concretely, we
observe an increase in uncertainty with contrast, and no convergence to the
prior for zero contrast. In order to fix this behaviour we propose a
modification to VAEs, giving rise to the Explaining-Away Variational
Autoencoder (EA-VAE). Our model incorporates a global multiplicative latent
variable $z$ in the generative model, whose posterior distribution is
simultaneously inferred with the traditional latent variables $\mathbf{y}$ by
the encoder (cf. Fig. 2 left and right). This was inspired by the Gaussian
Scale Mixture model (GSM), which has been shown to capture basic statistics of
natural image patches [23] and has previously been used to explain behavioral
and neural data related to both stationary and dynamic responses in visual
perception[24, 8]. In these models, the simultaneous high activation of a
large number of local variables in the VAE can be replaced in a simpler
fashion by an increase of a single multiplicative variable. Hence the term
explaining-away, in our EA-VAEs. This is closely related to the principle of
Occam’s razor [11] which favors choosing the simpler explanation to two
alternative competing hypotheses. We compared our model to classical VAEs in
the domain of natural images and other canonical computer vision domains such
as handwritten characters[25] and x-ray medical images[26]. Based on these
experiments, we systematically show that, contrary to standard VAEs, our model
is able to represent uncertainty in a meaningful and useful fashion.

Fig. 2: Comparison of the inference and reconstruction process in the VAE
(left) and EA-VAE (right). In the VAE a single pool of latent variables
$\mathbf{y}$ is inferred, while in the EA-VAE there is an additional global
latent variable $z$ which acts multiplicatively on the mixture of local
latents $\mathbf{y}$.

## Results

### Inference on natural images

We started by training a standard sparse VAE on gray-scale square patches of
natural images (Methods and Fig 2) and compared the resulting model to an EA-
VAE trained in the same way. We found that both the and EA-VAE could be
successfully trained on this task, showing a asimilar capacity to reconstruct
natural images (Fig 3(a)). Delving into the latent layer representations, we
found that, as had been shown for the VAE [21], the receptive fields of
neurons in the latent layer of EA-VAE also segregated into two types. One
subset was composed of units with localized receptive fields, which were
sensitive to specific orientations, while the remaining units presented noisy
and uninformative filters (Fig 3(b)). The system recruits as many units as it
needs given the effective dimensionality of the image space (Methods).
Importantly, the receptive fields of informative neurons in the EA-VAE
resemble those found in primary visual cortex [27] when employing a Laplacian
(but not a Gaussian) prior, in line with previous reports for VAEs[21], and
consistent with classical findings from sparse coding[10].

(a)

(b)

(c)

Fig. 3: a Test image patches and their respective reconstruction through the
trained VAE and EA-VAE. b Latent receptive fields in the VAE and EA-VAE. c
Statistical moments of the inferred latent posteriors in the optimized VAE
(left) and EA-VAE (right), as a function of patch contrast (Supporting
Information). Noise variance for the prior is shown in dotted black lines.

#### Uncertainty representation in the latent space

Next we analysed the properties of the VAE’s posterior representations, by
considering them in terms of their signal mean, signal variance, and noise
variance (Methods). Given the receptive fields previously shown, we know these
latent variables each represent the probability of a given orientation being
locally present in a portion of the image. As expected, we found that as
contrast increases, and different images become more distinct, the signal mean
and signal variance increase as well. This is a result of the centers of
individual posteriors separating from the prior and each-other, as expected
(Fig 3(c), left column, top and middle). Surprisingly however, noise variance,
encoding the remaining posterior uncertainty about the presence of each
orientation, does not decrease with contrast in the VAE (Fig 3(c), left
column, bottom). Moreover, as contrast tends to zero, and local orientations
become progressively harder to distinguish, uncertainty increasingly deviates
from that of the prior. We found that in this case the posterior of a blank
image (contrast zero), had a mean noise variance of 0.17, strongly deviating
from the unitary variance of prior.

EA-VAE’s posteriors also show increasing signal mean and variance as a
function of contrast, as the standard VAE. Critically, in contrast to VAEs,
EA-VAEs report a decreasing level of uncertainty for higher contrasts. We
found low noise variance for high contrast patches, which increases as
contrast becomes lower, finally converging to the prior noise variance at zero
contrast (Fig 3(c), right column). Indeed, the blank image mean uncertainty
has unitary standard deviation, matching the prior in the case of the EA-VAE,
but not in the VAE (Cf. Fig 3(c), bottom row).

### Inference in classical computer vision problems

To explore whether the observed properties of VAEs and EA-VAEs regarding
uncertainty representation were unique to the domain of natural images or may
extend to other domains, we next studied the MNIST [25] handwritten numbers
and ChestMNIST [26] chest x-rays datasets. Here again both types of models
were trained on each dataset and compared in a set of experiments exploring
uncertainty in a variety of scenarios.

#### Uninformative observations

In the case of MNIST and ChestMNIST, images have standardized contrast, which
means we cannot use this variable as a proxy for uncertainty as in the case of
natural images. We will instead resort to a number of experimental
configurations to show different properties of the proposed model. Firstly,
instead of a zero contrast image, we take the pixel-wise average image across
the dataset as an uninformative image in order to compare the models. As
expected, this image resembles a blurry, unrecognizable digit for MNIST (Fig
4(a)) and a slightly blurry x-ray for ChestMNIST (Fig 4(b)). When presented
with these averaged images, both the VAE and EA-VAE inferred a posterior with
mean close to that of the prior(not shown). The differences once again arise
when comparing uncertainty estimates. For the VAE trained on MNIST, the
average uncertainty, measured as the noise standard deviation of the posterior
along latent dimensions, was $u_{{}_{\text{VAE}}}$=0.20. This number was
comparable with those of the real handwritten digits, and considerably below
the unit standard deviation of the prior (Fig 4(a)). Notably, in the EA-VAE
the average uncertainty for this image was $u_{{}_{\text{EA-VAE}}}$=0.47,
doubling that of the VAE and halving the distance to that of the prior. While
less marked, a similar tendency emerged between the VAE and EA-VAE trained on
ChestMNIST (Fig 4(b)). Note that the average x-ray image resembles actual data
samples more closely than the average digit image, which may explain why the
effect is less pronounced.

(a)

(b)

Fig. 4: Latent uncertainty for uniformative images in the a MNIST and b
ChestMNIST domains. The distribution of reported noise standard deviation
(std.) in the VAE (green) and EA-VAE (orange) are presented as violin plots.
The uncertainty of an uninformative image from each domain is also presented
(full dot). Prior uncertainty is shown as a reference (dashed line).

(a)

(b)

(c)

Fig. 5: Latent uncertainty for morphed and corrupted images. a A
representative MNIST image of digit ‘zero’ is continuously morphed into a
‘one’ image. A sketch of the expected behaviour of the reported uncertainty is
presented below. b Results of the numerical experiments for the morphed
digits. Reported uncertainty (black line) for the VAE (left) and EA-VAE
(right) for 3 different combinations of digits. In red, quadratic function fit
employed to systematically assess the location the maximum. c Reported
uncertainty by the VAE (green) and EA-VAE (orange) for increasingly corrupted
ChestMNIST images either by Gaussian blurring (top) or additive pixel noise
(bottom).

#### Morphing and uncertainty

When dealing with data with categorical labels indicating classes such as
digits, one can compare the uncertainty for images belonging to one of the
classes with respect to continuously morphed images which lay in between two
of them. One would expect minimal uncertainty, that is low noise standard
deviation, at the images corresponding to actual digits, and higher noise
deviation for the in-between images. We conducted this numerical experiment
for each of the 45 pairs of different digits and plotted the uncertainty along
the interpolation path looking for its maximum (Fig 5(a)). We found that the
maximum fell inside a central window of half the interpolation path length in
all cases for the EA-VAE; while VAE satisfied this only in 17 cases (38%) (Fig
5(b) cf. left and right columns). This result remained consistent when varying
the window width (Details on computation and results of this experiment
available in the Supporting Information text).

#### Image corruption and uncertainty

We next looked at uncertainty representation as x-ray images fromChestMNIST
were progressively corrupted, rendering them increasingly uninformative. We
did so in two ways: by applying a Gaussian blur (Fig 5(c) top) and by adding
Gaussian pixel noise (Fig 5(c) bottom) (Details on computation in the
Supporting Information text). In both cases the EA-VAE reported increasing
uncertainty as the image quality deteriorated. This starkly contrasts with the
behavior of VAEs, where latent uncertainty remained almost constant at low
values even as the original image is no longer recognizable from the corrupted
one (cf. orange and green lines, Fig 5(c)).

#### Out of distribution detection

Finally, given that the EA-VAE presented a robust improvement in uncertainty
representation for uninformative, corrupted and interpolated images within
each dataset, we wondered whether it could also help with out-of-distribution
detection. Concretely we asked the question of whether models would report
particularly high uncertainty values when shown images from a entirely
different dataset to that used for training. Using the three previously
employed datasets and their corresponding VAE and EA-VAE trained models, we
compared internal representations and reported uncertainties for within
distribution vs. out of distribution images (Fig 6).

When looking at the models trained with MNIST and tested on ChestMNIST and
natural images we observe two important differences between the VAE and the
EA-VAE. Firstly, when analyzing the internal representation of the images by
projecting the posterior mean back to pixel space using the decoder, we note
how the VAE tends to interpret x-ray images and natural images as digits (Fig
6(a)). This is not the case with the EA-VAE, which reconstructs blurry images,
while also evidencing a marked increase in the reported uncertainty (Fig
6(b)). Indeed, the EA-VAE reports uncertainty levels much closer to the prior
when presented with out-of-distribution examples (further analysis in the
Supporting Information).

(a)

(b)

Fig. 6: Internal representations of the models and reported uncertainty for
both in-distribution and out-of-distribution examples when trained on MNIST.
a, examples of images from the same domain (blue panel) and different domains
(magenta panel) as inputs, and reconstructions by the VAE (middle row) and EA-
VAE (bottom row) respectively. b, uncertainty quantified by the posterior
noise standard deviation (std.) for individual within-distribution or out-of-
distribution images for both the VAE (green) and EA-VAE (orange). As a
reference, a dashed black line indicating prior uncertainty is presented.

## Discussion

The use of deep learning methods as modelling tools in neuroscience has proven
highly effective in recent years [16], and goal driven approaches in
particular continue to push boundaries in our understanding of the links
between perception and neural responses [2, 3, 4]. From a theoretical point of
view, the Bayesian framework provides a natural means to study optimal
perception in the brain [12]. Since exact Bayesian inference is often
intractable, two main approaches have consolidated over the years in order to
perform approximate inference beyond single point estimates of the posterior:
variational inference and sampling-based inference. While uncertainty
representation has been studied in the past in recurrent models of the cortex
implementing sampling-based inference [8], to our knowledge there is still an
important gap when it comes to a frequently used family of variational models,
such as VAEs. This is highly relevant given the tendency of deep CNNs to
suffer from calibration problems, resulting in inaccurate uncertainty reports
[17].

In this work, we went beyond point estimates, and dealt with representations
of latent uncertainty in models of perception obtained by variational deep
learning methods. Concretely, we studied VAEs both in the context of natural
image perception and computer vision. We present substantial evidence that
posteriors inferred by VAEs systematically struggle to represent uncertainty
in a meaningful way in cases of image corruption, interpolation, and out-of-
distribution detection. Motivated by generative models which employ
multiplicative variables such as the GSM [23], we incorporated a global
scaling latent variable to VAEs. In these models, the simultaneous high
activation of a large number of local variables in the VAE can be _explained
away_ by an increase in the global multiplicative variable in the EA-VAE. The
result of this modification is extremely robust. Results show that EA-VAEs
systematically outperform standard VAEs when it comes to uncertainty
representations. EA-VAE posteriors tend to more strongly revert to the prior
when images become less informative by either corruption or interpolation, and
report high uncertainty for out-of distribution examples.

Moreover, as shown in Ref. [8], artificial neural systems optimized for full
Bayesian inference beyond point estimates, for models with explaining away
variables such as the GSM, result in solutions with mean responses which
precisely recover divisive normalization, a canonical computational mechanism
in the cortex[28]. It is hence possible that as with divisive normalization in
the brain, explaining away variables in artificial neural architectures may
play an important role in multiple domains. Indeed, although originally
motivated by the study of cortical representations for visual perception in
the context of natural images, we have shown equally beneficial properties in
standard computer vision datasets such as MNIST and ChestMNIST. While further
work is required to assess these benefits in other non-visual modalities, we
believe EA-VAEs constitute a promising model type for a wide range of
scenarios. Inference models which can readily provide a meaningful
representation of uncertainty are crucial in a plethora of fields such as the
medical domain where optimal information fusion relies on weighing multiple
information sources by their relative uncertainties [29].

## Methods

### Variational Autoencoders

To obtain approximate solutions to the process of Bayesian probabilistic
inference, we trained VAEs to simultaneously learn the generative model
through a decoder, and the inference model through an encoder. For any input
$\mathbf{x}\in\mathbb{R}^{M}$, the encoder outputs a variational approximation
of the posterior $p(\mathbf{y}|\mathbf{x})$ through
$q_{\phi}(\mathbf{y}|\mathbf{x})$, where $\mathbf{y}\in\mathbb{R}^{D}$ is the
vector of latent features with an assumed prior $p(\mathbf{y})$. After taking
a sample $\mathbf{y}$ from $q_{\phi}(\mathbf{y}|\mathbf{x})$, an output
$\hat{\mathbf{x}}$ is reconstructed through the decoder by the approximate
likelihood $p_{\varphi}(\mathbf{x}|\mathbf{y})$ (Fig 2 left).

The objective function for learning the parameters of the encoder and decoder
is given by the Evidence Lower Bound (ELBO)[30]. This takes the form:

$\mathcal{L}_{VAE}(\mathbf{x},\varphi,\phi)=-\mathbb{E}_{q_{\phi}(\mathbf{y}|\mathbf{x})}\left[\log
p_{\varphi}(\mathbf{x}|\mathbf{y})\right]+\beta_{1}D_{KL}(q_{\phi}(\mathbf{y}|\mathbf{x})||p(\mathbf{y}))$
(1)

where the first term expresses the reconstruction error, and a second term
acts as regularization, penalising differences between the prior distribution
and the learnt approximate posterior of the latent variables. The $\beta_{1}$
hyperparameter determines the prior regularization weight.

### Explaining-away VAE

The EA-VAE incorporates to the standard VAE a multiplicative scalar latent
variable $z\in\mathbb{R}$, independent from the latent variables $\mathbf{y}$.
That is, any observation $\mathbf{x}$ is now assumed to have been generated by
combining spatial a set of latent features $\mathbf{y}$ via a function
$\mathbf{f}(\mathbf{y})$, further multiplied by a scalar variable $z$, and
corrupted by additive white noise $\boldsymbol{\eta}$:

$\mathbf{x}=z\,\mathbf{f}(\mathbf{y})+\boldsymbol{\eta}.$ (2)

Therefore, for an input $\mathbf{x}$, two separate encoders parameterize the
inference process through the approximate posteriors
$q_{\phi}(\mathbf{y}|\mathbf{x})$ and $q_{\psi}(z|\mathbf{x})$. For
reconstruction, both $\mathbf{y}$ and $z$ are sampled from their respective
posteriors before decoding. Note that $\mathbf{f}$ is analogous to the decoder
from the classic VAE architecture. The difference is that in the EA-VAE,
$\mathbf{f}(\mathbf{y})$ is further multiplied by $z$ before being corrupted
by noise, obtaining the reconstructed output $\hat{\mathbf{x}}$ (Fig 2 right).

Since $\mathbf{y}$ and $z$ are independent random variables, we write the
objective function for the EA-VAE as:

$\mathcal{L}_{EA-
VAE}(\mathbf{x},\varphi,\phi)=-\mathbb{E}_{q_{\phi}(\mathbf{y}|\mathbf{x})}\left[\log
p_{\varphi}(\mathbf{x}|\mathbf{y},z)\right]+\beta_{1}D_{KL}(q_{\phi}(\mathbf{y}|\mathbf{x})||p(\mathbf{y}))+\beta_{2}D_{KL}(q_{\psi}(z|\mathbf{x})||p(z)).$
(3)

### Trained models

#### VAEs and EA-VAEs for inference on natural images

(a)

(b)

(c)

(d)

Fig. 7: a Trained VAE model for inference on natural images. The encoder’s
input of size $M=1600$ is followed by two fully connected hidden layers with
softplus activation funtion. Two separate fully connected layers of $D=1800$
units encode the mean and variance of $q_{\phi}(\mathbf{y}|\mathbf{x})$. One
fully connected layer decodes the latent sample into a reconstructed patch. b
Trained EA-VAE model for inference on natural images. The described VAE model
in a was combined with an independent encoder composed of a fully connected
hidden layer with softplus activation that takes as input the patch’s pixel
intensity standard deviation $\sigma_{pix}^{(\mathbf{x})}$. Following that,
another fully connected layer encodes the variance of
$q_{\psi}(z|\mathbf{x})$. The posterior sampled $z$ value is globally
multiplied with the decoder’s output, resulting in the final reconstructed
patch. c VAE and d EA-VAE models trained on MNIST and ChestMNIST datasets.

In this work we focus on inference over gray-scale images. That is, the input
$\mathbf{x}$ corresponds here to a vector of intensity values for each of the
pixels in an image. We first implemented and trained a Sparse VAE [20]
identical to the one described in [21] (Fig 7(a)). In order to foster
sparsity, this model employs multivariate uncorrelated Laplace distributions
to parameterize the posteriors:
$q_{\phi}(\mathbf{y}|\mathbf{x})=\textrm{Laplace}(\mathbf{y};\boldsymbol{\mu}^{(\mathbf{x})},\boldsymbol{\sigma}^{2(\mathbf{x})})$,
with $\boldsymbol{\mu}^{(\mathbf{x})}$ and
$\boldsymbol{\sigma}^{2(\mathbf{x})}$ the respective mean and variance
vectors. These vectors are computed through the encoder for each specific
image $\mathbf{x}$. The latent prior $p(\mathbf{y})$ is defined from the same
family, with zero mean and unit variance for all spatial latent dimensions.
The likelihood for a given x is parameterized by a multivariate uncorrelated
Normal distribution
$p_{\varphi}(\mathbf{x}|\mathbf{y})=\mathcal{N}(\mathbf{x};Dec(\mathbf{x}),\boldsymbol{\eta}^{2})$
with $Dec(\mathbf{x})$ the mean (computed as the decoder output); and
$\boldsymbol{\eta}^{2}$ the variance (equal for all images). Thus, Eq. 1 takes
the form of

$\mathcal{L}_{VAE}(\mathbf{x},\varphi,\phi)=\sum_{i=1}^{M}(\hat{\mathrm{x}}_{i}-\mathrm{x}_{i})^{2}\\\
+\beta_{1}\sum_{j=1}^{D}\left(-1+|\mu^{(\mathbf{x})}_{j}|-\log{(\frac{\sigma^{(\mathbf{x})}_{j}}{\sqrt{2}})}+\frac{\sigma^{(\mathbf{x})}_{j}}{\sqrt{2}}e^{-\frac{\sqrt{2}|\mu^{(\mathbf{x})}_{j}|}{\sigma^{(\mathbf{x})}_{j}}}\right)$
(4)

The EA-VAE preserves the previous architecture but incorporates a scalar
multiplicative latent variable $z$. For natural images patches, this variable
is drawn from a Gamma distribution and represents in our model the contrast of
the patch, as in previous work modelling natural image patches employing
GSMs[31] (Fig 7(b)). A Gamma distribution can be parametrized in terms of the
shape and scale parameters $k$ and $\theta$ respectively. Here we have assumed
for simplicity that all posteriors
$q_{\psi}(z|\mathbf{x})=\textrm{Gamma}(z;k,\theta^{(\mathbf{x})})$ have a
fixed shape parameter $k=2$. Moreover, the proposed model infers the scale
parameter for each image from the across-pixel standard deviation of
intensities in $\mathbf{x}$: $\sigma_{pix}^{(\mathbf{x})}$. Since $k$ is
fixed, the posterior’s mean, mode and variance depend only on the inferred
scale parameter $\theta$. For the prior distribution, the shape parameter was
also fixed at $k=2$ and the scale parameter $\theta=\frac{1}{\sqrt{2}}$ was
chosen such that the prior had unit variance $s^{2}=k\theta^{2}=1$. $k=2$ was
chosen as the smallest integer number that qualitatively achieved a similar
distribution to the dataset’s $\sigma_{pix}$ distribution. Computationally,
the choice of a fixed integer $k$, simplifies the sampling procedure, since
one can always sample from the same standard Gamma function and then
reparametrize by the scale parameter. This enables simple gradient
computations across the sampling process. Furthermore, a constant shape
parameter k means that the posterior standard deviation is always proportional
to the posterior mean. This satisfies Weber’s Law of perception[32], where the
estimation error of a quantity is constant fraction of its mean. Under these
assumptions, Eq.3 can be expressed as

$\displaystyle\mathcal{L}_{EA-VAE}(\mathbf{x},\varphi,\phi)=$
$\displaystyle\sum_{i=1}^{M}(\hat{\mathrm{x}}_{i}-\mathrm{x}_{i})^{2}+\beta_{1}\sum_{j=1}^{D}\left(-1+|\mu^{(\mathbf{x})}_{j}|-\log{(\frac{\sigma^{(\mathbf{x})}_{j}}{\sqrt{2}})}+\frac{\sigma^{(\mathbf{x})}_{j}}{\sqrt{2}}e^{-\frac{\sqrt{2}|\mu^{(\mathbf{x})}_{j}|}{\sigma^{(\mathbf{x})}_{j}}}\right)$
(5)
$\displaystyle+\beta_{2}2(-\log\sqrt{2}-\log\theta^{(\mathbf{x})}+\sqrt{2}\,\theta^{(\mathbf{x})}-1),$

#### VAEs and EA-VAEs for handwritten characters and medical images

In contrast to the natural image models, the models trained on MNIST
handwritten digits dataset and ChestMNIST x-ray images consist of an encoder
and decoder composed by both fully-connected and convolutional layers (Fig
7(c) and 7(d)). Moreover, as is more standard in the ML field, the inferred
posteriors are parameterised as Multivariate uncorrelated Normal distributions
$q_{\phi}(\mathbf{y}|\mathbf{x})=\textrm{Normal}(\mathbf{y};\boldsymbol{\mu}^{(\mathbf{x})},\boldsymbol{\sigma}^{2(\mathbf{x})})$,
with $\boldsymbol{\mu}^{(\mathbf{x})}$ and
$\boldsymbol{\sigma}^{2(\mathbf{x})}$ the respective mean and variance
vectors. In the case of the EA-VAE we also have a multiplicative $z$ latent
viariable, whose posterior is given by
$q_{\psi}(z|\mathbf{x})=\textrm{Normal}(z;\nu^{(\mathbf{x})},\xi^{2(\mathbf{x})})$.
Priors $p(\mathbf{y})$ and $p(z)$ come from the same family as their
respective posteriors but with zero mean and diagonal unit covariance. As
before, the likelihood function assumes uncorrelated normally distributed
noise.

### Datasets

#### Inference on Natural Images

We used the same dataset as in [21], composed of natural image patches of
$40\times 40$ pixels each[33] ($6.4\times 10^{5}$ images for training and
validation and $6.4\times 10^{4}$ for test). The first set was further divided
into 80$\%$ for training and 20$\%$ for validation purposes. Image patches
already had already been preprocessed, whereby $100(1-\frac{\pi}{4})\%$ of the
higher frequency PCA components had been discarded. By means of this
filtering, the effective dimensionality of the image space had been reduced to
1256. Note that the patch size remains unaltered. Additionally, the
distribution of pixel intensities was rescaled such that it’s mean $\pm$ 3
deviations fit the range [0,1], to have the same range for all domains. This
was important for the out-of-distribution evaluations. Finally, the mean pixel
intensity of each patch was subtracted before the image enters the encoder
network.

#### Inference on the MNIST and ChestMNIST datasets

VAEs and EA-VAEs were trained using the handwritten numbers MNIST[25] and
x-ray images ChestMNIST[26] datasets separately.

Each original $28\times 28$ pixel image was resized to $32\times 32$ due to
the 3 layered convolutional encoder and decoder architecture, and the pixel
intensity distribution of the whole dataset was rescaled to fit the range
[0,1].

## Acknowledgments

This work was supported by Argentina’s National Scientific and Technical
Research Council (CONICET), which covered the salaries of all the authors. The
authors gratefully acknowledge NVIDIA Corporation for the GPU computing, and
the support of Universidad Nacional del Litoral (Grants CAID-
PIC-50220140100084LI, 50620190100145LI), ANPCyT (PICT-PRH-2019-00009,
2022-03-00593).

## References

  * [1] Anthony Zador, Sean Escola, Blake Richards, Bence Ölveczky, Yoshua Bengio, Kwabena Boahen, Matthew Botvinick, Dmitri Chklovskii, Anne Churchland, Claudia Clopath, et al.  Catalyzing next-generation artificial intelligence through neuroai.  Nature communications, 14(1):1597, 2023.
  * [2] Daniel LK Yamins, Ha Hong, Charles F Cadieu, Ethan A Solomon, Darren Seibert, and James J DiCarlo.  Performance-optimized hierarchical models predict neural responses in higher visual cortex.  Proceedings of the national academy of sciences, 111(23):8619–8624, 2014.
  * [3] Alexander JE Kell, Daniel LK Yamins, Erica N Shook, Sam V Norman-Haignere, and Josh H McDermott.  A task-optimized neural network replicates human auditory behavior, predicts brain responses, and reveals a cortical processing hierarchy.  Neuron, 98(3):630–644, 2018.
  * [4] Blake A. Richards, Timothy P. Lillicrap, Philippe Beaudoin, Yoshua Bengio, Rafal Bogacz, Amelia Christensen, Claudia Clopath, Rui Ponte Costa, Archy de Berker, Surya Ganguli, Colleen J. Gillon, Danijar Hafner, Adam Kepecs, Nikolaus Kriegeskorte, Peter Latham, Grace W. Lindsay, Kenneth D. Miller, Richard Naud, Christopher C. Pack, Panayiota Poirazi, Pieter Roelfsema, João Sacramento, Andrew Saxe, Benjamin Scellier, Anna C. Schapiro, Walter Senn, Greg Wayne, Daniel Yamins, Friedemann Zenke, Joel Zylberberg, Denis Therien, and Konrad P. Kording.  A deep learning framework for neuroscience.  Nature Neuroscience, 22(11):1761–1770, Nov 2019.
  * [5] Joel Dapello, Tiago Marques, Martin Schrimpf, Franziska Geiger, David Cox, and James J DiCarlo.  Simulating a primary visual cortex at the front of cnns improves robustness to image perturbations.  In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 13073–13087. Curran Associates, Inc., 2020.
  * [6] Tom Macpherson, Anne Churchland, Terry Sejnowski, James DiCarlo, Yukiyasu Kamitani, Hidehiko Takahashi, and Takatoshi Hikida.  Natural and artificial intelligence: A brief introduction to the interplay between ai and neuroscience research.  Neural Networks, 144:603–613, 2021.
  * [7] A Emin Orhan and Wei Ji Ma.  A diverse range of factors affect the nature of neural representations underlying short-term memory.  Nature neuroscience, 22(2):275–283, 2019.
  * [8] Rodrigo Echeveste, Laurence Aitchison, Guillaume Hennequin, and Máté Lengyel.  Cortical-like dynamics in recurrent circuits optimized for sampling-based probabilistic inference.  Nature neuroscience, 23(9):1138–1149, 2020.
  * [9] Daniel L. K. Yamins and James J. DiCarlo.  Using goal-driven deep learning models to understand sensory cortex.  Nature Neuroscience, 19(3):356–365, Mar 2016.
  * [10] Bruno A. Olshausen and David J. Field.  Emergence of simple-cell receptive field properties by learning a sparse code for natural images.  Nature, 381(6583):607–609, Jun 1996.
  * [11] David JC MacKay.  Information theory, inference and learning algorithms.  Cambridge university press, 2003.
  * [12] David C. Knill and Alexandre Pouget.  The bayesian brain: the role of uncertainty in neural coding and computation.  Trends in Neurosciences, 27(12):712–719, 2004.
  * [13] David J. C. MacKay.  Information Theory, Inference, and Learning Algorithms.  Copyright Cambridge University Press, 2003.
  * [14] Edwin T Jaynes.  Probability theory: The logic of science.  Cambridge university press, 2003.
  * [15] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [16] Santiago A Cadena, George H Denfield, Edgar Y Walker, Leon A Gatys, Andreas S Tolias, Matthias Bethge, and Alexander S Ecker.  Deep convolutional models improve predictions of macaque v1 responses to natural images.  PLoS computational biology, 15(4):e1006897, 2019.
  * [17] Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q Weinberger.  On calibration of modern neural networks.  In International conference on machine learning, pages 1321–1330. PMLR, 2017.
  * [18] Yarin Gal and Zoubin Ghahramani.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  In international conference on machine learning, pages 1050–1059. PMLR, 2016.
  * [19] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes, 2022.
  * [20] Gabriel Barello, Adam S. Charles, and Jonathan W. Pillow.  Sparse-coding variational auto-encoders.  bioRxiv, 2018.
  * [21] Ferenc Csikor, Balázs Meszéna, Bence Szabó, and Gergő Orbán.  Top-down inference in an early visual cortex inspired hierarchical variational autoencoder, 2022.
  * [22] Pietro Berkes, Gergő Orbán, Máté Lengyel, and József Fiser.  Spontaneous cortical activity reveals hallmarks of an optimal internal model of the environment.  Science, 331(6013):83–87, 2011.
  * [23] Martin J Wainwright and Eero Simoncelli.  Scale mixtures of gaussians and the statistics of natural images.  In S. Solla, T. Leen, and K. Müller, editors, Advances in Neural Information Processing Systems, volume 12. MIT Press, 1999.
  * [24] Dylan Festa, Amir Aschner, Aida Davila, Adam Kohn, and Ruben Coen-Cagli.  Neuronal variability reflects probabilistic inference tuned to natural image statistics.  Nature Communications, 12(1):3635, Jun 2021.
  * [25] Y. LECUN.  The mnist database of handwritten digits.  http://yann.lecun.com/exdb/mnist/, 1998.
  * [26] Xiaosong Wang, Yifan Peng, et al.  Chestx-ray8: Hospital-scale chest x-ray database and benchmarks on weakly-supervised classification and localization of common thorax diseases.  In CVPR, pages 3462–3471, 2017.
  * [27] Alexander G Dimitrov, Melissa A Sheiko, Jonathan Baker, and Shih-Cheng Yen.  Spatial and temporal jitter distort estimated functional properties of visual sensory neurons.  J Comput Neurosci, 27(3):309–319, April 2009.
  * [28] M Carandini and DJ Heeger.  Normalization as a canonical neural computation.  Nature Reviews Neuroscience, 13(1):51, 2012.
  * [29] María Agustina Ricci Lara, Candelaria Mosquera, Enzo Ferrante, and Rodrigo Echeveste.  Towards unraveling calibration biases in medical image analysis.  In Workshop on Clinical Image-Based Procedures, pages 132–141. Springer, 2023.
  * [30] C.M. Bishop.  Pattern Recognition and Machine Learning.  Information Science and Statistics. Springer, 2006.
  * [31] G Orbán, P Berkes, J Fiser, and M Lengyel.  Neural variability and sampling-based probabilistic representations in the visual cortex.  Neuron, 92(2):530–543, 2016.
  * [32] Ernst Heinrich Weber.  De Pulsu, resorptione, auditu et tactu: Annotationes anatomicae et physiologicae… CF Koehler, 1834.
  * [33] J H van Hateren and A van der Schaaf.  Independent component filters of natural images compared with simple cells in primary visual cortex.  Proc Biol Sci, 265(1394):359–366, March 1998.

## Supporting Information

### Models’ training

All models were trained in Python 3.11.3 with Pytorch 2.0.1 on a GPU NVIDIA
TITAN V. A 64GB RAM memory was available during training. An Adam optimizar
was used with a learning rate of $10^{-3}$ and weight decay of $10^{-5}$. The
fitting step was done each time a batch of size 128 passed through the
network.

#### VAE and EA-VAE trained with van Hateren

$\beta_{1}$ hyperparameter was fixed exploratorily such that the decoder’s
reconstructions did not present visibly strong deviations from the original
input. The selected choice was of $\beta_{1}=0.015$ and $\beta_{2}=0.1$ for
the VAE and EA-VAE trained on patches of natural images. Results are robust
with respect to changes in these hyperparameters, with no abrupt qualitative
changes.

The models were trained for 10000 epochs, and early stopping was employed so
that the parameters were finally set at the epoch with best performance on the
validation images. As a reference, the time taken to achieve each respective
best epoch in this setting was of $62.2$ hours for the VAE vs $67$ hours for
the EA-VAE. That is, the improvement in uncertainty representation in the EA-
VAE comes at an expense of roughly an additional 10% computational time.

#### VAEs and EA-VAEs trained with MNIST and ChestMNIST

The regularization hyperparameters were set to $\beta_{1}=5$ and $\beta_{2}=1$
for networks trained with MNIST dataset, and of $\beta_{1}=1$ and
$\beta_{2}=1$ for the networks trained with ChestMNIST.

Although the EA-VAE has $24\%$ more parameters than the VAE, all models took
around one hour to complete the 500 epochs. The final model parameters were
set to the epoch that achieved the lowest cost function on the validation set.
We note that the use of convolutions in this case drastically reduces training
times from the case of natural images employing all fully connected layers.

### Evaluation of VAEs and EA-VAEs trained on van Hateren

#### Latent receptive fields

The receptive field $\mathbf{F}_{j}$ associated to the $j$-th latent dimension
was computed through the Spike Triggered Average (STA) technique on the test
dataset for all $j=1,...,1800$. The STA was estimated as:

$\mathbf{F}_{j}=\underset{\left\\{\mathbf{x}\right\\}}{\textrm{mean}}\,\mu_{j}^{(\mathbf{x})}\mathbf{x}$

where $\mu_{j}^{(\mathbf{x})}$ is the $j$-th component of the mean of
$\mathbf{x}$’s inferred posterior $q_{\phi}(\mathbf{y}|\mathbf{x})$.

Latent units were discerned as informative or non-informative according to
their contribution to the dataset reconstruction. When computing the relative
change in the images reconstructed when all units are active, or when a
specific unit is replaced by its mean activation throughout the ensemble, it
is seen that the units associated to orientation-sensitive filters are
clustered with the most contribution, while the second group barely
contribute. All results shown in the main text were computed using only the
informative latent dimensions.

#### Uncertainty representation in the latent space

We quantified the response of neurons for different observations and levels of
uncertainty through the signal mean, signal variance and noise variance as a
function of contrast $z$. For this, given an ensemble of patches
$\left\\{\mathbf{x}\right\\}$, we start by inferring a contrast
$z^{*(\mathbf{x})}$ for each patch. In the case of the VAE, we infer contrast
as the across-pixel standard deviation of intensities:
$z^{*(\mathbf{x})}_{{}_{VAE}}=\sigma_{pix}^{(\mathbf{x})}$, while for the EA-
VAE model we use the encoded mean of the inferred posterior distribution of
$z$ as contrast: $z^{*(\mathbf{x})}_{{}_{EA-VAE}}=k\theta^{(\mathbf{x})}$

The signal mean is defined as the average distance of the posterior’s mean to
the prior’s mean (0 in this case) over all patches that have the same inferred
contrast $z$:

$\text{SM}(z)=\underset{\left\\{\mathbf{x}|z^{*(\mathbf{x})}=z\right\\}}{\textrm{mean}}\left\|\boldsymbol{\mu}^{(\mathbf{x})}\right\|$

Signal variance is computed as the variance of the posterior’s mean of patches
of same contrast, averaged across latent dimensions:

$\text{SV}(z)=\underset{j}{\textrm{mean}}\underset{\left\\{\mathbf{x}|z^{*(\mathbf{x})}=z\right\\}}{\textrm{Var}}\mu_{j}^{(\mathbf{x})}$

Finally, noise variance is defined as the posterior’s variance vector
$\boldsymbol{\sigma}^{2(\mathbf{x})}$ averaged across latent dimensions and
observations of equal inferred contrast $z$:

$\text{NV}(z)=\underset{j}{\textrm{mean}}\underset{\left\\{\mathbf{x}|z^{*(\mathbf{x})}=z\right\\}}{\textrm{mean}}\sigma_{j}^{2(\mathbf{x})}$

These quantities were computed for the test set, where no two patches have the
same exact contrast. Hence, these were binarized with small enough bins such
that no jumps were visible in the plots, and with error bars comparable to
line-widths. The binnarized expressions are as follows:

$\text{SM}(z_{i})=\underset{\left\\{\mathbf{x}|z^{*(\mathbf{x})}\in
B_{i}\right\\}}{\textrm{mean}}\left\|\boldsymbol{\mu}^{(\mathbf{x})}\right\|$
$\text{SV}(z_{i})=\underset{j}{\textrm{mean}}\underset{\left\\{\mathbf{x}|z^{*(\mathbf{x})}\in
B_{i}\right\\}}{\textrm{Var}}\mu_{j}^{(\mathbf{x})}$
$\text{NV}(z_{i})=\underset{j}{\textrm{mean}}\underset{\left\\{\mathbf{x}|z^{*(\mathbf{x})}\in
B_{i}\right\\}}{\textrm{mean}}\sigma_{j}^{2(\mathbf{x})}$

The associated errors can be expressed as:

$\epsilon_{{}_{\text{SM}}}(z_{i})=\frac{1}{\sqrt{N_{i}}}\sqrt{\underset{\left\\{\mathbf{x}|z^{*(\mathbf{x})}\in
B_{i}\right\\}}{\textrm{Var}}\left\|\boldsymbol{\mu}^{(\mathbf{x})}\right\|}$
$\epsilon_{{}_{\text{SV}}}(z_{i})=\frac{1}{\sqrt{D}}\underset{j}{\textrm{mean}}\sqrt{\frac{2}{(N_{i}-1)}}\underset{\left\\{\mathbf{x}|z^{*(\mathbf{x})}\in
B_{i}\right\\}}{\textrm{Var}}\mu_{j}^{(\mathbf{x})}$
$\epsilon_{{}_{\text{NV}}}(z_{i})=\frac{1}{\sqrt{N_{i}}\sqrt{D}}\sqrt{\underset{\left\\{j\cup\mathbf{x}|z^{*(\mathbf{x})}\in
B_{i}\right\\}}{\textrm{Var}}\boldsymbol{\sigma}^{2(\mathbf{x})}}$

Moreover, we selected the specific subset of patches with
$\sigma_{pix}^{(\mathbf{x})}=1$ and computed these quantities for them and for
copies of this images but synthetically modifying their contrast. The
resulting moments had the same behaviour as the ones found with the whole test
set with original contrast values.

It is worth mentioning that the natural scenes dataset used in this work has
patches with varied values of pixel intensity standard deviation
$z_{pix}^{(\mathbf{x})}$. As patches with high contrast live far away from the
origin in the in the N-dimensional pixel space (since
$z_{pix}^{(\mathbf{x})}=\sqrt{\frac{1}{N-1}\sum_{i=1}^{N}(x_{i}-\bar{x})^{2}}=\sqrt{\frac{1}{N-1}}\left\|\mathbf{x}\right\|$),
it is expected that the network will have more difficulty exploring farther
regions in the pixel space than the regions where low-contrast patches live.
Therefore, in evaluation instances, patches with contrast higher than
$\sigma_{pix}=1$ were discarded.

### Evaluation of VAEs and EA-VAEs trained on MNIST and ChestMNIST

#### Morphing and uncertainty

We conducted an experiment to study the uncertainty in the VAE and EA-VAE
latent space when morphing digits from one label to another. The first step
for doing this was selecting all the test images that belong to a specific
label A and computing the average coordinate in the latent space of these (Fig
S1a top panels). The output $\hat{\mathbf{x}}_{A}$ of the decoder from that
position is a representative image of the label A (Fig S1a bottom). The same
procedure is repeated with images belonging to label B obtaining
$\hat{\mathbf{x}}_{B}$. Then, morphing is done by linearly combining
$\hat{\mathbf{x}}_{A}$ and $\hat{\mathbf{x}}_{B}$ with weight $\lambda$
varying from 0 to 1:

$\hat{\mathbf{x}}_{AB,\lambda}=(1-\lambda)\hat{\mathbf{x}}_{A}+\lambda\hat{\mathbf{x}}_{B}$

We then inferred the corresponding posteriors taking these morphed images as
input, and computed the uncertainty of the network as the standard deviations
of the posterior averaged across latent dimensions:

$\text{noise
std.}(\lambda;AB)=\underset{j}{\textrm{mean}}\,\,\sigma_{j}^{(\hat{\mathbf{x}}_{AB,\lambda})}$

Through this procedure, uncertainty curves were computed for all the
combinations of morphing from labels $0$ to $9$ (Fig S1c black curves).

To evaluate if a certain curve has the expected behaviour of being maximally
uncertain at intermediate $\lambda$ values, we fitted a second order
polynomial function (Fig S1c red curves) and checked whether it’s curvature
was negative and the maximum was located inside a window of size $L$ centered
respect to $\lambda=0.5$, for different $L$ values (Fig S1b).

#### Image corruption

The image corruption experiment was done over the complete test set. Each
image was synthetically corrupted by applying a two-dimensional Gaussian
kernel of covariance $\eta_{b}I_{2x2}$. The resulting blurred image was
encoded into the models’ latent space, where the measure of uncertainty was
computed. We compute the uncertainty of the network for a given input
$\mathbf{x}$ as the standard deviations of the posterior
$q_{\phi}(\mathbf{y}|\mathbf{x})$ averaged across latent dimensions:

$\text{noise std.}=\underset{j}{\textrm{mean}}\,\,\sigma_{j}^{(\mathbf{x})}$

This process was repeated for all images, and the average uncertainty was
computed. This was repeated for different $\eta_{b}$ values until the images
were completely blurred out ($\eta_{b}=16$ equal to half the pixel width of
the images) We repeated the experiment with pixel-noise corruption. In this
case each pixel was intensity corrupted by additive white Gaussian noise of
null mean and variance $\eta_{p}^{2}$. The whole perturbed image was then
rescaled in intensity to preserve the original range of intensities.$\eta_{p}$
values between 0 and 5 were considered, going from no perturbation to a
strongly perturbed image.

#### Out of distribution experients

To study the capacity of the trained models to associate near-to-prior
uncertainty for images from previously unseen distributions, we compared the
latent representations of images in and out of the distribution used for
training. We computed the uncertainty of the network to a given observation
the same way as in the two previous experiments. We looked at: (1) models
trained with MNIST and tested on ChestMNIST and natural images (Fig 6 on main
text); (2) models trained with ChestMNIST and tested on MNIST and natural
images (Fig 2(a)); and (3) models trained with natural images and tested on
MNIST and ChestMNIST (Fig 2(b)). In all three cases, images from the out of
distribution datasets were resized to coincide the size of images from the
training set, and rescaled in pixel intensity to match the in-distribution
pixel intensity distribution.

We found the same behaviour as shown in the main text when training models
with x-ray images and testing on natural images (Fig 2(a), rightmost block).
Yet when testing on digits, both models seem to fare comparably, with reported
uncertainties that do not approach the prior, although the reconstructed
images do not collapse to within domain examples. When training with natural
images and testing on the other two (Fig 2(b)), the distribution of reported
uncertainties, which originally had two modes for within distribution images
corresponding to more informative (high contrast) and less informative (low
contrast) images, reduces to a single mode mode for VAEs and EA-VAEs. The
asymmetry in the results when original and target domains are inverted
suggests a hierarchy in the richness of the feature space learned from the
statistics of the images.

Fig. S1: Morphing experiments in MNIST. Centers for the distribution
corresponding to each digit are first found (a). For a given window size
around the middle point between the centers (b), the fraction of
interpolations where the maximally uncertain linear combination lies within
this window is computed via quadratic fitting (c).

(a)

(b)

Fig. S2: Out of distribution experiments. Same as in Fig. 6, now when training
on ChestMNIST (a), and natural images (b)