text:
\section{Analysis}
\label{section:analysis}

\subsection{The NFW halo}
\label{section_nfw_stable}

The DM halo is constructed using a Monte Carlo sampling
technique. First, for each particle, the three position coordinates in
spherical coordinates ($r,\ \theta,\ \phi$) are generated.  $r$ is
drawn from the density profile $\rho_{\rm NFW}$ using a standard
acceptance-rejectance technique, $\phi$ and $\cos(\theta)$ are drawn
from uniform distributions over the intervals $[0,2\pi]$ and $[-1,1]$,
respectively. Next, $v_{\rm r}$, $v_{\rm \theta}$ and $v_{\rm \phi}$
are drawn from the isotropic distribution function for the NFW model,
again with an acceptance-rejectance technique. This isotropic
distribution function was constructed from the NFW density profile
using the standard Eddington formula \citep{buyle07}. For each
particle, a symmetric partner was constructed with position
coordinates $(r,\ -\theta,\ -\phi)$ and velocity coordinates $(-v_{\rm
  r}, -v_{\rm \theta}, -v_{\rm \phi})$. This drastically improved the
stability of the central parts of the halos. The very inner part of
the steep cusp of the NFW model is populated by relatively few
particles, destroying its spherical symmetry and introducing
unbalanced angular momenta. This initial deviation leads to the
ejection of particles from the cusp and triggers a more widespread
dynamical response of the DM halo, over time erasing the inner
cusp. Introducing the partner particles, cancelling out the angular
momenta and increasing the symmetry of the particles' spatial
distribution, greatly alleviates these problems. Such techniques for
constructing ``quiet'' initial conditions have been applied before
with great success, see e.g. \citep{selwood86}. The
  improvement of the stability of the DM halo in simulations with a
  ``quiet'' start over simulations without a ``quiet'' start is
  illustrated in the top panel of Fig.  \ref{fig:NFWstability} where
  the density distribution of both haloes at $z=0$ is plotted as red
  and green dots, respectively.

\begin{figure}
\centering \includegraphics[width=0.48\textwidth]{figures/compDensity.eps}
 \caption{The density profile of the N03 NFW halo for different
   simulations: in the upper panel only DM was included, in the central panel
   DM, gas is included but star formation was turned off. The bottom
   panel shows the results of a simulation with DM, gas and star
   formation. \label{fig:NFWstability}}
\end{figure}

First, to test the stability of the NFW halos, we ran several simulations for the N03 and N05 mass models:
\begin{description}
 \item[\textbf{Run 1:}] only DM
 \item[\textbf{Run 2:}] DM and gas but no star formation
 \item[\textbf{Run 3:}] DM and gas and star formation
 \end{description}
For these test simulations an $n_{\rm SF}$ of $0.1\ \rm cm^{-3}$
\citep{katz96} and $\epsilon_{\rm FB}$ of $0.1$ \citep{thornton98} was
used.

Fig. \ref{fig:NFWstability} shows the density profile of the test
simulations for the N03 mass model. From the upper panel, it is evident that the DM density of the
DM-only simulation remains stable and cusped until
the end of the simulation. The simulations presented in the middle and bottom panels,
show a clear conversion of the cusp into a core over time. Moreover, the width of the core depends on the mass of
the system, with more massive halos having larger cores.

Our simulations largely confirm the results from \cite{read05}, where a rapid removal of gas results in a conversion from cusp to core as stated first by \cite{navarro96}.
As gas cools and flows into the halo, the center of the dark matter
halo is adiabatically compressed. Without star formation, the central gas pressure builds up,
eventually stops further inflow, and even makes the gas re-expand somewhat. This re-expansion
happens rapidly enough for the DM halo to respond non-adiabatically:~the central DM density
experiences a net lowering and the cusp is transformed into a core. With star formation turned on,
feedback is responsible for a fast removal of gas from the central parts of the
DM halo, with the same effect:~a conversion from a cusp to a core.

Unlike us, \cite{governato10} found that the density threshold
  for star formation needed to be high enough for a cusp-to-core
  conversion to occur. Only for $n_{\rm SF} \gtrsim 10$~cm$^{-3}$ does
  supernova feedback lead to sufficient gas motions to flatten the
  cusp in their simulated dwarfs, which are taken from a larger
  cosmological simulation. In contrast, in our more idealized,
  initially spherically symmetric setup, even a low density threshold
  leads to sufficient gas outflow for the cusp to flatten.

\subsection{Star formation histories}
\label{section_SFH}

\begin{figure}
 \centering \includegraphics[width=0.5\textwidth]{figures/SFR.eps}
 \caption{Top panel:~the SFR of several N07 models as a function of time. Bottom panel:~the stellar
   mass as a function of time.\label{fig:SFR}}
\end{figure}

In Fig. \ref{fig:SFR}, we show the star-formation histories (SFHs) of
different realizations of the N07 mass model. Also, in table \ref{table_finalprop}, the starting time of star formation is tabulated along with the final total stellar mass. Several conclusions can be drawn:
\begin{itemize}
 \item The delay between the start of the simulation and the start of the
first star-formation event is an increasing function of $n_{\rm SF}$. This appears
logical:~it takes longer for the gas to collapse to higher densities and ignite star formation.
Comparing different mass models, star formation starts earlier in more massive models for a given
$n_{\rm SF}$. This is most likely due to the more massive models having steeper gravitational
potential wells, increasing their ability to compress the inflowing gas.
\item If $n_{\rm SF}$ is increased while  $\epsilon_{\rm FB}$ is kept fixed, more stars are formed
(e.g. going from the green to the blue curve or similarly from the cyan to the magenta curve in Fig.  \ref{fig:SFR}). This is because gas collapses
to higher densities and the feedback is no longer able to sufficiently heat and expel this gas and to interrupt star formation.
\item Related to the previous point, the SFR also becomes more rapidly varying if $n_{\rm SF}$ is increased while
$\epsilon_{\rm FB}$ is kept fixed. The reason is that in the small high-density star-forming regions, feedback can only
locally interrupt star formation during short timespans. At lower $n_{\rm SF}$, star formation is more widespread, leading
to more global behavior:~as supernovae go off, star formation can be completely halted.
\item Increasing $\epsilon_{\rm FB}$ while $n_{\rm SF}$  is kept fixed leads to a decrease in star formation (e.g. going from the blue
to the cyan curve in Fig. \ref{fig:SFR}). This is because once feedback is strong enough, it is able to extinguish
star formation, even at high gas densities.
\item The most low-mass models fail to form stars for high $n_{\rm SF}$ values. E.g.
no stars form in the N03 models for $n_{\rm SF}>0.1$~cm$^{-3}$. This is due to the masses of these models being too small for gas to collapse to densities where stars can be formed. This point is further elaborated in the next paragraph.
\end{itemize}

\subsection{Density distribution of the ISM}
\label{section_ISM}

\begin{figure}
 \centering \includegraphics[width=0.5\textwidth]{figures/N03densityComp.eps}
 \caption{The density distribution of the ISM at different times for the least massive galaxy, N03, with different density threshold and a fixed feedback efficiency of 0.1.\label{fig:N03}}
\end{figure}

In Fig. \ref{fig:N03}, the density of the ISM is plotted as a function of radius. For the N03 model in the left panel a density threshold of $0.1$~cm$^{-3}$ was used while for the model in the right panel, the density threshold was set to a value of $6$~cm$^{-3}$. The red points show the gas distribution at the moment just before the start of star formation in the case of $n_{\rm SF}$ = $0.1$~cm$^{-3}$. Since up to that moment, all models have experienced the same evolution, there is no difference between the red points in both panels.
As can be seen in the left panel, the gas density in this N03 model reaches the star-formation threshold and star formation occurs. Moreover, the influence of supernova feedback can be seen in the green and blue points, where gas expands to larger radii and lower densities after having been heated. As is clear from the right panel, for $n_{\rm SF}$ = $6$~cm$^{-3}$ the gas simply keeps falling in. It will continue to do so during the first 4 Gyr until the built-up central pressure causes the gas to re-expand again. No stars are formed during the course of this simulation.

As the density threshold is increased to higher values, star formation tends to occur more and more in small collapsed clumps. This becomes clear when comparing the panels from Figs. \ref{fig:N03} and \ref{fig:N07}. The latter shows the gas density distributions of two N07 models with $n_{\rm SF}$ = $6$~cm$^{-3}$ and $n_{\rm SF}$ = $50$~cm$^{-3}$. While the $n_{\rm SF}$ = $50$~cm$^{-3}$ model only exhibits star formation in a small number of discrete high-density clumps, the  $n_{\rm SF}$ = $6$~cm$^{-3}$ model lacks such well-defined clumps and star formation occurs more widespread.

\begin{figure}
 \centering \includegraphics[width=0.5\textwidth]{figures/N07densityComp.eps}
 \caption{The density distribution of the ISM at different times for the N07 model, with different density thresholds and a fixed feedback efficiency of 0.7. \label{fig:N07}}
\end{figure}

\subsection{Scaling relations}
\label{section_scaling_relations}
In this section we discuss the properties of each of our models and draw some conclusions regarding the influence of the  $n_{\rm SF}$ and $\epsilon_{\rm FB}$ parameters on the models. An overview of some basic properties can be found in Table \ref{table_finalprop}.
\begin{table*}
 \begin{minipage}[ht]{\linewidth}
  \begin{center}
   \caption{Final properties of our large set of simulations. Columns: (1) model number (see Table \ref{Models_table}), (2) density treshold for star formation, (3) feedback efficiency, (4) final stellar mass, (5) starting time of star formation, (6) half-light radius, (7) mean surface \textbf{brightness} within the half-light radius, (8) central one dimensional velocity dispersion, (9) mass-weighted metallicity, (10) central surface brightness, (11) S\'ersic parameter, (12) circular velocity.} \label{table_finalprop} \input{table_finalprop}
  \end{center}
 \end{minipage}
\end{table*}

\subsubsection{Half-light radius $R_{e}$}
The half-light radius, or effective radius, denoted by $R_{e}$, encloses half of
a galaxy's luminosity. In panel a.) of Fig. \ref{fig:bigFig},
$R_{e}$ is plotted as a function of the $V$-band magnitude.
The following trends can be observed in this figure:
\begin{itemize}
 \item For a fixed $n_{\rm SF}$, the effective radius varies only very slightly throughout the $\epsilon_{\rm FB}$-range and this without a clear trend between $R_{e}$ and $\epsilon_{\rm FB}$. However, for a fixed $n_{\rm SF}$ and dark-matter mass the stellar mass and consequently the luminosity decrease with increasing $\epsilon_{\rm FB}$. This is due to star formation being shut down more rapidly when feedback is more effective. As a result, galaxies tend to have higher stellar densities for smaller $\epsilon_{\rm FB}$.
 \item For a fixed $\epsilon_{\rm FB}$ of 0.1, an increase of $n_{\rm SF}$ from $0.1$ to $6$~cm$^{-3}$ results in a decrease of the effective radius. This is due to the size of the region where the SFC are fulfilled, which is much smaller for $n_{\rm SF}$ = $6$~cm$^{-3}$ than for $n_{\rm SF}$ = $0.1$~cm$^{-3}$, and the feedback is too weak to overcome this. In the case of an increase of $n_{\rm SF}$ from $6$ to $50$~cm$^{-3}$, the effective radius increases which is caused by the higher star formation peaks resulting in more supernovae explosions which redistribute the gas more efficiently.
 \item The simulations with high density threshold, $n_{\rm SF} > 0.1\ \rm cm^{-3}$, and high feedback efficiency, $\epsilon_{\rm FB} > 0.1$, have effective radii which are in agreement with the observations.
\end{itemize}
\textit{From this scaling relation we can constrain the $\epsilon_{\rm FB}$-parameter to be higher then 0.1 to produce galaxies with effective radii in agreement with observations of dwarf galaxies.}

\subsubsection{The fundamental plane}
The fundamental plane (FP) is an observed relation between the effective radius, $R_{\rm e}$, the mean surface brightness within the effective radius, $I_{\rm e}$, and the central velocity dispersion, $\sigma_{\rm c}$ of giant elliptical galaxies. It is a linear relation, given by
\begin{equation}
 \log( R_{\rm e} ) = -0.629 - 0.845 \log( I_{\rm e} ) + 1.38 \log( \sigma_{\rm c} ),
\end{equation}
between the logarithms of these quantities \citep{burstein97}. In panel b.) of Fig. \ref{fig:bigFig}, we plot the ''vertical'' deviation  of the simulated galaxies from the giant galaxies' FP.

Dwarf galaxies generally lie above the FP in this projection. This is thought to be a consequence of their having shallower gravitational potential wells than giant galaxies. This, together with
the feedback, results in more diffuse systems. Models with a high star-formation threshold in combination with a low supernova feedback turn out to be very compact. They actually populate
the FP at low luminosities. However, this region of the three-dimensional space spanned by $\log(R_{\rm e})$, $\log( I_{\rm e} )$, and $\log( \sigma_{\rm c} )$ is observed to be devoid of galaxies.
Hence, models with low stellar feedback, $\epsilon_{\rm FB}$ up to $0.3 $, and high density thresholds, $n_{\rm SF} > 0.1$~cm$^{-3}$, can be rejected.

\begin{figure*}
\begin{minipage}[ht]{\linewidth}
\begin{center}
 \centering

 \includegraphics[width=0.95\textwidth,clip]{figures/bigFig.eps}
 \caption{Some scaling relations and the surface brightness parameters as a function of the magnitude. In a.), the half-light radius $R_{\rm e}$ is plotted, b.) shows the vertical deviation of the simulated dwarf galaxies from the giant galaxies' FP, in c.) the $V-I$ color is plotted, d.) shows the iron content $[Fe/H]$. In panel e.) and f.), the S\'ersic index $n$ and central surface brightness $\mu_{0}$ are plotted. All these quantities are plotted against the $V$-band magnitude, except the FP which are plotted as a function of the $B$-band luminosity. The models with a density threshold of $6 \rm ~cm^{-3}$ and $50 \rm ~cm^{-3}$ are represented by blue-green diamonds and yellow-red triangles, respectively, where the colorscales represent a varying feedback efficiency.
 For each color, the datapoints are connected by a line showing the mass evolution of the models. In the case of $n_{\rm SF}=0.1\ cm^{-3}$, represented by the black line, the models from N03 until N09 are plotted. In the cases of higher densities, represented by the colored lines, the datapoints are from models N05-N09.
 Our models are compared with observational data
   obtained from \citet{derijcke05}, \citet{graham03}
   , LG data come from \citet{peletier93},
   \citet{peletier93}, \citet{irwin95}, \citet{saviane96},
   \citet{grebel03}, \citet{mcconnachie06},
   \citet{mcconnachie07}, \citet{zucker07}, Perseus data from
   \citet{derijcke09}, Antlia data from
   \citet{castelli08}. For the $[Fe/H]-M_{V}$ plot, data from \citet{grebel03}, \citet{sharina08} and \citet{lianou10} was used, the yellow and magenta dots represent data from dSph and dIrr galaxies, respectively. \label{fig:bigFig}}
 \end{center}
 \end{minipage}
\end{figure*}

\subsubsection{Color $V-I$}
Fig. \ref{fig:bigFig}, panel c.) shows the $V-I$ color in function of the $V$-band magnitude.
The color scatter between the different models is rather small.  The observed galaxies follow a mass-metallicity relation so the metallicity generally increases with the galaxy (stellar) mass, resulting in increasing $V-I$ values for increased galaxy mass. Within the relatively small mass range covered by the models, color
is only a very weak function of stellar mass. For a fixed feedback efficiency, when increasing the density threshold the $V-I$ also increases slightly resulting in bluer galaxies for the models with low density threshold. This is due to the effect that stars are formed in more metal enriched regions in the models with high density threshold.
When the density threshold is kept constant and only the feedback efficiency is increased the $V-I$ slightly decreases, so the models get slightly bluer due to a dilution of the gas when it is more spread out by supernovae explosions.

\subsubsection{Metallicity}
In  panel d.) of Fig. \ref{fig:bigFig} a plot of iron content [Fe/H] as a function of the $V$-band magnitude is shown. The mass-weighted value of [Fe/H] is a measure of the metallicity of a galaxy. The yellow and magenta dots represent observational data from  dwarf spheroidal and dwarf elliptical galaxies and irregular dwarf galaxies, respectively.
Some general conclusions we can take away from this figure are:
\begin{itemize}
 \item Low-mass models with low density threshold, $n_{\rm SF} \approx 0.1$~cm$^{-3}$, and low feedback, $\epsilon_{\rm FB} \approx 0.1$, keep forming stars throughout cosmic history
and do not expel enriched gas. As a consequence, they turn out to be too metal rich, compared with observed dwarf galaxies. Models with higher $n_{\rm SF}$ compare much
more favorably with the data in this respect.
  \item For a fixed $n_{\rm SF}$, increasing $\epsilon_{\rm FB}$, produces more metal poor galaxies. This is likely due to the fact that the increased feedback
  extinguishes star formation more rapidly and disperses the metal enriched gas more widely.
 \item Increasing $n_{\rm SF}$ at fixed $\epsilon_{\rm FB}$ and fixed mass, results in an increase of the metallicity and of the stellar mass when going from $n_{\rm SF} = 0.1$~cm$^{-3}$ to $n_{\rm SF} = 6$~cm$^{-3}$. A further increase of $n_{\rm SF}$ at fixed $\epsilon_{\rm FB}$, up to $n_{\rm SF} = 50$~cm$^{-3}$, has a much smaller impact on metallicity and stellar mass. The former is likely due to more vigorous star formation in less easily dispersable high density regions.
 \end{itemize}

\subsubsection{Surface brightness profiles}
We fitted a S\'ersic profile, of the form
\begin{equation}
 I(R)=I_{0}e^{-\left(\frac{R}{R_{0}}\right)^{1/n}},
\end{equation}
to the surface brightness profiles of the simulated galaxies.
The S\'ersic parameter $n$ and the central surface brightness $\mu_{0}$ are plotted respectively in the panels e.) and f.) of Fig. \ref{fig:bigFig} as a function of the $V$-band magnitude.
\begin{itemize}
 \item For a fixed $n_{\rm SF}$, when increasing the $\epsilon_{\rm FB}$, there is a weak trend for the S\'ersic parameter $n$ and the central surface brightness to decrease. More vigorous
feedback appears to result in more diffuse dwarf galaxies, as one would expect.
 \item As an echo of the $R_{\rm e}-M_{\rm V}$ relation, simulations with high density threshold, $n_{\rm SF} > 0.1\ \rm cm^{-3}$, and low feedback efficiency, $\epsilon_{\rm FB} = 0.1-0.3$, are systematically too compact, with $\mu_0 \sim 20$~mag~arcsec$^{-2}$, compared with the observations.
 \item The models with high density thresholds and strong feedback are in general agreement with the observations.
\end{itemize}

\subsubsection{The Tully-Fisher relation}
\begin{figure}
 \centering \includegraphics[width=\columnwidth]{figures/TF_FJ.eps}
 \caption{The top panel shows the Tully-Fisher relation between the circular
   velocity and the luminosity in the B-band. The full gray line shows the TF relation for early type galaxies, the dashed gray line is the TF relation of spiral galaxies as determined by \citet{derijcke07}. The lower panel shows the Faber-Jackson relation between the velocity
   dispersion and the luminosity in the B-band.}
 \label{fig:TF_FJ}
\end{figure}

Panel a.) of Fig. \ref{fig:TF_FJ} shows the B-band Tully-Fisher
relation (TFR) between the circular velocity, denoted by $V_{c}$, and the luminosity in the B-band, $L_{\rm B}$. The simulations are compared with observational data and with the Tully-Fisher relation for early-type (full gray line) and for spiral galaxies (dotted gray line) that was determined by \cite{derijcke07}.
All simulations predict that the TFR becomes substantially shallower in the dwarf regime, below luminosities of the order of $L_{\rm B} \sim 10^7\,L_{\odot, \rm B}$.  This can be seen as a consequence of the
very steep $M_{\rm star}-M_{\rm halo}$ relation in the dwarf galaxy regime (see paragraph \ref{subsubMstMh}). For a fixed $n_{\rm SF}$, an increase in feedback efficiency does not influence $V_{c}$ very much since there are so few stars that $V_c$ is set by the dark-matter halo. The effect on the stellar mass, and consequently on $L_{\rm B}$, is, however,
quite large. Therefore, increasing $\epsilon_{\rm FB}$ at fixed $n_{\rm SF}$ and dark-matter mass causes galaxies to shift leftwards in panel a.) of Fig. \ref{fig:TF_FJ}. Except for this effect, once $n_{\rm SF}$ and $\epsilon_{\rm FB}$ are raised above their minimum values of $0.1$~cm$^{-3}$ and $0.1$, respectively, there is no significant differences between the TFRs traced by the different series of models.

\subsubsection{The Faber-Jackson relation}
The Faber-Jackson (FJR) relation, plotted in panel b.) of Fig. \ref{fig:TF_FJ} is the relation between the stellar central velocity dispersion and
the luminosity in the $B$-band. The stellar central velocity dispersion is a projection of the velocity dispersion along the line of
sight. This is measured by fitting an exponential function to the dispersion profile and retaining the maximum of the function as the
central value.

From this figure we see:
\begin{itemize}
 \item For a fixed $n_{\rm SF}$, when increasing the $\epsilon_{\rm FB}$, the velocity dispersion decreases first after which it settles around a value which depends on the dark-matter mass of the model.
 \item For a fixed $\epsilon_{\rm FB}$, when increasing $n_{\rm SF}$, only a minor influence on the velocity dispersion is observed.
\end{itemize}

\subsubsection{The $M_{\rm star}$-$M_{\rm halo}$ relation} \label{subsubMstMh}
In Fig. \ref{fig:MstarMhalo}, the $M_{\rm star}$-$M_{\rm halo}$ relation of the simulations at $z=0$ is plotted.
We can make similar conclusions here as were made in the SFH section:
\begin{itemize}
 \item If $n_{\rm SF}$ is fixed, the stellar mass will decrease if the $\epsilon_{\rm FB}$ is increased. This is what was expected because with more
feedback the gas is distributed over a larger area and
the infall of the gas to the appropriate density threshold will take
longer.
 \item If $\epsilon_{\rm FB}$ is fixed, for increasing $n_{\rm SF}$, the stellar mass increases too. When feedback is very small, the gas density will stay high and the star formation will not be interrupted, resulting in a high stellar mass. The effect is smaller for higher feedback.
\end{itemize}

In Fig. \ref{fig:MstarMhalo}, our different sets of
models are found to be in agreement with the results from
the Aquila simulation where a density threshold of 10 $\rm cm^{-3}$ and a feedback efficiency of 0.7 was used.
While the initial conditions of our dwarf galaxy simulations are admittedly quite simplified, they do have high spatial resolution
and realistic implemented physics. It is therefore encouraging that they compare favorably with cosmological simulations
like the Aquila simulation, which have cosmologically well motivated initial conditions but in which dwarf galaxies are very close to the resolution limit \citep{sawala11b}.
However it is impossible by further tuning of the feedback efficiency and/or the density threshold to reproduce the trend that was derived by \cite{guo10}.

By increasing the density threshold and feedback efficiency, the stellar mass is reduced by almost two orders of magnitude, but there still remains a difference of many orders of magnitude between our simulations and the $M_{\star}$-$M_{halo}$ relation from \cite{guo10}. It is also interesting to notice that although our models do not reproduce the relation, they do have a very similar slope.

\section{Numerical details}
\label{section:numerical_details}

We use a modified version of the Nbody-SPH code {\sc Gadget-2}
\citep{springel05}. The original {\sc Gadget-2} code was extended with
star formation, feedback and radiative cooling by \cite{valcke08}.
While the initial conditions of the simulations are cosmologically
motivated (see below), we do not perform full cosmological
simulations. Our approach yields a high mass resolution at
comparatively low computational cost. Still, previous work by
\cite{valcke08}, \cite{valcke10} and \cite{schroyen11} has shown that
with this code realistic dwarf galaxies, following the known
photometric and kinematic scaling relations, can be produced.
We set up the simulations using 200,000 gas particles and 200,000 DM particles.
Depending on the model's total mass, this results in gas particle masses in the range of
$350-2,620 M_{\odot}$ and DM particle masses in the range of $1,650-12,380 M_{\odot}$.
We use a gravitational softening length of 0.03~kpc.

Our results are visualized with our own software package HYPLOT. This
is freely available from SourceForge\footnote{http://sourceforge.net/projects/hyplot/}
 and is used for all the figures
in this paper.

\subsection{Initial conditions}

Our models are set up, as in \cite{valcke08,valcke10,schroyen11}, with a spherically symmetric dark matter halo
and a homogeneous gas cloud. This gas cloud has a density of $5.55~\rho_{\rm crit}$, with $\rho_{\rm crit}$ the critical density of the universe at the halo's formation redshift, here taken to be $z_{\rm c}=4.3$. This is equivalent with a number density for the gas of 0.0011 hydrogen atoms per cubic centimeter.
We use a flat $\Lambda$-dominated cold dark matter cosmology with the following cosmological parameters:
$h = 0.71, \Omega_{\rm tot} = 1, \Omega_{\rm m} = 0.2383, \Omega_{\rm DM} = 0.1967$. The baryonic mass fraction will be the difference between $\Omega_{\rm m}$ and $\Omega_{\rm DM}$, in practice it will have a value that is 0.2115 times that of the dark-matter.
At the start of the simulations the gas particles are initially at rest, their initial metallicities are set to $10^{-4}\ Z_{\odot}$ and their initial temperature is $10^4$ K.
The dark matter halo has a NFW density profile \citep{NFW}:
\begin{equation}
 \rho_{\rm NFW}(r) = \frac{\rho_{\rm s}}{(r/r_{\rm s})(1+r/r_{\rm s})^{2}}
\end{equation}
where $\rho_{\rm s}$ and $r_{\rm s}$ are, respectively, the
characteristic density and the scale radius. In order to fix the values
of these parameters, we use the correlation between them found by
\cite{wechsler02} and \cite{gentile04}, which makes the NFW density
distribution essentially a one-parameter family of the dark matter
virial mass, $M_{\rm DM}$. The relations we use for $\rho_{\rm s}$,
$r_{\rm s}$ and the concentration parameter $c$ (=$r_{\rm max}/r_{\rm
  s}$) are~:
\begin{eqnarray}
 c &\simeq& 20 \left(\frac{M_{\rm DM}}{10^{11}M_{\sun}}\right)^{-0.13}
 \\ r_{\rm s} &\simeq& 5.7 \left(\frac{M_{\rm
     DM}}{10^{11}M_{\sun}}\right)^{0.46}~{\rm kpc} \\ \rho_{\rm s}
 &\simeq& \frac{101}{3}\frac{c^{3}}{\ln(1+c)-c/(1+c)}\rho_{\rm crit}.
\end{eqnarray}
Here, $r_{\rm max}$ is the halo's
virial radius. At $r_{\rm max}$, the DM halo is truncated and the
density drops to zero, so the entire mass $M_{\rm DM}$ is situated
inside the radius $r_{\rm max}$.

\subsection{Criteria for Star formation}
\label{subsection:SFC}
Star formation is assumed to take place in cold, dense, converging and
gravitationally unstable molecular clouds \citep{katz96}. Gas
particles that fulfill the star formation criteria (SFC) are eligible
to be turned into stars. These SFC are:
\begin{eqnarray}
 \rho_{\rm g} & \geq & \rho_{\rm SF} \\ T & \leq & T_{\rm c} = 15000 \mathrm{K}
 \\ \vec{\nabla}.\vec{v} & \leq & 0,
\end{eqnarray}
with $\rho_{\rm g}$ the gas density, $T$ its temperature and $\vec{v}$
its velocity field. $\rho_{\rm SF}$ is the density threshold for star
formation. We employ a Schmidt law \citep{schmidt59} to convert gas
particles that fulfill the SFC into stars:
\begin{equation}
\frac{\mathrm{d}\rho_{\rm s}}{\mathrm{d}t} =
-\frac{\mathrm{d}\rho_{\rm g}}{\mathrm{d}t} =
c_{\star}\frac{\rho_{\rm g}}{t_{\rm g}}, \label{cstar}
\end{equation}
with $\rho_{\rm s}$ the stellar density and $c_\star$ the
dimensionless star formation efficiency. The timescale $t_{\rm g}$ is
taken to be the dynamical time for the gas $1/\sqrt{4 \pi G \rho_{\rm
    g}}$. Here, we choose $c_{\star}=0.25$. \cite{stinson06} showed that
the influence on the mean SFR of the value of $c_{\star}$ with values in
the range of 0.05 to 1 is negligible. Lowering $c_{\star}$ reduces the
star formation efficiency as well as the amount of supernova feedback, causing
more particles to fulfill the density and temperature criteria. This
compensates for the lower value of $c_{\star}$, producing a SFR which is
roughly independent of $c_\star$.

\cite{revaz09} also investigated the influence of $c_\star$ by varying it between the values of 0.01 and 0.3. They concluded that the star formation history is mainly determined by the initial total mass with a minor influence of $c_\star$. Self-regulating models, in which star formation occurs in recurrent bursts due to the interplay between cooling and supernova feedback, were achieved for $c_{\star}\sim0.2$.
            Such models best resemble real dwarf galaxies.

\subsection{Feedback}
We consider feedback from star particles by supernova Ia (SNIa),
supernova II (SNII) and stellar winds (SW). They deliver energy and
mass to the ISM and enrich the gas. Feedback is
distributed over the gas particles in the
neighborhood of the star particle according to the SPH smoothing kernel. Each star particle represents a
single-age, single-metallicity stellar population (SSP). The stars within each SSP
are distributed according to a Salpeter initial mass function~:
 \begin{equation}
\Phi(m)\mathrm{d}m = Am^{-(1+x)}\mathrm{d}m,
\end{equation}
with $x=1.35$ and $A=0.06$. The limits for the stellar masses are
$m_{\mathrm{l}} = 0.01~\mathrm{M_{\odot}}$ and $ m_{\mathrm{u}} =
60~\mathrm{M_{\odot}}$. The energy release of a SN is set to $E_{\rm
  tot} = 10^{51}$~erg and that by a SW to $E_{\rm tot} = 10^{50}$~erg
\citep{thornton98}. The actual energy injected into the ISM is
implemented as $\epsilon_{\rm FB} \times E_{\rm tot}$, where
$\epsilon_{\rm FB}$ is a free parameter.

\subsection{Cooling}
Metallicity-dependent radiative cooling is implemented using the cooling curves
from \cite{sutherland93}. With this recipe it is possible
to cool gas to a minimum temperature of $10^{4}$~K.
We also implemented the \cite{maio07} cooling curves, making it possible for particles to cool below $10^{4}$~K.

\subsection{Production runs}

In Table \ref{Models_table}, we give an overview of the parameters
that were used to set up the models. A benefit of our code is that we can retain the same initial conditions
and easily adapt our parameters to perform a detailed parameter
survey. In the remainder, we will quantify the density threshold by $n_{\rm SF}$ expressed in hydrogen ions
per cubic centimeter (so $\rho_{\rm SF} = 1~{\rm amu} \times n_{\rm SF}$). At the start of the simulations,
the models only contain
gas and dark matter. During the first few $10^8$~years, the gas collapses
in the gravitational potential well of the DM. The simulations run for
12.22 Gyr, till $z=0$.

\begin{table}
 \caption{Details of the basic spherical dwarf galaxy models that were
   used in the simulations. Initial masses for the DM halo and gas
   are given in units of $10^{6} \mathrm{M_{\odot}}$, radii in
   kpc.\label{Models_table}} \input{Models_table}
\end{table}

In the literature, a large variety of values for the density threshold can be found. \cite{stinson06} use a low density
threshold of $0.1$~cm$^{-3}$ while \cite{governato10} use a high density threshold of $100\ \rm cm^{-3}$ which, these authors
argue, is a better representation of the conditions in star-forming regions in real galaxies. The simulations of \cite{sawala11a}
have been performed with a density threshold of $10\ \rm cm^{-3}$. In this paper, we increase the density threshold from $n_{\rm SF}=0.1$~cm$^{-3}$, over
$n_{\rm SF}=6$~cm$^{-3}$ to $n_{\rm SF}=50$~cm$^{-3}$. For the fiducial series of low-density threshold simulations, we matched the $n_{\rm SF}=0.1$~cm$^{-3}$
with a feedback efficiency of $\epsilon_{\rm FB}=0.1$. For the intermediate-density threshold simulations, with $n_{\rm SF}=6$~cm$^{-3}$, we varied the
feedback efficiency between $\epsilon_{\rm FB}=0.1$ and $0.9$. Finally, for the high-density threshold simulations, with $n_{\rm SF}=50$~cm$^{-3}$, we varied the
feedback efficiency between $\epsilon_{\rm FB}=0.3$ and $0.9$.

\section{Introduction}

Dwarf galaxies are the most common type of galaxy in the local
universe but also the faintest and least easy to observe.  In the
$\Lambda$CDM cosmology, our universe consists of matter, both
luminous and dark, and dark energy, which is responsible for the
accelerating expansion of the universe. Galaxies form when gas
collapses in dark matter halos. Baryons, be it in the form of gas, dust
or stars, are the most accessible form of matter, emitting radiation
over the whole electromagnetic spectrum. Dark matter, on the other
hand, as it only interacts gravitationally, is much more difficult to
``observe''.

There have been many attempts to estimate dark halo masses and
mass-to-light ratios for galaxies and clusters of galaxies from direct
observations. These include methods that make use of gravitational
lensing \citep{mandelbaum06, liesenborgs09},
dynamical modeling of the observed properties of a kinematical tracer
such as stars or planetary nebulae \citep{kronawitter00, derijcke06, napolitano11, barnabe09}.
One thing virtually all these works have in
common is the relatively limited size of the data set they are based
on. \cite{guo10} determined the halo mass as a function of stellar
mass for a large sample of galaxies using a statistical analysis of
the Sloan Digital Sky Survey, which yields the stellar masses, and the
Millennium Simulations, which yield the dark-matter masses. In the
range of the most massive halos and bright galaxies, the derived
$M_{\rm star}$-$M_{\rm halo}$ relation, which is of the form $M_{\rm
  star} \propto M_{\rm halo}^{0.36}$, is found to be in good
agreement with gravitational lensing data \citep{mandelbaum06}. Below
a halo mass of $M_{\rm halo} \sim 10^{11.4}$~$M_\odot$, this relation
becomes much steeper:~$M_{\rm star} \propto M_{\rm
  halo}^{3.26}$. \cite{guo10} extrapolate the latter relation into the
dwarf regime, where $M_{\rm halo} \lesssim 10^{10}$~$M_\odot$. This
leads then to the prediction that faint dwarf galaxies with stellar
masses of the order of $M_{\rm star} \sim 10^6$~$M_\odot$ should live
in comparatively  massive $M_{\rm halo} \sim 10^{10}$~$M_\odot$
dark-matter halos.

The \cite{guo10} $M_{\rm star}$-$M_{\rm halo}$ relation was compared
with that found in simulations of dwarf galaxies
\citep{valcke08,stinson07,stinson09,governato10,pelupessy04,mashchenko08}
by \cite{sawala11a} and \cite{sawala11b}. They found that simulated
dwarf galaxies had stellar masses that were at least an order of
magnitude higher at a given halo mass than predicted by
\cite{guo10}. There could be several causes for numerical dwarf
galaxies to be overly prolific star formers:
\begin{itemize}
\item The star formation efficiency could be too high because of an
  underestimation of the feedback efficiency. \cite{stinson06}
  investigated the influence of the feedback efficiency on the mean
  star formation rate (SFR). The general trend they have observed was a
  decrease of the mean SFR when increasing the feedback efficiency.
\item \cite{stinson06} also reported finding a decreasing mean SFR with
  increasing density threshold for star formation. Recently, high
  density thresholds for star formation have come in vogue, see
  e.g. \cite{governato10}.
\item Dwarf galaxies, due to their low masses, are expected to be
  particularly sensitive to reionisation. Not properly taking into
  account the effects of reionisation may lead to an overestimation of
  the gas content of dwarfs and an underestimation of the gas cooling
  time.
\item Dwarf galaxies are metal poor and hence also dust poor. This
  lowers the production of H$_2$ molecules and causes poor
  self-shielding of molecular clouds \citep{buyle06} which
  could be expected to inhibit star formation. Not taking these
  effects into account will lead to an overestimation of the
  SFR \citep{gnedin09}.
\end{itemize}

Using the high values for the density threshold above which gas
particles become eligible for star formation, denoted by $\rho_{\rm SF}$, as promoted by \cite{governato10}, in combination with
radiative cooling curves that allow the gas to cool below $10^4$~K
\citep{maio07}, makes the gas collapse into small, very dense and cool clouds
before star formation ignites. If the supernova feedback
$\epsilon_{\rm FB}$, defined as the fraction of the average energy
output of a supernova that is actually absorbed by the interstellar
medium (ISM), is too weak to sufficiently heat and/or disrupt such a
star-forming cloud, one can consequently expect the mean SFR to be
very high, leading to overly massive (in terms of $M_{\rm star}$)
dwarfs. Therefore, one could hope to remedy this situation by
increasing $\epsilon_{\rm FB}$ accordingly. In that case, a
correlation between $\epsilon_{\rm FB}$ and $\rho_{\rm SF}$ would be
expected to exist.

In the present paper, we analyze a large suite of numerical
simulations of isolated, spherically symmetric dwarf galaxies in which
we varied both the feedback efficiency $\epsilon_{\rm FB}$ and the
density threshold $\rho_{\rm SF}$. Our goal is to investigate {\em (i)}
if such a correlation between $\epsilon_{\rm FB}$ and $\rho_{\rm SF}$
exists and, if it exists, how to break it, {\em (ii)} which
$\epsilon_{\rm FB}$/$\rho_{\rm SF}$-combinations lead to viable dwarf
galaxy models in terms of the observed photometric and kinematic
scaling relations, and {\em (iii)} how well these models approximate
the aforementioned $M_{\rm star}$ -$M_{\rm halo}$ relation.

In section \ref{section:numerical_details}, we give more details about the numerical methods that
are used in our code. An analysis of the simulations is given in
section \ref{section:analysis}, where some details are given of the NFW halo that is used
for the simulations and a large set of scaling relations are
plotted comparing our models to observations. In section \ref{section:results} we discuss the obtained results and
 conclude.

\section*{Acknowledgements}
We thank Volker Springel for
making publicly available the {\sc Gadget-2} simulation code and Till Sawala for making available to us his
data from the Aquila simulation. We also would like to thank the anonymous referee for his/her stimulating remarks
that have greatly improved the manuscript.

\input{article_bibliography}

\label{lastpage}

\end{document}

\section{Discussion and conclusions}
\label{section:results}

\begin{figure}
 \centering
 \includegraphics[width=\columnwidth]{figures/Mstar_Mhalo_2C.eps}
 \caption{ The stellar mass versus the DM halo mass, plotted in
   comparison with the models by \citet{sawala11b}. The gray dots show data from gravitational lensing from \citet{mandelbaum06}. The black
   line is the trend for this relation that was determined by
   \citet{guo10}. \label{fig:MstarMhalo}}
\end{figure}

\subsection{Cusp to core}
Whether the halo density profile is cusped or cored has been a point of discussion for quite some time. Observationally, evidence for cored DM profiles is found \citep{gentile04}, but from cosmological DM simulations a cusped density profile is deduced \citep{NFW, moore96}. The inherent limitation due to the angular resolution of the observations is ruled as a cause of the observed flat density profiles by \cite{deblok02}. \cite{gentile05} also excluded the possibility of non-circular gas motions which might result in a rotation curve that is best fitted by a cored halo, while the dark matter halo actually has a cuspy profile. However, from the simulation point of view, \cite{mashchenko06} mentioned a natural transition of a cusp to a flattened core when the dark matter halo is gravitationally heated by bulk gas motions.

Our simulations are set up with a cusped NFW halo in agreement with cosmological simulations. The infall of gas causes an adiabatic compression of the dark halo. When gas is evacuated from the central
regions, be it by a fast re-expansion as the gas pressure builds up or by supernova feedback, the dark-matter halo reacts non-adiabatically and kinetic energy of the gas is transferred to the dark matter. This results in a flattening of the central density and so the cusp is converted into a core. We can conclude that the conversion of the cusped halo density profile to a cored profile is realized by the removal of baryons
from the galaxy center \citep{read05}, whether this is due to a re-expansion of the gas or by feedback effects or by another process.

\subsection{Degeneracy}
By increasing both the  density threshold and the feedback efficiency, the simulated galaxies  move along the observed kinematic and photometric scaling relations.
These two parameters, the feedback efficiency $\epsilon_{\rm FB}$ and the density threshold $n_{\rm SF}$, correlate with each other and an increase of the one can be counteracted by an increase of the other, resulting in galaxies with similar properties. To be more specific:~the individual galaxies are drastically different for different parameter values but they all line up along the same scaling relations and
can therefore be seen as good analogs of observed dwarf galaxies.

The feedback efficiency quantifies the fraction of the $10^{51}$~ergs of energy that are released during a supernova explosion and thermally injected into the ISM. For each value of the density threshold we can determine the feedback efficiency range for which the models are in agreement with the observations, although we are not able to deduce a unique $n_{\rm SF}$/$\epsilon_{\rm FB}$-combination which would be the ``correct`` representation of the physical processes that happen in galaxies.

For a certain density threshold, a lower limit of the corresponding $\epsilon_{\rm FB}$-parameter can be determined from the effective radius: the galaxies become too centrally concentrated when the feedback is too low.
From the scaling relations we cannot deduce an upper limit for the $\epsilon_{\rm FB}$-parameter, but one could argue that the ISM cannot receive more energy than there is released by the supernova explosion, resulting in a maximal value for the feedback efficiency of 1.

In the case of a density threshold of $n_{\rm SF}=0.1$~cm$^{-3}$, the models are generally in good agreement with the observations besides the somewhat high metallicities. This is also the reason why the feedback efficiency was not varied in this case. If we compare the high density threshold models, $n_{\rm SF}>0.1$~cm$^{-3}$, with the observations we can conclude that the feedback efficiency should be larger then $\sim 0.3$. For a density threshold of $n_{\rm SF}=6$~cm$^{-3}$, we prefer a value of 0.7 for the feedback. Similarly we prefer a feedback efficiency of 0.9 in the case of a density threshold of $n_{\rm SF}=50$~cm$^{-3}$

The fact that different $n_{\rm SF}$/$\epsilon_{\rm FB}$-combinations result in simulated galaxies with properties that are in agreement with the observations invokes a warning for future simulations and indicates that there is still some work left to determine the density of the star forming regions and the fraction of supernova energy that is absorbed by the ISM, quantities which are hard to determine observationally.

There are however other parameters that might influence the starformation rate and our degeneracy, which are not investigated here:
\itemize{
  \item{Given the fact that the star-formation efficiency $c_\star$ was found by other authors not to have a significant impact on stellar mass, we did not investigate it in detail in this paper.}
  \item{The choice of the IMF, for which in our simulations a Salpeter IMF is used, determines the mass distribution of stars. The fraction of high-mass stars  influences the number of SNIa and SNII explosions and as a consequence it will influence the amount of feedback and the chemical evolution. However, given
the large number of IMF parameterizations available in the literature, testing them is a very daunting task which falls outside the scope of this paper. Moreover,
part of the IMF-variation is quantified approximately by the variation in $\epsilon_{\rm FB}$ which we do investigate. }
  \item{There are other possible feedback implementations, next to the release of feedback energy as thermal energy to the gas. It also could be released as kinetic energy by kicking the gas particles or by blast-wave feedback \citep{mayer08}.}
\item Other implementations of star formation, e.g. based on a subgrid model of H$_2$-formation \citep{pelupessy04}, are possible.
}

\subsection{The dwarf galaxy dark-matter halo occupancy}
To conclude, Fig. \ref{fig:Mst_Mh_f} shows the models which best agree with the observations for each density threshold that was used in our analysis.
Increasing $n_{\rm SF}$ together with $\epsilon_{\rm FB}$ leads to a strong reduction, of almost two orders of magnitude, of the stellar mass, especially in the most massive models.
However, with the physics included in our simulations, we are unable to reproduce the $M_{\rm star}-M_{\rm halo}$ relation of \cite{guo10}. Surprisingly, the best models trace a $M_{\rm star}-M_{\rm halo}$ relation
with a slope that is similar to that of the relation of \cite{guo10}. Our simulations are in agreement with results from  cosmological simulations, which have, however, much lower spatial resolution in the dwarf regime \cite{sawala11b}. We did not explore yet higher values for $n_{\rm SF}$ and $\epsilon_{\rm FB}$ because it is clear from Fig. \ref{fig:Mst_Mh_f} that the reduction of $M_{\rm star}$ stagnates for
high $n_{\rm SF}$-values. Moreover, to compensate for the high density threshold, an unphysical large value for $\epsilon_{\rm FB}$, higher than 1, would be required. Thus, we arrive at $(n_{\rm SF}=6\,{\rm cm}^{-3},\, \epsilon_{\rm FB}\sim 0.7)$ and $(n_{\rm SF}=50\,{\rm cm}^{-3},\, \epsilon_{\rm FB}\sim 0.9)$
as the models which are in best agreement with the observed photometric and kinematical scaling relations and with the $M_{\rm star}-M_{\rm halo}$ relation derived directly from cosmological simulations.

\begin{figure}
 \centering \includegraphics[width=\columnwidth]{figures/Mstar_Mhalo_2C_final.eps}
 \caption{The $M_{\rm star}-M_{\rm halo}$ of our best models for different density threshold compared to the relation of \citealt{guo10}, other simulations from \citealt{sawala11b} and observations from \citealt{mandelbaum06}.
 \label{fig:Mst_Mh_f}}
\end{figure}

While it appears impossible to place isolated dwarf galaxies on the $M_{\rm star}-M_{\rm halo}$ relation of \citet{guo10}, it is possible to envisage external influences that may further reduce
$M_{\rm star}$, as already mentioned in the Introduction:
\begin{itemize}
 \item Not properly taking into
  account the effects of reionisation may lead to an overestimation of
  the gas content of dwarfs and an underestimation of the gas cooling
  time. However, even taking into account reionisation, the dwarf galaxies simulated by \citet{sawala11a}
  had much too high stellar masses.
\item At a given gas density, the star-formation efficiency of dwarf galaxies could be lower than that of more massive stellar systems because of their lower metallicity and hence lower dust content.
  This could be mimicked by reducing the star-formation efficiency parameter $c_\star$ (see eq. (\ref{cstar})) in the dwarf regime. However, \citet{stinson06} have shown that, because of self-regulation, the star-formation rate is very insensitive to this parameter:~varying $c_\star$ between 0.05 and 1 left the mean star-formation rate virtually unchanged.
\item External processes such as ram-pressure stripping and tidal stirring may lead to a premature cessation of star formation and hence lower $M_{\rm star}$ \citep{mayer06}. However, these processes are
only effective if the gravitational potential wells of dwarf galaxies are sufficiently shallow and if they are stripped early enough in cosmic history, before they converted their gas into stars. It is
unclear whether these constraints are met. In \citet{derijcke10}, and references therein, it was argued that the number of red-sequence, quenched dwarf galaxies increased significantly over the last half of the Hubble time and that the dwarf galaxies now residing in the Fornax cluster were accreted less than a few crossing times age (i.e. less than a few Gyr). This timescale would have left dwarf galaxies ample time to form stars before entering the cluster.
\end{itemize}


text:
\section{Introduction}\label{sec:intro}
In natural language processing (NLP), the notion and construction of an \textit{embedding} (i.e., a mapping of a linguistic object such as a word, phrase, sentence, or an entire document to a vector in Euclidean space) is the essential link for making precise the features of language that we hope to learn or understand with statistical or machine learning tools. \textit{Word embeddings}, arguably introduced in \citet{deerwester1990indexing}, have grown greatly in popularity since their utility in downstream tasks was demonstrated by \citet{collobertDownstream}. Specifically, the Word2Vec algorithm \citep{mikolovWord2Vec} greatly influenced the use of word embeddings by providing a fast and effective unsupervised approach to constructing word embeddings. This was quickly followed by the GloVe algorithm \citep{penningtonGlove}, which had performance comparable to Word2Vec. More modern deep learning word embedding algorithms, such as ELMo \citep{petersElmo}, BERT \citep{DevlinBERT}, and GPT-3 \citep{brownGPT3}, have further pushed the performance of word embeddings to state-of-the-art levels on a variety of downstream tasks such as masked language modeling, part-of-speech tagging, analogy completion, and text generation. Though these powerful techniques have already demonstrated far-reaching societal impacts, it remains important to understand the nature of the embeddings that underlie their success.

All the aforementioned embedding techniques employ some strategy to generate word embeddings from a corpus based on some intuited feature of natural language, theoretically giving the word embeddings relational meaning. Some examples of such strategies include the order of words in a sentence, how often pairs of words appear near one another, and the number of times certain words occur in various documents.
However, embeddings are {\em not} ultimately evaluated on how well they represent the intuited features of natural language that they are designed to learn---NLP algorithms are typically only judged by their success on downstream tasks. In fact, due to a lack of any precise mathematical formulation of theoretical features that govern the meaning of natural language, it is not clear that embeddings have any significance beyond auxiliary data features that are useful for training NLP models on downstream tasks (such as those tested by the GLUE benchmark of \citealp{wang2019glue}, evaluating on sentiment classification, semantic equivalence evaluation, text similarity evaluation, and recognition of textual entailment, among others).

We argue that to begin to understand the natural language features that embeddings represent, we must first mathematically formalize the precise features of natural language that we conjecture to make inferences on.  Second, we must be able to generate synthetic natural language data that exhibit the precise features that have been formulated (i.e., there must exist a generative model).  And third, we must be able to demonstrate theoretical and/or empirical statistical consistency of estimation procedures designed to learn the features.  Without any understanding of how natural language can be generated, attempting to gain insight into learned features---much less gain statistical guarantees on estimation properties for ``true" underlying linguistic features---is hopeless.

The contributions of our paper are as follows:
\begin{itemize}
\item We consider a theoretical formulation of the skip-gram algorithm used for training unsupervised word embeddings, as in Word2Vec.  It is established in \citet{levy2014neural} that the skip-gram algorithm minimizes its loss function at the pointwise mutual information (PMI), but the statistical properties of the estimated embedding features are not investigated.  We investigate the statistical properties via simulation studies by proposing a copula-based statistical model for natural language data that truly has a given PMI matrix as a feature of the generative model.
\item We provide a solution to the problem of constructing a singular value decomposition (SVD) of the PMI matrix from real natural language text data, which exhibits many infinite-valued components. \citet{levy2014neural} propose ad hoc truncation rules, but we adapt recent developments in missing values SVD (MVSVD) algorithms that rely on expectation-maximization (EM) based imputation and estimation for matrices with missing components.  Moreover, the right and left singular vectors resulting from the SVD of a PMI matrix can be used as embedding vectors with comparable meaning to the embedding vectors trained from the skip-gram algorithm and with similar performance in training on downstream tasks but with improved interpretability. This is implicit in the fact that matrix decompositions such as SVD are relatively well-understood mathematically, unlike unsupervised learning algorithms.
\item The copula-based statistical model that we propose is motivated by linguistic literature on theoretical distributional features of data from natural language, namely that ordered word frequencies in text corpora are often Zipfian distributed.  We illustrate that it is possible to build generative models for natural language data that can be used to study theoretical features not limited to the PMI matrix.
\item While the construction of embedding vectors in NLP contexts is almost exclusively driven by training models for optimized performance on downstream tasks, our work forces the question of whether embedding vectors have inherent meaningfulness for representing features of natural language.
\end{itemize}
The organization of our paper is as follows.  An overview of existing language modeling approaches is given in Section \ref{sec:prevwork}. This is followed by a brief discussion of the theoretical properties of the Word2Vec algorithm and our proposed model for natural language in Section \ref{sec:w2vgenest}. In Section \ref{sec:PMIestimation}, we discuss various estimators that can act as alternatives to skip-gram embeddings, examining previously used proposals as well as our new algorithm. Section \ref{sec:simstudy} provides simulation studies to analyze Word2Vec and its alternatives in terms of their abilities as statistical estimators. Finally, Section \ref{sec:realdata} presents the performance of these algorithms on a standard sentiment analysis task to illustrate that similar behavior as statistical estimators yields similar performance in downstream tasks. The code to reproduce the results of this paper can be found at \url{https://github.com/neil-dey/word-embeddings-as-estimators}.

\section{Existing Studies}\label{sec:prevwork}

The analysis of text data and word embeddings are largely an unexplored avenue of research in the statistics literature. There do, however, exist papers in the statistics literature on the related topics of sentiment analysis \citep{taddy2013multinomial, wang2016classification}, topic modeling \citep{airoldi2016improving, li2021topic, chen2022learning, wu2022sparse}, and text classification \citep{zhou2016priors, xia2021intentional}.

Previous studies in the statistics and linguistics literature more closely related to our work include the following.  In \citep{dai2022coupled}, an approach is proposed for sentence generation in next-word prediction, and it is shown that the generation error can be sufficiently well estimated. Recently, Dai et al. developed and proposed a concept of U-minimal sufficient embeddings, a class of embeddings that minimizes a predefined U-embedding loss. They then provide a theory determining when U-minimal sufficient embeddings can be recovered from text \citet{dai2022embedding}. Related to our direct application, work by \cite{yan2021rare} demonstrates how the sparsity of a co-occurrence matrix can be reduced by grouping rare words together in aggregate features. However, their technique assumes additional side information which encodes word similarity, namely a hierarchical clustering of GloVe embeddings. If our generative models were to reduce their complexity using such a technique, we would most likely require the side information used in comparing semantic similarity not to be generated from an arbitrary embedding. Finally, word embeddings were recently extended in \cite{feng2022word} to include information regarding grammar patterns. They additionally argued that multi-word patterns follow a Zipfian distribution, and pattern lengths are inversely related to their frequency. The commonality of both text and language patterns following Zipfian distributions allows for an easy extension of our copula text generation procedure to include identified multi-word phrases.

A great deal more work has been done beyond the scope of the statistics and linguistics literature in the area of natural language generation and language modeling.
A standard model for text generation is the $n$-gram model, in which words are generated using an $n$th order Markov chain.
These $n$-gram models are popular due to their simple nature allowing for ease of statistical analysis; however, a naive $n$-gram model, creating transition probabilities by using observed frequencies from a training corpus, often fails to capture desirable linguistic features.
To remedy this limitation, a variety of methods have been proposed. A common approach to creating more sophisticated $n$-gram models is smoothing---adjusting maximum likelihood estimates for probabilities by making the distribution of word occurrence probabilities more uniform \citep{chen1998Smoothing}.
Examples include Good-Turing estimation \citep{good1953Smoothing}, Jelinek-Mercer smoothing \citep{Jelinek1980InterpolatedEO, brownEntropy}, and Bayesian smoothing \citep{nadasSmoothing, mackayModel}, among others.
However, most attempts at smoothing still fail to take into account how similar words occur in similar contexts; class-based $n$-gram models \citep{brownClassModel} address this by separating words into distinct classes based on the frequency of co-occurrence with other words.
The primary limitation of such an approach is in the difficulties encountered when words can have disparate meanings and appear in wildly different contexts.

Other approaches more sophisticated than the $n$-gram model also exist. Neural network approaches to language modeling, such as those inspired by \citet{bengioNeuralModel}, address many of these issues and perform admirably in generating reasonably likely text (as measured in \citet{bengioNeuralModel} by the perplexity metric).
However, jumping to neural networks again decreases the interpretability of the model and makes proving theoretical results difficult. Another approach is found in the log-linear models of \citet{minhModel} and \citet{aroraPMI}. These models do provide sophisticated statistical models for language that yield interesting theoretical results but are tailored to capturing specific linguistic features and do not generalize easily to multiple features.

\section{Statistical Framework}
\label{sec:w2vgenest}
\subsection{Word2Vec and Pointwise Mutual Information}\label{sec:w2vpmi}

Given a training corpus, \citet{levy2014neural} proved that the Word2Vec algorithm, using a skip-gram algorithm with one negative sample, implicitly factors the empirical PMI matrix, given by
\begin{equation*}
    \operatorname{PMI}(w, c) = \log \frac{\Pr(w, c)}{\Pr(w)\Pr(c)}
\end{equation*}
where $w$ and $c$ together form a word-context pair, $\Pr(w, c)$ is the probability of drawing $(w, c)$ as a word-context pair from the corpus, and $\Pr(w)$ and $\Pr(c)$ are the probabilities of drawing $w$ and $c$ respectively from the corpus. That is, Word2Vec generates a matrix $W$ of word embeddings and a matrix $C$ of context embeddings, and its objective function attains a global maximum when $W$ and $C$ are such that $WC^\top = \text{PMI}$. More generally, Word2Vec with $k$ negative samples factors the empirical Shifted PMI (SPMI) matrix
\begin{equation*}
    \operatorname{SPMI} = \operatorname{PMI} - \log(k)\cdot\vb*{J}
\end{equation*}
where $\vb*{J}$ is the all-ones matrix. Note that neither the empirical PMI nor SPMI matrix is guaranteed to have only finite entries. In a finite corpus, most words do not co-occur with each other, leading to $\Pr(w, c) = 0$ for any such non-co-occurring pair and hence $\log \Pr(w, c) = -\infty$. The presence of many $-\infty$ entries in the empirical PMI matrix limits the mathematical and statistical techniques that can be used to interpret the matrix directly, as most linear algebra techniques require the entries of a matrix to come from a ring or field, but $\mathbb{R}\cup\{\pm\infty\}$ fails to even form a semigroup. Thus, a mathematical analysis of the PMI matrix and Word2Vec will require a deeper understanding of these infinite entries.

\subsection{Two Settings for Natural Language Modeling}\label{sec:nlpdistbndisc}
We identify two mutually exclusive possibilities for the nature of a language that produces infinite entries for an empirical PMI matrix: a \textit{sparse} setting and a \textit{dense} setting.
In both settings, we interpret the observed corpus as a \textit{sample} which is obtained from a hypothetical, infinite \textit{population} of text data. This population-sample framework is a fundamental building block of statistical reasoning.
In any statistical inference problem, the goal is to learn something about the population, e.g., the population mean.
However, it is impossible to actually observe the entire population, so we use the statistical sample as a representative and finite subset of the population on which calculations can be carried out in a feasible manner, e.g., computing the sample mean.
We then use the sample metric as an estimate of the unknown population metric of interest.

We now carry this intuition to text data and, more specifically, co-occurrence and PMI, under sparse and dense settings.
In practice, the \textit{sample} co-occurrence matrix, i.e., the co-occurrence matrix computed from the observed corpus, is almost always sparse, with most of its entries being zero.
This implies that the corresponding terms of the PMI matrix are $-\infty$.
However, this does not necessarily mean that the \textit{population} version of the co-occurrence matrix and the PMI matrix suffers from the same issues.

In the sparse setting, there do indeed occur words (say $w$ and $c$) that never appear in the same context in the population itself; therefore, any corpus will have an infinite entry for $(w, c)$ in its corresponding empirical PMI matrix. On the other hand, the dense setting allows for any two words to appear with each other in \textit{some} context in the population (though this context may be very rarely seen); thus, for any two words, there is a positive probability of their co-occurrence. In other words, there is a corpus that contains a finite entry in the empirical PMI matrix for the word-context pair.
Note that even under the dense setting, we are likely to observe many zero entries in the observed co-occurrence matrix since the sample might not be large enough for observed frequencies to be non-zero even though the underlying probabilities are.

Both of these settings have an underlying notion of the ``truth" of the language that Word2Vec intends to train on. However, discussing ``truth" in a statistical sense is hard, if not impossible, without any understanding of how the data itself is generated. Indeed, the lack of a data-generating model for NLP data contributes to the need for the use of performance on downstream tasks for various NLP methods. To draw an analogy, suppose that a new method for regression was proposed. Rather than measuring the efficacy of this new method using downstream tasks, the designers could either directly prove results or at least test the method using simulated data $(X_i, Y_i)_{i=1}^n$ from $Y_i = f(X_i) + \varepsilon_i$ for some function $f$, covariates $X$, and random errors $\varepsilon_i$ from some distribution. With a data-generating model for natural language, we can proceed similarly with NLP methods.

Thus, choose some data-generating model for natural language. Then in the dense setting, the sparsity (in the sense of having many infinite entries) of the empirical PMI is an artifact of sampling from the model rather than a consequence of theoretical sparsity. On the other hand, in the sparse setting, sparsity is an inherent attribute of language (i.e. the data-generating model) itself. A mathematically rigorous definition of these two settings is provided in Section \ref{sec:theoframework}.

Under a sparse model, many traditional mathematical approaches to interpreting the PMI are fruitless, since infinite entries are inherent to the data itself: There necessarily exists a pair $(w, c)$ such that $\Pr(w, c) = 0$, and thus $\log \Pr(w, c) = -\infty$ in both the population and empirical PMI matrix. Indeed, the objective function of Word2Vec itself (without normalization of word and context vectors) would require these words to have infinitely large dot products with each other, so even direct analysis of the algorithm proves difficult. It is due to these complications that the remainder of this paper primarily focuses on a dense model for language.
\subsection{A Text Generation Framework}\label{sec:theoframework}
In light of the need for a data-generating model for natural language, we provide the following theoretical framework: Let $\mathscr{C}$ be an infinite sequence of tokens given by $(t_i)_{i=-\infty}^\infty$. We assume that $\mathscr{C}$ is ``stationary" in the sense that any corpus (i.e. a finite substring from $\mathscr{C}$) we draw would form a representative sample---compare this to a non-stationary $\mathscr{C}$ that adjoins a fantasy textbook with a medical textbook, where which section of $\mathscr{C}$ we draw from to sample a corpus would matter.
This is inspired by the statistical framework of time series analysis, where the observed data is assumed to be a finite, continuous subset of an infinite temporal sequence. Every finite sequence of tokens can be assigned some probability of occurrence based on their frequency in $\mathscr{C}$; the $n$-gram model assigns the probability of the sequence of tokens $(w_i)_{i=1}^m$ to be
\begin{equation*}
    \Pr(w_1, \ldots, w_m) = \prod_{i=1}^m \Pr(w_i \mid, w_{i-(n-1)}, \ldots, w_{i-1}).
\end{equation*}
Hence, an $n$th order Markov chain can be used to generate data if the transition probabilities are specified (e.g. via a transition probability matrix or tensor); for simplicity, we restrict ourselves in this paper to a unigram (first order) model. We can thus formally define the notions of dense and sparse models for language, as discussed in Section \ref{sec:nlpdistbndisc}: If the associated Markov chain is strongly connected, then the language model is dense; otherwise, it is sparse.

It then still remains to answer the question of how the transition probabilities should be specified to match the distributional properties of natural language.  A natural distributional property that the model's transition probabilities should match is the marginal distribution for word frequencies (i.e. the probability of occurrence for individual words). It is well-established that words in the English language tend to follow a Zipfian distribution with Zipf parameter approximately $1$ \citep{sanchezZipf}. To illustrate this phenomenon, we display the empirical frequencies of words from the Brown Corpus \citep{Francis1979} and overlay the expected frequencies from a Zipfian distribution in Figure \ref{fig:zipf_fit}.
\begin{figure}[htb]
    \centering
    \includegraphics[width=0.75\textwidth]{Images/zipf_fit.png}
	\caption{The empirical word frequencies from the Brown Corpus and expected word frequencies from a $\operatorname{Zipf}(1, V)$ distribution, where $V\approx 10^4$. Infrequent words (with $<10$ occurrences) have been omitted.}
    \label{fig:zipf_fit}
\end{figure}
Thus, in a fixed vocabulary of $V$ unique words, the distribution of the word probabilities should follow a $\operatorname{Zipf}(1, V)$ distribution (though our methodology can be extended to any distribution). This thus gives us the marginal distributions for a $V \times V$ word co-occurrence matrix; however, this does not give all the information necessary to construct the co-occurrence matrix in its entirety.

The problem of constructing a bivariate distribution when only the marginals are known is solved using \textit{copula} distributions. A function ${C:[0, 1]^d\rightarrow[0, 1]}$ is a copula function if it is a joint cumulative distribution function of a $d$-dimensional random vector on the unit cube $[0,1]^{d}$ with uniform marginals. With this definition, copulas are often used in the study of multivariate data. For example, the Fr\'{e}chet-Hoeffding Theorem \citep{sklarprobmetric}
gives sharp bounds the behavior of the joint random variables, with the upper bound attained by co-monotone random variables. Another use of copulas can be found in approximating the expected value of a response function $g$ applied to a random vector $(X_1, \ldots, X_d)$ \citep{mcNeil2005}; one can simulate draws of $(X_1, \ldots, X_d)$ (whose joint distribution may be unknown) by instead drawing samples $(U_1, \ldots, U_d)$ from the (known) copula and then simply applying the inverse marginal distribution functions to each component, resulting in i.i.d. samples with the correct joint distribution. Because of the utility of copulas in studying multivariate data, they have been applied to a variety of areas, including quantitative finance \citep{li2000,ibragimovCopula}, engineering \citep{thompson2011,zhang2019}, and medicine \citep{med1,med2}, among others.

For our purpose of completing the co-occurrence matrix, the relevant result from the literature on copulas is Sklar's Theorem \citep{sklar1959}, which states that every continuous multivariate cumulative distribution function $F_{\vb*{X}}(x_1, \ldots, x_d)$ can be expressed in terms of its marginal cdfs $F_i(x_i)$ and a unique copula function $C$:
\begin{equation*}
    F_{\vb*{X}}(x_1, \ldots, x_d) = C(F_1(x_1), \ldots, F_d(x_d)).
\end{equation*}

Though Sklar's Theorem only holds for continuous distributions, many discrete multivariate distributions can still be closely approximated using copula functions. In particular, given Zipfian marginals with a Zipf parameter near $1$, a Gaussian copula with correlation matrix $R$,
\begin{equation*}
    C_R(\vb*{u}) = \Phi_R(\Phi^{-1}(u_1), \Phi^{-1}(u_2)),
\end{equation*}
 where $\Phi$ denotes a normal cumulative distribution function, does well empirically in generating a bivariate Zipfian distribution with the desired marginal distributions. This thus allows us to construct a dense co-occurrence matrix, and hence (by normalizing rows to sum to unity) a transition probability matrix to generate a corpus.

Note that though this data-generating model is limited as-is, it easily extends to other, more complex situations. The current model simply generates words via a first-order Markov chain; a straightforward extension would be to use a higher-order chain. This of course comes at a cost of having to work with tensors rather than matrices and greatly increases the computational cost of computing copulas. Another extension would be to use more sophisticated copulas than the Gaussian copula, allowing for more complex structures in a language to be replicated in the model; see \citet{ibragimovCopula} for an in-depth treatment of appropriate copulas for heavy-tailed distributions such as the Zipfian distribution. Different copulas also allow for a choice between the dense and sparse settings referenced in Section \ref{sec:nlpdistbndisc}.

\section{Embedding Algorithms to Estimate the SPMI Matrix}\label{sec:PMIestimation}

In our theoretical framework introduced in Section \ref{sec:theoframework}, a corpus is an observed finite sample from a population of infinite length. In light of this, the empirical SPMI matrix from a corpus is simply a point estimate for the SPMI matrix associated to the population. The only reason Word2Vec factors an empirical SPMI matrix rather than the population SPMI matrix is due to the limitation of having to train on a finite corpus. Though impossible in practice, Word2Vec would ideally factor the population SPMI matrix, as this would allow for the most generally applicable word embeddings. Thus, in order to theoretically understand Word2Vec, we need methods that are easier to analyze mathematically but still approximate Word2Vec in its ability to implicitly factor the population SPMI matrix as well as in its performance on downstream tasks.

\subsection{Truncating the SPMI}
As stated in section \ref{sec:w2vpmi}, the objective function of Word2Vec with $k$ negative samples aims to find matrices $W$ and $C$ such that $WC^\top = \text{SPMI}$. Hence, one can bypass the need to perform the Word2Vec algorithm by directly factorizing the empirical SPMI matrix; the simplest way to do so would be to use singular value decomposition (SVD) on the matrix
\begin{equation*}
    \operatorname{SPMI} = U\Sigma V^\top
\end{equation*}
and use $U\Sigma^{1/2}$ and $V\Sigma^{1/2}$ as our word and context embeddings respectively. However, the infinite entries of the empirical PMI matrix prevent SVD from being used directly.

To circumvent the problem of the infinite entries of the empirical SPMI matrix, \citet{levy2014neural} instead perform SVD on the empirical \textit{positive} SPMI (SPPMI) matrix
\begin{equation*}
    \operatorname{SPPMI}(w, c) = \max(\operatorname{SPMI}(w, c), 0)
\end{equation*}
which truncates all negative entries to $0$. Though the SPPMI metric is shown to perform well empirically on semantic similarity tasks, it has the significant downside of ``throwing away" a large amount of data regarding negative associations. Indeed, Table 2 of \citet{levy2014neural} illustrates that as the number of negative entries increases (e.g. via an increase in negative sampling for the Word2Vec procedure), performance of SVD on the SPPMI generally decreases. Similarly, Table 1 of \citet{levy2014neural} demonstrates that as the number of negative entries increases, SVD of the SPPMI matrix continues to deviate further from the empirical PMI. Thus, a better approach to dealing with missing co-occurrence frequencies is necessary.

\subsection{Missing Values SVD}\label{sec:mvsvd}
In the dense setting, it makes sense to look into MVSVD algorithms in order to work around the sparsity of the empirical PMI matrix. Several such algorithms are presented in \citet{Kurucz2007}. One particular algorithm, yielding an EM approach to MVSVD, is shown in Algorithm \ref{alg:EMMVSVD}. This algorithm essentially imputes missing values by minimizing the Frobenius norm between the imputed matrix and the given matrix with missing values on the non-missing values.

\begin{algorithm}[!htbp]
\caption{EM-MVSVD Algorithm from \citet{Kurucz2007}}\label{alg:EMMVSVD}
\begin{algorithmic}[!tbp]
    \Require $W$, a matrix with missing values
    \Require $d$, the number of singular values to keep
    \State $R \gets \{(i, j) \mid W_{ij} \text{ is missing}\}$
    \For {$(i, j)\in R$}
        \State $W_{ij} \gets$ initial guess
    \EndFor
    \State $U, \Sigma, V^\top \gets \operatorname{SVD}(W, d)$
    \State $\widehat{W}^{(0)} \gets U\Sigma V^\top$
    \For {$(i, j) \in R$}
        \State $W_{ij} \gets \widehat{W}^{(0)}_{ij}$
    \EndFor
    \For{$t = 1, 2, 3, \ldots$}
        \If{converged}
            \State \Return $U$, $\Sigma$, $V^\top$
        \EndIf
        \State $U, \Sigma, V^\top \gets \operatorname{SVD}(W, d)$
        \State $\widehat{W} \gets U\Sigma V^\top$
        \State $\lambda \gets \arg\,\min \sum\limits_{(i,j)\not\in R} \qty[W_{ij} - (\lambda \cdot \widehat{W}_{ij} + (1-\lambda) \cdot \widehat{W}^{(t-1)}_{ij})]^2$
        \State $\widehat{W}^{(t)} \gets \lambda \cdot \widehat{W} + (1-\lambda) \cdot \widehat{W}^{(t-1)}$
        \For{$(i, j)\in R$}
            \State $W_{ij} \gets W^{(t)}_{ij}$
        \EndFor
    \EndFor
\end{algorithmic}
\end{algorithm}
\begin{algorithm}[tbp]
\caption{DD-MVSVD Algorithm}\label{alg:DDMVSVD}
\begin{algorithmic}
    \Require $W$, a matrix with missing values
    \Require $\widetilde{W}$, a matrix approximating the ``true" matrix from a distribution
    \Require $d$, the number of singular values to keep
    \State $R \gets \{(i, j) \mid W_{ij} \text{ is missing}\}$
    \For {$(i, j)\in R$}
        \State $W_{ij} \gets$ initial guess
    \EndFor
    \State $U, \Sigma, V^\top \gets \operatorname{SVD}(W, d)$
    \State $\widehat{W}^{(0)} \gets U\Sigma V^\top$
    \For {$(i, j) \in R$}
        \State $W_{ij} \gets \widehat{W}^{(0)}_{ij}$
    \EndFor
    \For{$t = 1, 2, 3, \ldots$}
        \If{converged}
            \State \Return $U$, $\Sigma$, $V^\top$
        \EndIf
        \State $U, \Sigma, V^\top \gets \operatorname{SVD}(W, d)$
        \State $\widehat{W} \gets U\Sigma V^\top$
        \State $\lambda \gets \arg\,\min \sum\limits_{(i, j)\not\in R} \frac{\qty[\widetilde{W}_{ij} - (\lambda \cdot \widehat{W} + (1-\lambda) \cdot \widehat{W}^{(t-1)})_{ij}]^2}{\vphantom{\widetilde{\widetilde{W}}}\widetilde{W}_{ij}}$
        \State $\widehat{W}^{(t)} \gets \lambda \cdot \widehat{W} + (1-\lambda) \cdot \widehat{W}^{(t-1)}$
        \For{$(i, j)\in R$}
            \State $W_{ij} \gets W^{(t)}_{ij}$
        \EndFor
    \EndFor
\end{algorithmic}
\end{algorithm}

These methods typically excel in the case of data that is missing completely at random; this is not the case for the empirical PMI matrix, where the smallest values are the entries most likely to be missing. Indeed, this limitation is a special case of the primary downside of the EM-MVSVD algorithm: It only aims to minimize the error on the known, non-missing values, so no information regarding the missing values is incorporated. As a result, if the matrix entries come from a known distribution, the EM-MVSVD algorithm may yield a factorization that is incompatible with the generating distribution on the missing values. This thus motivates the need to incorporate distributional information into the EM-MVSVD algorithm in certain problems.

Thus, our proposed methodology to factorize the population SPMI matrix is thus as follows: Given a corpus of finite length with $V$ unique words, compute the empirical co-occurrence matrix; sort the rows and columns by marginal frequencies and then normalize the matrix to form a bivariate probability distribution. Now identify each word with its rank ($1$ through $V$, with $1$ being most frequent and $V$ being least frequent); this allows us to compute the correlation terms to use in the Gaussian copula. The copula then gives us a dense estimate for the ``true" co-occurrence probabilities, which can then be transformed into an estimate for the population SPMI matrix with no missing values.
We then use this estimate as the input $\widetilde{W}$ to Algorithm \ref{alg:DDMVSVD}, a \textit{distribution-dependent} MVSVD algorithm. This modifies EM-MVSVD (Algorithm \ref{alg:EMMVSVD}) by minimizing a Chi-Square Goodness-of-Fit statistic with the estimated population SPMI as opposed to merely trying to match the empirical SPMI on non-missing entries.

\section{Simulation Study}\label{sec:simstudy}
To study the effectiveness of these methods in estimating the population SPMI matrix, we generated text and compared the matrices generated by Word2Vec, SVD of the SPPMI, EM-MVSVD (Algorithm \ref{alg:EMMVSVD}), and DD-MVSVD (Algorithm \ref{alg:DDMVSVD}) to the population SPMI matrix.

To do so, we read in words from the Brown Corpus and sampled 500 words uniformly from the set of unique words; analogously to our proposed methodology from Section \ref{sec:mvsvd}, we then created a positive transition probability matrix to generate text as well as the associated dense SPMI matrix for this set of words via a Gaussian copula.
We then ran Word2Vec, EM-MVSVD, DD-MVSVD, and SVD on the empirical SPPMI matrix; each algorithm was asked to create 100-dimensional word embeddings based on a shift of $10$ negative samples. The factorizations produced were then multiplied back together to form an estimate of the population SPMI matrix; the root mean square errors (RMSE) of these algorithms compared to the population SPMI matrix are shown in Figure \ref{fig:sppmi_performance}.

Evidently, the SVD of the empirical SPPMI matrix is an increasingly worse approximation to the population PMI as the corpus size increases, and it yields approximations far worse than any of the other algorithms. On the other hand, the two MVSVD algorithms perform essentially identically to Word2Vec in terms of RMSE from the population SPMI matrix. This thus demonstrates, without the need to check performance on downstream tasks, that MVSVD algorithms will yield word embeddings much more similar to Word2Vec (and thus perform similarly to Word2Vec) than SVD on the empirical SPPMI matrix does.
\begin{figure}[tbp]
    \centering
    \includegraphics[width=0.75\textwidth]{Images/ppmi.png}
    \caption{The RMSE of Word2Vec, the MVSVD algorithms, and SVD on the SPPMI matrix with respect to the population SPMI matrix. The RMSEs plotted are the average over 20 independent trials. Error bars are not clearly visible, as all standard errors are $<0.5$.}
    \label{fig:sppmi_performance}
\end{figure}

\section{Real Data Analysis}\label{sec:realdata}
To provide more evidence of our claim from the end of Section \ref{sec:simstudy}, we now compare the performance of Word2Vec, the MVSVD algorithms, and SVD on the SPPMI matrix on a sentiment analysis task. We trained a simple sentiment analysis classifier on the IMDB movie review data set \citep{maas-EtAl:2011:ACL-HLT2011}, which is comprised of 50,000 anonymized movie reviews split evenly between bad reviews (1-4 star ratings) and good reviews (7-10 star ratings). The goal of the task is to train a model to predict how a review rated any given movie when given that review's text. We randomly split the data set into test-train partitions, trained each embedding technique using the training partition, then trained Bidirectional Long Short-Term Memory Networks (BiLSTMs) on movie sentiment analysis. This experiment was repeated twenty times.
The input to each two-layer BiLSTM (with 100 units) was a sequence of embeddings corresponding to the words contained in the review.  Each model's output was a binary classification layer using a binary-cross entropy loss function and a sigmoid activation function. To keep computations feasible, we removed stopwords as well as words with fewer than 300 occurrences, reducing the corpus to 3104 distinct words. Additionally, we zero-padded or truncated all reviews to 500 words. The use of truncation avoided the problem of vanishing or exploding gradients in long BiLSTM chains without severely compromising the data, as 95\% of reviews required no truncation.

\begin{table}[tbp]
    \footnotesize
    \centering
    \begin{tabular}{r c | c c c c c c c c c c}
        \hline
         & &  \multicolumn{2}{c}{One-Hot} & \multicolumn{2}{c}{Word2Vec} & \multicolumn{2}{c}{SPPMI} & \multicolumn{2}{c}{EM-MVSVD} & \multicolumn{2}{c}{DD-MVSVD} \\
         & & Acc. & S.E. & Acc. & S.E. & Acc. & S.E. & Acc. & S.E. & Acc. & S.E. \\
         \hline
         \multirow{3}{*}{\parbox{1.07cm}{Negative Samples}}
         & 1 & 0.86 & 0.0047 & 0.79 & 0.0072 & 0.82 & 0.0087 & 0.82 & 0.0070 & 0.82 & 0.0146 \\
         & 2 & 0.86 & 0.0065 & 0.80 & 0.0138 & 0.78 & 0.0103 & 0.81 & 0.0096 & 0.84 & 0.0034 \\
         & 5 & 0.86 & 0.0028 & 0.79 & 0.0173 & 0.70 & 0.0133 & 0.80 & 0.0109 & 0.82 & 0.0114\\
         \hline
    \end{tabular}
    \caption{Model accuracies in positive/negative sentiment analysis on the IMDB data ste across multiple levels of negative sampling.}
    \label{tab:lstm_results}
\end{table}

Table \ref{tab:lstm_results} shows the performance of each model across five distinct embedding algorithms and negative sampling levels of one, two, and five. In addition to the previously discussed algorithms, we included a one-hot encoding of the input words followed by a dense layer with 100 nodes, which serves the purpose of acting as a benchmark for roughly how well the ``perfect" embedding overfits to the data can do in this task. All algorithms produced 100-dimensional embeddings.

As expected, the one-hot embedding performed the best, achieving 86\% accuracy regardless of the amount of negative sampling. We see that the MVSVD algorithms perform comparably to Word2Vec across all levels of negative sampling, with DD-MVSVD performing the best among the three. We see that these algorithms still perform decently in comparison to the one-hot embeddings, achieving accuracies of roughly 80\%. This is in stark contrast to the performance of SVD on the SPPMI, which matches the MVSVD algorithms at one negative sample, but quickly decreases in accuracy to 70\% as the number of negative samples increases. This phenomenon occurs because the sparsity of the SPPMI matrix increases rapidly as the number of negative samples increases, so the embeddings contain less and less information to allow the BiLSTM to make a well-tuned classifier.

\section{Conclusions}
Through the use of our data-generating model for natural language, we find that the MVSVD algorithms perform similarly to Word2Vec in estimating the population PMI matrix, whereas SVD on the SPPMI matrix does not.
This is also reflected in the algorithms' performance on the downstream task of sentiment analysis, and perhaps to others as well.
As such, the MVSVD algorithms can be seen to be quite reasonable approximations to Word2Vec that still remain tractable to analyze mathematically.

In the future, we hope to develop a theory for our proposed DD-MVSVD algorithm, such as proving the convergence properties of the algorithm. A similar EM-type algorithm for matrix completion is studied in \citet{mazumder2010}, and theoretical convergence results are provided; we believe that following a similar path as in \citet{mazumder2010}---though extending their proofs to our DD-MVSVD algorithm---is likely nontrivial due to the additional SVD steps in our algorithm.

We also hope to be able to determine what other population parameters of natural language, if any, are correlated with performance on various downstream tasks of interest. For example, word embeddings generated by GloVe \citep{penningtonGlove} target a log-count matrix in a similar way that Word2Vec targets a PMI matrix; thus, this log-count matrix may also be considered as a population parameter for a natural language. There are, of course, other parameters derived from languages that are not based on frequencies of words: capturing grammatical structures such as adjective ordering or dominant word order, as well as measures based on semantic similarity, would also be of great interest. As part of this investigation, we believe that examining constructed languages prior to natural languages would be a reasonable course of action; whereas structures that may yield interesting population parameters in natural language are difficult both to identify and to control for, it is easier to account for such problems in simpler, constructed languages. Even exceedingly simple languages with less than a dozen words and simple grammatical structures may be able to offer insight into such parameters. More complex constructed languages exist, such as Toki Pona \citep[][with only 14 phenomes and 120 words, but a well-developed grammatical structure allowing for full communication]{tokipona}, which may act as a bridge between the study of completely artificial constructed languages and natural languages.

We also plan to investigate other more complex NLP methods and build upon our generating model to allow for more expressive modeling of language. For example, a simple extension of our basic model would be to use a more general $n$-gram model for some $n > 1$ rather than limiting ourselves to the unigram Markov chain. Additionally, we would like to investigate how other choices of copulas other than the Gaussian copula impact the estimation of various population parameters. We are also interested in studying more modern, non-static embedding techniques, such as ELMo \citep{petersElmo} and BERT \citep{DevlinBERT}, and investigating whether or not these can be interpreted as statistical estimators.

\subsubsection*{Acknowledgements}
Jonathan Williams's and Neil Dey's work in this publication was partially supported by the National Heart, Lung, And Blood Institute of the National Institutes of Health under Award Number R56HL155373.
Srijan Sengupta's work was partially supported by an NIH R01 grant 1R01LM013309.
Matthew Singer's work was partially supported by a grant from the Pacific Northwest National Laboratory Distinguished Graduate Research Program.
The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health or the Pacific Northwest National Laboratory.

redpajama_set_name: RedPajamaArXiv