\section{Introduction}
Hyperpersonalization and nudging applications are frequently laced with problems caused by the non-deterministic nature of human behavior.
Self-report data in particular often runs the risk of producing unclear or inaccurate results, as an individual's perception of their own behaviors may vary drastically from one day to the next due to an infinite number of unknown external factors.
As a result, it becomes difficult to decipher whether stochasticity in self-report results is due to the input features themselves, or some sort of unknown, outside influences that caused a participant to recall a distorted view of their behaviors.
Being able to account for and adapt to this uncertainty is crucial for the development of resilient personalization approaches in a number of real-life applications, see e.g.~\cite{mills2022personalized}.

This paper addresses the issue of modeling human behavior based on observations collected from survey questionnaires, which is a particularly challenging task due to the fundamental stochasticity of the given data.
In addition to the non-deterministic nature of human behavior, this data collection via surveys often introduces additional noise to the recorded responses.

Learning stochastic processes with conventional data approximation methods~--- e.g. statistical inference, dictionary decomposition, machine learning, etc.~--- is known to be a highly challenging task, see e.g.~\cite{angluin1988learning, natarajan2013learning, raychev2016learning, aswani2018inverse}.
Typically in cases where the stochasticity is sufficiently low, the given data is treated as a perturbation of the underlying ``ground truth.''
While such an approach is widely utilized in some practical applications (such as image segmentation, face recognition, asset pricing, etc.) it is not a universal solution, as it fails to generalize to all real-world scenarios.

In particular, the issue of human feedback is not handled well by such methods since ``human-level stochasticity'' generally cannot be viewed as a slight perturbation of the underlying manifold.
For instance, consider a user leaving feedback (like/dislike) on a piece of media content.
The outcome of this interaction depends on a number of external factors that affect the mood of the user, which in turn affect the user's rating.
Thus, fully predicting the feedback of a user is generally an unfeasible task since the same user might provide different ratings to the same piece of content under the influence of different external factors.

This issue particularly plagues machine learning algorithms, as they are often designed to treat the training data as ground truth.
Indeed, a typical machine learning algorithm assumes that human behavior, or at least the feedback received from human behavior, is deterministic.
Even though there has been a lot of development in robust machine learning (see e.g.~\cite{song2022learning} and the references therein), when faced with fundamentally stochastic data, most conventional approaches perform unsatisfactorily~\cite{erev2017anomalies}.

In this paper we propose a method of modeling human behavior data by using statistical modeling and further actualize the predictions by implementing a domain expert knowledge in the form of behavioral science research that is relevant to the setting at hand.

\subsection{Behavioral Science and Behavior Modeling}
Behavioral science is a multidisciplinary field that leverages psychology, sociology, economics, and other various disciplines to understand human behavior and how individuals make decisions in the real world.
It functions as a diverse field of study with numerous applications in areas such as public health, finance, and marketing.
Behavioral research serves as a method to examine and understand a multitude of human behaviors and the different factors that influence our decisions, attitudes, and beliefs.
Behavioral scientists conduct research and run experiments to better understand human actions and why people do the things they do.
In the social sciences, meta-analyses are becoming increasingly popular as methods for examining data from a number of independent studies of the same subject in order to determine overall trends.
Meta-analyses involve statistical analysis that combines the results of multiple scientific studies.
The accumulation of results across studies is the only known solution to the problem of sampling error inherent in small sample studies~\cite{hunter1991meta}.
Advantages of meta-analysis include, among others, illuminating trends across independent studies, maintaining statistical significance, minimizing wasted data, and finding moderator variables~\cite{rosenthal2001meta}.
Further improving meta-analytic methods by improving the transparency and reproducibility of meta-analyses is becoming a priority~\cite{lakens2017examining}.
Since 2018, more than 32,500 meta-analyses in behavioral science have been made available (Google Scholar search of “meta-analysis” + “behavior” + “choice”).

Our main goal in building a stochastic model of human feedback is to evaluate the extent to which identical demographic inputs can result in different outcomes in food preferences, due to the stochasticity and lack of predictability inherent in each of the demographics (as identified by the behavioral science literature).
Variations in human preferences are due to individual factors (e.g., age, self-efficacy), external factors (e.g., context, day), engagement factors (e.g., exposures over time), and more.
Every single factor plays a role in the stochasticity of human preferences, but we cannot know what all the factors are nor the size of the role(s) they play.
In this work we model human feedback as food preferences, with demographic variables serving as the sources of stochasticity.

\subsection{Related Work}
The focus of this paper is to enhance the conventional statistical learning approach by employing domain expert knowledge in the form of behavioral science research.
Such an approach is regularly considered in relevant literature, see e.g.~\cite{peysakhovich2017using, turgeon2020tutorial, hagen2020can, mac2020artificial}.

We are addressing the task of modeling human behavior.
Such a task is considered in Choice Prediction Competition\footnote{\url{https://cpc-18.com/}} and its follow-up literature, e.g.~\cite{erev2010choice, plonsky2018and, bourgin2019cognitive}.

We are using NHANES dataset to learn food preferences from the given demographic data.
This dataset is regularly used in behavioral literature, see e.g.~\cite{zadshir2005prevalence, block2004foods, mozumdar2011persistent, storz2021diet, storz2022shiftwork}.

We leverage behavioral science studies to construct a non-deterministic model of human behavior.
A similar approach is taken in~\cite{capote2012stochastic}, where the authors propose a non-deterministic model of human behavior during train evacuation processes; however, they do not take into account any of the individuals' information, and thus their model predicts an average behavior rather than that of a given person.

In particular, we use the Theory of Planned Behavior (TPB) to model non-determinicity of human decisions.
Similar approach was utilized in~\cite{dunn2011determinants, ajzen2015consumer}, where the authors employed the TPB in the context of food preferences.

\section{Problem Formulation}
We construct a non-deterministic model to approximate human decision-making on an example of food choice selection.
Specifically, we utilize NHANES dataset, examined in Section~\ref{sec:nhanes}, to model the participants' food choices based on their demographic information.
We use the dataset to train a statistical model as stated in Section~\ref{sec:data_modeling} and employ the behavioral insights outlined in Section~\ref{sec:besci_insights} and summarized in Table~\ref{tab:besci_metaanalyses}.

\subsection{NHANES Dataset}\label{sec:nhanes}
In this paper we utilize data from the Centers for Disease Control and Prevention’s National Health and Nutrition Examination Survey (NHANES) repository~\cite{nhanes}.
We extract participants' demographic features and food preferences from NHANES 2017--2018 survey\footnote{\url{https://wwwn.cdc.gov/nchs/nhanes/continuousnhanes/default.aspx?BeginYear=2017}} to construct a non-deterministic model that predicts the participants' food choice selections based on their demographic features.
NHANES is a large, biennial, stratified, multistage survey conducted by the Centers for Disease Control and Prevention and is one of the largest and most important cross-sectional studies conducted in the United States in terms of participant size, scope, ethical diversity, and free data accessibility~\cite{leroux2019organizing}.

For this study, we use data from various modules within the NHANES, including demographic and questionnaire data.
Specifically, we use the Diet Behavior and Nutrition questionnaire (DBQ) that asks participants about their eating styles over the past week, available for download at~\url{https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/P_DBQ.XPT} with a detailed description of the data collection at~\url{https://wwwn.cdc.gov/nchs/data/nhanes/2017-2018/questionnaires/DBQ_J.pdf}.
The data contains 29 demographic and 46 response features for the 15,560 participants.
Demographic data includes gender, age, marital status, race/ethnicity, education level, and household income.
Each participant is asked how many times they ate out (fast-food, restaurants, vending machines, etc.) vs how many times they ate at home, out of total $21$ meals (assuming breakfast, lunch, and dinner every day).

\subsection{Behavioral Expert Knowledge}\label{sec:besci_insights}
The Theory of Planned Behavior (TPB) has been the behavioral science theory of study in at least 56 meta-analyses (Google Scholar search of published article titles including all the words “meta-analysis” + “Theory of Planned Behavior”). TPB, an extension of the Theory of Reasoned Action~\cite{ajzen1985intentions}, is a psychological theory developed by Icek Ajzen~\cite{fishbein1977belief} as a model to predict and understand an individual's behavior.
The TPB argues that an individual's behavioral intention is collectively shaped and influenced by three main components: an individual’s attitude, societies’ subjective norms, and their perceived control of the behavior~\cite{ajzen1985intentions}.
TPB has been leveraged to successfully predict and explain a variety of health-related behaviors from physical activity to smoking~\cite{wing2009predicting, norman1999theory}.
Food choice is a complex human behavior that is influenced by a wide range of individual and environmental factors.

For the purposes of this paper, we choose to leverage the TPB as this theory offers a strong theoretical basis for trying to understand the role and influence of certain factors on health-related behaviors like food choice.
Namely, we utilize the available behavioral research (listed below) for the demographic data used in our study (gender, age, marital status, race/ethnicity, education level, and household income).

\subsubsection{Gender}
Research has highlighted the relationship between gender and dietary choices.
For example, in a study that examined gender differences in fruit and vegetable consumption through the lens of the TPB, it was found that women were more likely to consume fruits and vegetables when compared to men~\cite{emanuel2012theory}.

\subsubsection{Age}
Age predicts dietary preferences and choices.
For example, teenagers and younger adults are less likely to consume health promoting foods~\cite{mcdermott2015theory}.
A systematic review and meta-analysis on the association between variables specified by the TPB (i.e., Attitude, Social Norms, Perceived Behavioral Control) and discrete food choice behaviors showed stronger pooled mean effect sizes of the association between intentions and behaviors as age went up.
Being in the age group 18-29 was associated with 1) less frequent deviation from a choice (i.e., stronger intention-behavior association) and 2) when there was deviation, less of it (i.e., lower heterogeneity).

\subsubsection{Marital Status}
Some research has explored the role of marital status in dietary behaviors and food choices~\cite{konttinen2021sociodemographic, lee2013understanding}.
In a study that was conducted to understand the influence of sociodemographic factors on an individual’s healthy eating behaviors in restaurants using both the TPB and the Health Belief Model, it was shown that singles were more willing than married couples to patronize a restaurant that offered healthy food choices~\cite{lee2013understanding}.
Marital status comprised the following categories: married or living with a partner, widowed/divorced/separated, and never married.

\subsubsection{Race/Ethnicity}
Increasing evidence indicates that race plays a critical role in diet-related disparities, with racial and ethnic minority groups (i.e. African American, Hispanic, American Indian/Alaska natives) typically having poorer nutrient profiles and dietary behaviors/patterns compared to the majority population groups (e.g., white/European Americans~\cite{satia2009diet}).
Explanations for this phenomenon can be attributed to social disadvantage, and it is well documented that minority racial groups tend to have higher levels of unemployment and lower income~\cite{barr2014health}.

\subsubsection{Education}
Prior research has explored the relationship between educational attainment and diet quality, showing individuals with higher educational attainment adhering to healthier diets~\cite{backholer2016association, darmon2015contribution}.

\subsubsection{Household Income}
Nutrient-rich food tends to cost more per calorie compared to food with lower nutritional values~\cite{darmon2015contribution}.
Individuals from high socioeconomic status (SES) groups tend to consume healthier and more expensive diets and individuals from low SES groups tend to select cheaper, more energy dense diets lacking in fruits and vegetables.
Lower SES is associated with less disposable income, thus creating a barrier to adhering to a healthy diet~\cite{darmon2015contribution, bertoni2011multilevel}.

\begin{table*}
    \centering
    \caption{Demographic Determinants of Choice Uncertainty}
    \begin{tabular}{|c|c|c|c|}
        \hline
        \textbf{Demographic} & \multicolumn{3}{|c|}{\textbf{Dichotomized Choice Uncertainty}}
        \\\cline{2-4}
        \textbf{Determinants} & \textbf{\textit{Direction}} & \textbf{\textit{Strength}} & \textbf{\textit{Meta-analyses}}
        \\\hline
        Gender & \multirow{2}{*}{Positive} & \multirow{2}{*}{Very small} & \multirow{2}{*}{\cite{mcdermott2015theory, mcdermott2015theory2, riebl2015systematic, steinmetz2016effective}}
        \\(male vs. female) &&&
        \\\hline
        Age & \multirow{2}{*}{Negative} & \multirow{2}{*}{Small} & \multirow{2}{*}{\cite{mcdermott2015theory, mcdermott2015theory2, mceachan2011prospective}}
		\\(high vs. low) &&&
        \\\hline
        Marital status & \multirow{2}{*}{Positive} & \multirow{2}{*}{Small} & \multirow{2}{*}{\cite{lee2013understanding}$^{\mathrm{a}}$}
        \\(married vs. single) &&&
        \\\hline
        Race/ethnicity & \multirow{2}{*}{None} & \multirow{2}{*}{None} & \multirow{2}{*}{\cite{li2019socioeconomic}}
        \\(majority vs. minority) &&&
        \\\hline
        Education & \multirow{2}{*}{Positive} & \multirow{2}{*}{Very small} & \multirow{2}{*}{\cite{steinmetz2016effective, li2019socioeconomic}}
        \\(high vs. low) &&&
        \\\hline
        Household income & \multirow{2}{*}{None} & \multirow{2}{*}{None} & \multirow{2}{*}{\cite{li2019socioeconomic}}
        \\(high vs. low) &&&
        \\\hline
        \multicolumn{4}{l}{$^{\mathrm{a}}$Not a meta-analysis. None were found that explored the effect of marital status.}
    \end{tabular}
    \label{tab:besci_metaanalyses}
\end{table*}

We are interested in modeling choice uncertainty as it relates to food choices, where uncertainty is how often does someone deviates from a choice. See Table~\ref{tab:besci_metaanalyses} for a summary of demographic determinants of choice uncertainty.
In the provided table, the ``Direction'' refers to whether a positive or negative association was found between demographic variables and uncertainty in food choices.
``Strength'' refers to the magnitude of the positive or negative association found between demographic variables and uncertainty in food choices.
And ``None'' indicates no moderating effects of demographics on healthy food choices were found, or that the literature was sufficiently mixed to warrant an absence of definitive direction and strength.

\section{Modeling Food Choices}\label{sec:data_modeling}
We begin this section by motivating the need for a behavioral data simulator.
If someone is deciding whether or not to eat at home, then there is some probability that at that particular instant in time they might eat at home.
Denote this probability as $p\in[0,1]$.
Let $X\sim Bern(p)$ (i.e., $\mathbb{P}(X=1)=p$ and $\mathbb{P}(X=0)=1-p$) which represents the person's decision of where to eat; that is, if $X=0$ then the person will eat at home and if $X=1$ then the person eats elsewhere.

In our dataset, we have the number of times that each person ate somewhere other than home over 7 days (the corresponding probability mass function will be referred to as the ``population PMF'').
We will assume that everyone made 3 decisions on where to eat each day, for a total of 21 decisions.
Let $x_{n,i}$ represent person $n$'s decision for the $i$-th meal.
So, $x_{n,i}=1$ means the person did not eat at home and $x_{n,i}=0$ means the person are at home or did not eat meal $i$.
These are samples from $X_{n,i}\sim Bern(p_{n,i})$, where $p_{n,i}$ is the likelihood of person $n$ eating out for meal $i$.
Let $Y_n$ represent the number of times person $n$ ate out over these 21 meals.
This means $Y_n=\sum_{i=1}^{21}X_{n,i}$.

Modeling the number of times each person will eat out boils down to modeling each of their meal decisions.
In the dataset, we only have one observation for each person, which makes it impossible to accurately model the probabilities $p_{n,i}$ for each person.
Instead we will model the ``average'' person and then we will modify the probabilities using either the data or expert knowledge.
Modifications can be done as follows: if we estimate the $X_i\sim Bern(p_i)$ for meal $i$, we can regress $p_i$ to or away from 0.5 if there is more or less uncertainty in this person's choices based on their demographics, respectively.

In order to simulate the data, we begin by identifying $M$ peaks and plateau-like events in the population PMF (denoted by $h$), call them $\{k_1,...,k_M\}$.
In particular, $M$ is chosen to allow the simulated distribution to approximate the population distribution.
We will create $M$ Binomial random variables, $X_1,...,X_M$, where $X_i\sim Binom(n_i,p_i)$ which is the sum of $n_i$ i.i.d. $Bern(p_i)$ random variables.
Next, we specify the variance of $X_i$ for each $i$, denoted $\sigma_i^2$.
Recall that if $X\sim Binom(n,p)$, then $E(X)=np$, $Var(X)=np(1-p)$, and the maximum of the PMF is at $np$.
We use these facts to specify that $n_ip_i=k_i$ and $n_ip_i(1-p_i)=\sigma_i^2$.
This specifies
\[
    p_i=1-\frac{\sigma_i^2}{k_i}
    \text{ and }
    n_i=\frac{k_i}{p_i}.
\]
Finally, we want to add together the PMFs of $X_1,...,X_M$ in a way so that the sum is also a PMF.
We can use linear regression to solve the following system, which is almost certainly over-specified:
\begin{align*}
    \sum_{j}w_j\binom{n_j}{k_i}p_j^{k_i}(1-p_j)^{n_j-k_i}&=h(k_i)
    \quad\forall i
    \\
    \sum_{i}w_i&=1
\end{align*}
Let $Y=\sum_{i}w_iX_i$, then $Y$ is a random variable that approximates the population's distribution of counts.

We can modify the probabilities of $X_i$ based on the uncertainty of someone's decision using expert knowledge and, in turn, this modifies the distribution of $Y$.
Let $\alpha\in (0,1)$ represent the uncertainty change.
For each $i$, $\hat{p}_i = (1\pm\alpha)\cdot(p_i-0.5)+0.5$ where a ``$+$'' indicates less uncertainty and a ``$-$'' indicates more uncertainty (obviously constraining $\hat{p}_i$ to be in $[0,1]$).
For example, according to Table~\ref{tab:besci_metaanalyses}, males have more uncertainty in decisions than females.
Let $\hat{p}^{m}_i = (1-\alpha)\cdot(p_i-0.5)+0.5$, $\hat{p}^{f}_i = (1+\alpha)\cdot(p_i-0.5)+0.5$, $X^m_i\sim Binom(n_i,\hat{p}^{m}_i)$ and $X^f_i\sim Binom(n_i,\hat{p}^{f}_i)$, where ``$m$'' indicates male and ``$f$'' indicates female.
Then $Y^{m}=\sum_{i}w_iX^m_i$ and $Y^{f}=\sum_{i}w_iX^f_i$ can be used to represent the counts for males and females, respectively.
Furthermore, we can use $X^{\cdot}_i$ or the corresponding Bernoulli random variables to simulate actual decisions.

\section{Experiment Setup}
In this section we present our non-deterministic simulator that predicts a person's food preferences from their demographic features.
\footnote{All relevant code can be found at \url{https://github.com/acstarnes/wain2022_food_choice_modeling}.}
The demographic information is given as the participant's gender, age, marital status, race/ethnicity, education, and income.
Food preferences are represented by the number of meals the participant ate out versus the number of meals the participant ate at home.

We first extract the relevant data from NHANES dataset, as outlined in Section~\ref{sec:nhanes} and preprocess the obtained data in the following way:
\begin{itemize}
	\item Continuous demographic features (age, household income) are normalized to be in $[0,1]$;
	\item Categorical demographic features (gender, race/ethnicity, marital status) are parameterized via one-hot encoding;
	\item Ordinal demographic features (education) are assigned to equidistant points in $[0,1]$;
    \item The food preference vector is normalized to form a probability distribution over the number of times chosen to eat out.
\end{itemize}
We drop any missing values from the original dataset and obtain a set of $7,784$ datapoints.
Training and test sets are created with a random 67\%-33\% split, respectively.
We then train a statistical model via the procedure stated in Section~\ref{sec:data_modeling} and further enhance the prediction results with the behavioral science insights discussed in Section~\ref{sec:besci_insights}.

We begin by plotting the probability distributions from the training and test sets in Figure~\ref{fig:simulation}, where we identify key points of support in the training distribution at $k_0=0$, $k_1=2$, $k_2=4$, $k_3=5$, $k_4=7$, $k_5=10$, $k_6=14$, and $k_7=21$.
We set the initial variances to 1 and modify them to more accurately match the training distribution.
Specifically, we use the following $\sigma_0^2=0.2$, $\sigma_1^2=1$, $\sigma_2^2=0.8$, $\sigma_3^2=0.1$, $\sigma_4^2=0.4$, $\sigma_5^2=0.1$, $\sigma_6^2=0.1$, $\sigma_7^2=0.1$.
Once we have the $k$'s and $\sigma^2$'s, we model the training distribution as in Section~\ref{sec:data_modeling}, which can be seen in Figure~\ref{fig:simulation}.
We use histogram intersection (HI) as our measure of accuracy and find that the training HI is $93.8\%$ and the test HI is $92.1\%$ (recall that HI values are between 0 and 1, with 1 indicating the distributions are the same).
The underlying binomial distributions, centered around the $k$ values, are shown in Figure~\ref{fig:binomial_distributions}.

We now will modify the binomial probabilities based on Table~\ref{tab:besci_metaanalyses}.
For gender, there is more uncertainty with males than with females.
As such, we need to move all of the binomial probabilities closer to 0.5 for males and further from 0.5 for females.
Using a scaling factor of 0.1, the resulting probabilities for males are $p_0^m=0.23$, $p_1^m=0.5$, $p_2^m=0.77$, $p_3^m=0.878$, $p_4^m=0.899$, $p_5^m=0.941$, $p_6^m=0.944$, $p_7^m=0.946$ and for females are $p_0^f=0.17$, $p_1^f=0.5$, $p_2^f=0.83$, $p_3^f=0.962$, $p_4^f=0.987$, $p_5^f=1$, $p_6^f=1$, $p_7^f=1$.
Plotting the new distributions for male and female along with the distributions from the testing set can be seen in Figure~\ref{fig:gender}.
Additionally, the test accuracy for female is $87.4\%$ and for male is $91.0\%$.
We repeat this for marital status using a scaling factor of 0.15, the resulting distributions can be seen in Figure~\ref{fig:marital_status}.
Here, the test error when marital status is married is $86.3\%$ and when marital status is single is $91.2\%$.
The modifications made to the probabilities based on expert knowledge follow the true test distributions with reasonable accuracy.

Since the simulated distributions based on expert knowledge roughly follow the test distributions that the simulator was never trained on, data sampled from these distributions are reasonable approximations to the real data. We can use them to generate an unlimited number of new data points. One way to do this is simply by sampling from the distributions in order to obtain samples for gender and number of times eaten out. Alternatively, we can make 21 samples from the underlying binomial distributions in order to simulate individual meal choices and see the impact that gender has on them. Finally, we can simulate data based on features that are not available in our dataset by using expert knowledge and modifying the distributions as we did above with gender.

\begin{figure}
    \centering
    \includegraphics[width=0.7\linewidth]{base_simulation.png}
    \caption{Probability mass distributions based on counts in train and test sets along with the simulated distribution}
    \label{fig:simulation}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.8\linewidth]{underlying_binomials.png}
    \caption{Binomial distributions centered around $k$ values with widths based on $\sigma^2$ values}
    \label{fig:binomial_distributions}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.86\linewidth]{gender_male_modification_two_plots.png}
    \caption{Modified distributions for males and females along with the corresponding distributions for the test sets}
    \label{fig:gender}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.86\linewidth]{marital_status_Single_modification_two_plots.png}
    \caption{Modified distributions for marital status along with the corresponding distributions for the test sets}
    \label{fig:marital_status}
\end{figure}

\section{Discussion}
In this paper we construct a non-deterministic model predicting human food preferences based on demographic information.
Our model is based on the combination of real data obtained from the NHANES dataset and the domain expert knowledge obtained from behavioral science studies.

Our model can be used as a behavioral data simulator to generate an arbitrary amount of food choice probabilities for given demographic features, either real or synthesized.
Such a simulator is of interest to machine learning practitioners who work with conventional approaches and thus require vast amounts of data unavailable in behavioral surveys.

While the current work is rather preliminary, it suggests a number of immediate promising directions by utilizing the data obtained from our simulator in machine learning applications dealing with behavioral data.
One area of particular interest in this direction is the training of reinforcement learning agents on feedback obtained from agent-human interactions, which, for instance, includes content recommendation algorithms.

\subsection*{Limitations}
Dichotomizing determinants of food uncertainty into low vs. high, positive vs. negative, etc. undermines the vast amount of behavioral science literature that repeatedly demonstrates the individual variation in human behavior.

The irony is not lost on these authors that our stochastic model, although a step up from previous deterministic models, is a general model that averages the wide variety of human behavior into a simplified reduction, all in the effort to train a reinforcement learning agent how to personalize.
Over time and with additional research, our goal is to increase the complexity of our stochastic model of human behavior in a way that better appreciates each individual’s variation in preferences.

It would not be feasible nor possible to hold all determinants constant, so even if we controlled all that we could humanly control, each individual will possibly make a different choice in the same context under the same circumstances every time.

\bibliographystyle{abbrv}


\section{Introduction}

Encoding the statistics of a given dataset and generating data points that match those statistics are crucial issues in machine learning.
Many different strategies  have been devised.
Among the most popular solutions are Variational Autoencoders (VAEs)~\cite{kingma2018glow}, Generative Adversarial Networks (GANs)~\cite{goodfellow2020generative}, Energy Based Models (EBMs)~\cite{lecun2006tutorial} or their close relatives, Diffusion Models~\cite{sohl2015deep,ho2020denoising}.
The latter
are becoming increasingly popular nowadays due to their amazing ability to generate photorealistic images~\cite{dhariwal2021diffusion}.

EBMs offer several fundamental advantages over their competitors due to their simplicity:  A single neural network is involved in training, which means that fewer parameters need to be learned and training is less costly.
They are also particularly appealing for interpretive applications: Once trained, the energy function can be analyzed with statistical mechanics tools~\cite{decelle2021restricted}, or shallow EBMs can serve as an effective model to ``learn" something from the data. EBMs have been exploited for instance to infer the three dimensional structure~\cite{morcos2011direct} or building blocks~\cite{tubiana2019learning} of proteins, to generate artificial pieces of genome~\cite{10.1371/journal.pgen.1009303}, for neuroimaging~\cite{HJELM2014245}, simulation of complex wavefunctions in quantum many-body physics~\cite{doi:10.1126/science.aag2302,melko2019restricted}, or to impute missing input~\cite{stairs}.
Furthermore, EBMs naturally implement compositionability,
i.e. the possibility of combining information from different trained models to generate new composite data~\cite{hinton1999product,du2020compositional}.

EBMs have been around for a very long time, starting with the well-known single-layer~\cite{ackley1985learning}, Restricted ~\cite{smolensky1986information} or the Deep~\cite{salakhutdinov2009} Boltzmann machines, to modern applications where EBMs are encoded in Deep Convolutional Neural Networks (ConvNets)~\cite{nijkamp2019learning,Nijkamp_Hill_Han_Zhu_Wu_2020,du2019implicit,du2020compositional}, which show very similar performance in image generation as state-of-the-art GAN methods~\cite{du2019implicit}.
Despite their many advantages, the direct training of EBMs is hard and unstable in practice, a problem directly related to the approximations made to estimate the  negative term of the gradient~\cite{decelle2021equilibrium}.
As it is intractable, it must be approximated by sampling methods, typically Markov Chain Monte Carlo (MCMC) or Langevin dynamics, and such sampling is required each time parameters are updated.
However, practical applications use a number of iterations far shorter than the mixing time of the algorithm, required to ensure an ergodic exploration of the model phase space and a reliable estimation of the log-likelihood gradient.
To circumvent this ``equilibration'' problem, several strategies have been proposed, mainly relying either on good initializations of the MCMC chains to reduce the equilibration time, such as the (persistent) contrastive divergence method~\cite{hinton2002training,tieleman2008training} and its refinements~\cite{du2019implicit}, or on perturbative methods to approximate the gradient~\cite{gabrie2015training,kim2016deep,kumar2019maximum}.
\\
A very successful alternative to avoid altogether worrying about thermalization of the sampling processes during training
has been to give up encoding the empirical data statistics in the Boltzmann weight of a given energy model, as it is the normal EBMs' modus operandi.
The idea is instead to adjust the parameters of that model to reproduce a dynamical path able to generate new data points~\cite{sohl2015deep}, using an idea inherited from out-of-equilibrium thermodynamics.
That is, sacrificing the usual equilibrium-like interpretability of the EBMs for the sake of a fast and high quality generation performance. This idea is at the core of the so-called diffusion models~\cite{ho2020denoising}.

Furthermore, recent experimental studies
in several EBM families have shown that such a dynamic path is also printed in EBMs when non-convergent samplings are used to train the models, and that this
can be easily exploited to generate high-quality samples
in short times both in sampling and training~\cite{nijkamp2019learning,Nijkamp_Hill_Han_Zhu_Wu_2020,decelle2021equilibrium,muntoni2021adabmdca,liao2022gaussian}.  It was shown  that this optimal regime for sample generation of EBMs is achieved by forcing out-of-equilibrium sampling of the model during training, specifically by always initializing the Markov chains with random conditions and iterating the MCMC process for only a few steps $k$.
The key observation was that particularly high-quality samples are generated when sampling the trained model for exactly the same time and initial distribution that was used to train the model~\cite{nijkamp2019learning,decelle2021equilibrium}.
The apparition of this strong non-monotonic performance in the sample generation can be directly related to the MCMC mixing time, which typically increases with learning progress~\cite{decelle2021equilibrium}.
In fact,  the same work showed that most of the instabilities and difficulties typically encountered in training RBMs are directly related to this out-of-equilibrium regime, since all these effects disappear when the number of Gibbs steps used to compute the gradient exceeds the thermalization time,
e.g. when an incredibly long number of steps per iteration is considered~\cite{decelle2021equilibrium}.

In this paper, we provide a first-principle analytical explanation of these previous experimental observations.
In our theoretical framework, we show that the Gibbs-Boltzmann distribution of models trained in the out-of-equilibrium regime poorly reproduce the target statistics of the dataset. We also show that non-convergent samples can be used to accurately reproduce certain statistics of the dataset, but only through a very specific dynamic path fixed by the training choice.
We first prove these results for generic EBMs. Then,
we analyze two simple EBMs, the Gaussian model and the Bernoulli Boltzmann machine in the early stages of training, where the effects introduced into the model by the non-convergent samplings can be analytically described along the entire training dynamics.
Finally, we illustrate this phenomenon by performing real learning in the inverse Ising problem.

\section{Preliminaries: Maximum-likelihood training of Energy-Based Models}
\label{sec-preliminaries}

Suppose $\bx$ is the value taken by a data point.
We consider an arbitrary energy-based model (EBM) in order to model its probability in terms of a Gibbs-Boltzmann distribution:
\begin{equation}
\label{eq:ptheta}
p_{\bt} (\bx) = \frac{e^{-E_{\bt} (\bx)}}{Z(\bt)} \quad \text{with} \quad Z(\bt) = \int d\bx\ e^{-E_{\bt} (\bx)},
\end{equation}
where ${E_{\bt} (\bx)}$ is a dimensionless energy function that depends on a set of parameters ${\bt=\lbrace \theta_i \rbrace}$.
We let $E_{\bt} (\bx)$ be an arbitrarily complex function (e.g. parameterized by a ConvNet or a single-layer Boltzmann machine), the only requirement
being that $p_{\bt} (\bx)$ must exist together with its moments as functions of
their conjugate variables $\bm{\theta}$.

Suppose now that we want to find a good set of parameters $\bt$ that minimizes the distance between $p_{\bt}(\bx)$ and the empirical distribution $\pemp(\bx)$ of a given set of training examples.
This is usually done by maximizing the log-likelihood:
\begin{equation*}
\like (\bt) = \av{\log p_{\bt}(\bx)}_{\pemp}=-\av{E_{\bt} (\bx)}_{\pemp}-\log Z(\bt),
\end{equation*}
where $\av{\cdot}_{p}$ denotes the average over a given distribution $p(\bx)$.
This function is intractable in general because so is $Z(\bt)$: as a consequence, $\like (\bt)$ is typically maximized through a (possibly stochastic) gradient ascent process. Its derivative with respect to parameter $\theta_i$ is
\begin{equation}
\label{eq:gradEBM}
\frac{\partial\like (\bt)}{\partial \theta_i} =\av{f_{\theta_i}(\bx,\bt)}_{ p_{\bt}}- \av{f_{\theta_i}(\bx,\bt)}_{\pemp},
\end{equation}
with the local `forces' in the phase space of parameters\\
\begin{equation}\label{eq:fdef}
f_{\theta_i}(\bx,\bt)\equiv\frac{\partial E_{\bt}(\bx)}{\partial \theta_i}.
\end{equation}
Similarly, the Hessian matrix  is given by the components
\begin{eqnarray*}
\label{eq:hessian}
\nonumber
& H_{ij}(\bt)
\equiv \frac{\partial^2\like (\bt)}{\partial \theta_i \partial \theta_j}
=\av{\frac{\partial f_{\theta_j}(\bx) }{\partial  \theta_i}}_{\! p_{\bt}} - \av{\frac{\partial f_{\theta_j}(\bx) }{\partial  \theta_i}}_{\!\pemp\!}
\\
&\!-\!\av{f_{\theta_i}(\bx,\bt) f_{\theta_j}(\bx,\bt)}_{\! p_{\bt}}\!
+\!\av{f_{\theta_i}(\bx,\bt)}_{\! p_{\bt}\!}\av{f_{\theta_j}(\bx,\bt)}_{\! p_{\bt}\!}\!.\!
\end{eqnarray*}
This means that in models where all $f_{\theta_j}(\bx,\bt)$ are a function only of $\bx$ and not of the $\bt$
(e.g. $E_\bt$ is a linear function of the parameters),
the optimization problem is always convex (i.e. the Hessian is negative semi-definite). This occurs because the covariance matrix is positive semi-definite.
We will later analyze various models that satisfy this condition.

In practice the gradient \eqref{eq:gradEBM} is approximated by a stochastic sampling process (MCMC or Langevin dynamics). Thus, it is correctly estimated as long as this stochastic process converges to equilibrium.
However, such a sampling process must be repeated each time a gradient ascent step is performed.
As a consequence, ensuring convergence at each training step is not only tedious, but also infeasible in terms of the time required to train an adequate EBM.
A shortcut was proposed by Hinton: a few Gibbs sampling steps are sufficient if the sampling process is initialized wisely.
This is the idea underlying the recipe of contrastive divergence (CD), where the chains are always initialized at the training samples in the minibatch,
or of its refinement, persistent CD (PCD), where the configurations reached after a sampling process are kept as initial configurations for the subsequent sampling process after an update of the parameters.

The effects introduced into EBMs by the lack of convergence of these Markov chains at each update of the training, have gone largely unnoticed for decades, and only recently has it become clear that they introduce memory effects on the training process into models.
Even more interesting is the observation that
this issue can be exploited to our advantage:
non-persistent short MCMC runs are
the most efficient strategy to train good generative EBMs,
admittedly at the cost of forgetting the link between the equilibrium properties of the family of models parametrized by $\bt=\lbrace \theta_i \rbrace$,
and the statistics of the dataset encoded in $\pemp$.
In the next section, we provide an analytical explanation of these results.

\section{Finite-time sampling: The master equation} \label{sec:master_eq}

The implications of using short, non-convergent MCMC runs to estimate the gradient of Eq.~\eqref{eq:gradEBM}
in the training of EBMs can be formalized and proved generically as follows. For simplicity, we assume that the number of possible states of our model is finite; generalization to continuous variables is straightforward and this case will be considered in an example later.
We also assume continuous time, and provide the discrete time description in the Appendix~\ref{ap:discrete}.

We consider a MCMC sampling for a given model ${E_\bt}$, and denote $p_a(k)$ as the probability of being at the $a$-th possible configuration at time $k$. From now on, we denote the {\em sampling} time by $k$ to distinguish it from $t$, the {\em training} time.
Firstly, the {\em master equation} describes the time evolution of the probability of the system to occupy each of the states $a$
\begin{eqnarray}
\frac{\der p_a(k)}{\der k} = \sum_{b (\neq a)}w_{b \to  a} \, p_b{(k)}
-\sum_{b (\neq a)}w_{a \to  b} \, p_a{(k)},
\label{eq:mastereqcon}
\end{eqnarray}
where ${w_{a\to b}}$ is the transition probability of going from state $a$ to $b$ for that model.
These transition probabilities are independent of $k$ as the process is by definition Markovian.
We introduce the vector notation $\bm p{(k)}$, whose $a$-th component is $p_a{(k)}$. The evolution of the probabilities can then be described by a {\em stochastic evolution} matrix $\bU$ ($\U_{ab}=w_{b \to a} \: \text{for } a\neq b$ and $\mathcal{U}_{aa}
=-\sum_{b (\neq a)} w_{a \to  b}$),
as
\begin{equation}
\label{eq:evolp}
\dot{\bm p}{(k)}=\bU {\bm p}{(k)}
\quad \Rightarrow \quad
{\bm p}{(k)}
=e^{k{ \bU }}\,{\bm p}_0 ,
\end{equation}
where $\bm p_0$ represents the initial distribution (i.e. ${\bm p}(0)$).

Secondly, if the Markov chain is reversible, or equivalently satisfies detailed balance (which is the typical situation in EBMs), then the matrix $\bU$ is guaranteed to have a spectral decomposition.
Therefore, $\bm p(k)$ can be explicitly obtained from \eqref{eq:evolp} through a spectral expansion.
Let $\{\lambda_\alpha\}$ be the eigenvalues of $\bU $, and $\bm u_{\alpha}$ the associated eigenvectors.
Rewriting the initial distribution $\bm p_0$ in this basis,
we have
\begin{equation}
\bm p_0 = \textstyle\sum_\alpha c_\alpha {\bm u}_\alpha
\quad \Rightarrow \quad
{\bm p}(k)
=\textstyle\sum_\alpha c_\alpha e^{\lambda_{\alpha}k} \bm u_{\alpha},
\end{equation}
with the coefficient $\lbrace c_\alpha \rbrace$ encoding the memory of the initial condition of the sampling run.
The probability $\bm{p}{(k)}$ must be normalized at all $k$ (i.e. $\sum_a p_a{(k)}=1$), hence at least one eigenvalue must be $0$, and the rest must be necessarily negative for the long-time probabilities to remain bounded.

Thirdly, when $\bU$ is also irreducible (i.e. any configuration can be reached from any other in a finite number of steps), it follows from the Perron-Frobenius theorem that the $0$ eigenvalue must be unique (namely, $\lambda_0$)~\cite{keizer1972solutions}.
Furthermore, one can always choose $c_0=1$,
which imposes from the long-time normalization that ${\sum_a u_{0,a}=1}$, and consequently for the initial normalization that all the eigenvectors with negative eigenvalues must fulfill $\sum_a u_{\alpha,a}=0$ with $\alpha\neq0$. Under these constraints, we have then
\begin{equation}
\label{eq:probev}
 {\bm p}{(k)}=\bm u_{0}+\sum_{\alpha>0} c_\alpha\, e^{-k| \lambda_\alpha|}\bm u_{\alpha}.
\end{equation}
An analogous expression is obtained for discrete sampling steps in Appendix~\ref{ap:discrete}.
By construction,  it is guaranteed that the stochastic dynamics converges to the distribution $\bm u_0$ after a long time. Moreover, the {\em mixing} time is $\kappa_\mathrm{mix}\!=\!1/| \lambda_1|$, with $\lambda_1$ the second smallest absolute eigenvalue. Up to this point, the discussion was completely generic.

Fourthly, we can now explicitly focus on EBMs by forcing the long-time distribution ${\bm u}_0$ to coincide with the Gibbs distribution~\eqref{eq:ptheta}, i.e. $u_{0,a}= p_{\bt}(\bx^{(a)}) \equiv \exp\caja{-E_\bt (\bx^{(a)})}/Z(\bt)$, where  $\bx^{(a)}$ are the values of the variables in configuration $a$.
To be consistent, ${\bm u}_0$ must be the solution of the steady-state condition $\bU (\bt)\,\bm p_{\bt}=\bm 0$ (cf.~\eqref{eq:evolp}). It is in particular directly satisfied if we impose the {\em detailed balance} condition $w_{a \to  b}\, e^{-E_\bt(\bx^{(a)})}=w_{b \to  a}\,e^{-E_\bt(\bx^{(b)})}$
to hold, which has previously been imposed for the spectral expansion to exist. Standard choices for the update rules, such as the Metropolis-Hastings or the Heat Bath dynamics satisfy detail balance.

Finally, under this choice of transition probabilities, and the absence of some pathological cycles in the dynamics, we have the guarantee that the sampled
distribution will converge to the equilibrium one with increasing $k$ as
\begin{equation}
  {\bm p}{(k)}=\bm p_{\bt}+\bm \Lambda_{\bt}\paren{k,\bm p_0},
  \label{eq:dynU}
\end{equation}
where $\bm\Lambda_\bt$ is the term $\sum_{\alpha > 0}$ in Eq.~\eqref{eq:probev}, a decay function of $k/\kappa_\mathrm{mix}$ at long times.
One cannot compute this function $\bm\Lambda_\bt$ in general, but we know it depends only on
the initial conditions $\bm p_0 \leftrightarrow \lbrace c_\alpha \rbrace$, the run duration $k$, and the evolution matrix ${\bU (\bt)}$. Then, the expected value for $\bm f_{\bt}(\bx)\equiv\bm\nabla E_{\bt}(\bx)$ in the log-likelihood gradient (cf.~Eqs.~\eqref{eq:gradEBM}-\eqref{eq:fdef}), estimated after a sampling time $k$, is given by
\begin{eqnarray}
\label{eq:finitek}
\av{\bm f_{\bt}}_{k, \bm p_0}&\!\equiv\!&\textstyle\sum_a \bm f_{\bt}\textstyle\left(\bx^{(a)}\right)\, p_a {(k)}\\\nonumber  &\!=\!& \av{\bm f_{\bt}}_{ p_{\bt}}\!+\!\textstyle\sum_a \bm f_{\bt}\left(\bx^{(a)}\right)\Lambda_{\bt,a}\paren{k,\bm p_0}.
\end{eqnarray}
In other words, by approximating the gradient with a non-convergent MCMC sampling, one computes the real model average plus a correction that depends on $k$, $\bm p_{0}$, and the model parameters $\bm \theta$.
Therefore, at each training update, the update rule of Eq.~\eqref{eq:gradEBM} is replaced by
\begin{eqnarray*}
\label{eq:gradEBM-fink}
\left.\frac{\partial\like (\bt)}{\partial \theta_i}\right|_{k,\bm p_0}\!=\!\av{ f_{\theta_i}}_{ p_{\bt}\!}\!-\!\av{ f_{\theta_i}}_{\pemp\!}\!+\!\sum_a\! f_{\theta_i\!}\textstyle
\left(\bx^{(a)}\right)\textstyle\Lambda_{\bt,a}\!\paren{k,\bm p_0}.
\end{eqnarray*}
Moreover, this implies
that the fixed point $\bt=\bt^*$ of this dynamics (i.e. $\nabla \mathcal{L}=0$), if it can be reached,
satisfies
\begin{equation}
\label{eq:fixed-point-EBM}
    \av{\bm f_{\bt^*}}_{\pemp}\!=\!\av{\bm f_{\bt^*}}_{ p_{\bt^*}}\!+\!\textstyle\sum_a \bm f_{\bt^*}\left(\bx^{(a)}\right)\Lambda_{\bt^*\!,a}\!\paren{k,\bm p_0}.
\end{equation}
\vspace{-0.5cm}
\begin{theorem}
    \label{thm:markov}
    We consider
    a training of an EBM that follows the gradient in Eq.~\eqref{eq:gradEBM} estimated with nonconvergent, irreducible, and reversible Markov chains.
    If it reaches a dynamic fixed point where $\bm \nabla \mathcal{L}=\bm 0$, then, the trained model generates samples consistent with the statistics of the training set -- given by $\av{\bm\nabla E_{ \bt^*}}_{k,\bm p_0}\!=\!\av{\bm\nabla E_{ \bt^*}}_{p_\mathcal{D}}$ --, but only if the samples are drawn via the same identical stochastic process used for the training: i.e., the same update rules, the same initialization, and the same Markov chain length.
\end{theorem}
The proof results from combining the Eqs.~\eqref{eq:finitek}-\eqref{eq:fixed-point-EBM} using a model with parameters ${\bt^*}$ trained with a given sampling dynamics and $(\bm{p}_0,k)$. Now, the samples generated with this model, with ${\bU({\bt^*})}$, an arbitrary initial distribution ${\bm q_0}$ and after a time $k^\prime$, must satisfy the following statistics:
\begin{align}
\label{eq:errorfinitek}
&\av{\bm f_{\bt^*}}_{k^\prime \! , \bm q_0}
=\av{\bm f_{\bt^*}}_{\pemp}+\bm{D}_{\bt^*}(k^\prime,\bm q_0,k,\bm p_0),&
&
\end{align}
with the mismatch given by Eq.~\eqref{eq:probev}:
\begin{align*}
    &
    \bm{D}_{\bt^*\!}(k^\prime,\bm q_0,k,\bm p_0)\!
    = \! \textstyle\sum_{a} {\bm f}_{\bt^*\!,a} \left[ \Lambda_{\bt^*\!,a}(k^\prime \!,{\bm q}_0)\!
    -\! \Lambda_{\bt^*\!,a} (k,{\bm p}_0) \right]
    \\
    &\!=\!\textstyle\sum_{a,\alpha}\! \bm f_{\bt^*}(a) \paren{c_\alpha({\bm q_0})e^{-k^\prime|\lambda_\alpha|}\!- \! c_\alpha(\bm p_0)e^{-k| \lambda_\alpha|}} u_{\alpha,a}.&
\end{align*}
and the model-dependent $\lbrace \lambda_\alpha (\bt^*),\bm u_\alpha (\bt^*) \rbrace$.
This implies that, if we choose
the same sampling time ${k^\prime \!=k}$,
the same  initialization probabilities ${\bm q_0= \bm p_0}$,
and the same dynamic rules, as used during learning, then we will match perfectly the empirical statistics. This statement holds for all derivatives of the energy function $E_{\bt^*}$, with respect to the model parameters.
Also, a long sampling ($k^\prime\!\to\!\infty$) might sharply differ from empirical expectation when the training $k$ is much shorter than the mixing time $\kappa_{\text{mix}}$. This mismatch stems from the remaining factors ${\lbrace -c_\alpha({\bm p}_0) \, \exp^{-k \vert \lambda_\alpha \vert} \rbrace}$ in Eq.~\eqref{eq:errorfinitek}, which would vanish in the opposite limit $k \gg \kappa_{\text{mix}}$.

The results of Theorem~\ref{thm:markov} also hold before reaching the fixed point [i.e. for trained nonconvex EBMs or unfinished trainings of convex EBMs] but with several limitations.
\begin{theorem}
    \label{thm:markovbis}
    We consider a non-convergent EBM trained
    to a point where $\bm\nabla\mathcal{L} (\bt^\dagger)=\bm\varepsilon \neq \bm 0$.
    Then, the samples generated by the model $\bt^\dagger$ (with the same stochastic process used for training --the same $k$, the same initialization $\bm p_0$ and update rules) reproduce the statistics of $\av{\bm\nabla E_{ \bt^\dagger}}_{p_\mathcal{D}}$ of the training set with an error $\bm\varepsilon$. Thus, the error committed after time $k$ decreases (in average) as the training approaches its fixed point. Moreover, for a conveniently small $\bm \varepsilon$ (which depends on $\bm p_0$, $k$ and $\bm f_{\bt^\dagger}$), it is possible to find a sampling time $k^\dagger\neq 0$ at which the error is zero. This $k^\dagger$ moves toward the length $k$ of the Markov chains used for training as the training approaches its fixed point.
\end{theorem}
This second theorem is straightforward from the definition of the gradient and Eq.~\eqref{eq:errorfinitek}.
Then, upon  sampling model $\bt^\dagger$ with the same stochastic matrix $\bU (\bt^\dagger)$, the absolute error of $f_{\bt{^\dagger}}$ averaged on samples generated at a time $k^\prime$ from the initial distribution $\bm{q}_0$ is:
\begin{eqnarray*}
\textstyle\left|\av{\bm f_{\bt^{\dagger}}(\bx)}_{k^\prime,\bm p_0} \! -\av{\bm f_{\bt^\dagger}(\bx)}_{\pemp}\right| \! =\left|\bm\varepsilon + \bm{D}_{\bt^\dagger}(k^\prime,\bm q_0,k,\bm p_0)\right|,
\end{eqnarray*}
and $\bm{D}_{\bt^\dagger}(k^\prime,\bm q_0,k,\bm p_0)$ vanishes for $k^\prime=k$ and $\bm q_0=\bm p_0$. The second statement stands from the fact that $\bm{D}_{\bt^\dagger}$ is a smooth function of $k^\prime$ that (for $k$ short compared to $\kappa_\mathrm{mix}$) changes sign around $k^\prime=k$.
This means that the error of $\bm f_{\bt})$ will be strictly zero at $k^\dagger(\varepsilon)\neq k$ for $\bm\varepsilon$ sufficiently small,
which progressively moves towards $k$ as the learning progresses ($\bm\varepsilon \to 0)$.
In a typical MCMC scheme, time is always an integer variable, so this optimum might not be easy to access. However, in Langevin dynamics this is possible, as we show in a later example.

\section{Solvable models}
In the next sections we will examine two simple classes of EBMs for which:
\textit{(i)}~the correction
$\bm \Lambda_{\bt}(k,\boldsymbol{p}_0)$
introduced by the non-convergent Markov chains can be explicitly computed and \textit{(ii)},~its effect on the inferred model parameters can be explicitly quantified along the entire learning process.
Specifically, we focus on the inverse problem for families of exponentials with quadratic arguments:
Given a dataset drawn from an empirical distribution $\pemp(\bx)$, we wish to infer the best parameters  ${\bt}\equiv\{\bJ,\bh\}$ of the target distribution
with an energy function of the type
\begin{align*}

    \textstyle E_{\bt} (\bx)&=\frac{1}{2} \textstyle\sum_{ij} x_i J_{ij} x_j + \sum_i h_i x_i.

\end{align*}
For simplicity we consider $h_i=0$ in the following, but a generalization with external fields/biases is straightforward.
As in the general case in Sec.~\ref{sec-preliminaries}, we estimate the parameters by maximizing the log-likelihood of the model evaluated on the dataset, which reads
\begin{equation}
    \like (\bt) = -\frac{1}{2} \sum_{ij} J_{ij}\av{x_i x_j}_{\pemp} - \log{Z(\bt)}.
\end{equation}
Again, the classical approach to maximizing log-likelihood is based on calculating its gradient
\begin{align}
    \frac{\partial \like (\bt)}{\partial J_{ij}} &= -\paren{\av{x_i x_j}_{\pemp} - \av{x_i x_j}_{p_{\bt}} }\!, \label{eq:gradJ}
\end{align}
The average over ${\bm p}_\bt$
is usually calculated using the heat-bath MCMC algorithm. In the next two sections, we analytically investigate the effect of replacing the equilibrium averages with estimates obtained by a non-convergent MCMC process.
We  first consider the Gaussian model and then the high-temperature regime of the inverse Ising problem.

\section{Exact results for the Gaussian model}

In the Gaussian case, we consider data points $\bx \in \mathbb{R}^D$. To be well-defined, the coupling matrix $\bJ$ must be positively definite, namely all its eigenvalues must be strictly positive. Under this condition, and given the convexity of this likelihood function, we can explicitly compute the correlation function, and the stationary solution of Eq.\eqref{eq:gradJ} reads:
\begin{equation*}
    \av{\bm x \bm x^{\text{T}}}_{p_{\bt}}= \bm J^{-1}
    \; \stackbin[k\to\infty]{[\nabla \mathcal{L}=0]}{\Longrightarrow} \;

    \bm J^{*} = \av{\bm x \bm x^{\text{T}}}_{\pemp}^{-1} .
\end{equation*}
Despite the maximum likelihood point can be found analytically, we now suppose to compute it through a gradient ascent method, with a suitable learning rate. Moreover, the gradient can be projected on the eigenbasis of $\bJ$. Denoting {$\lbrace \bm{u}^\alpha \rbrace$} the eigenvectors, we have that $J_{ij} = \sum_\alpha u_i^\alpha J_\alpha u_j^\alpha$. Following~\cite{decelle2017spectral,decelle2018thermodynamics,decelle2021restricted}, we can write a set of evolution equations for the eigenvalues {$\lbrace J_\alpha \rbrace$} and for the rotation {$\Omega_{\alpha \beta}$} of the eigenvectors
\begin{align}
    \frac{\der J_\alpha}{\der t} &= -\hat{c}_{\alpha \alpha} + 1/J_\alpha, \label{eq:grad_Ja}\\
    \Omega_{\alpha \beta} & {\equiv} \bm{u}^\alpha \frac{\der\bm{u}^\beta}{\der t} = \frac{\langle x_\alpha x_\beta \rangle_{p_\bt} - \hat{c}_{\alpha \beta}}{J_\alpha - J_\beta}, \label{eq:grad_ua}
\end{align}
where we defined $\hat{c}_{\alpha \beta} = \sum_{i,j} u_i^\alpha u_j^\beta \av{x_i x_j}_{\pemp}$ and $x_\alpha = \sum_i x_i u_i^\alpha$, the projection on the eigenmodes of $\bJ$ of the correlation matrix and $\bx$, respectively.
In this formulation, we see that the eigenvectors have to diagonalize the correlation matrix to cancel the rotation,
i.e. each vector
$\bm{u}^\alpha$ will align itself with the principal directions of the dataset, and
the eigenvalues will eventually adjust to $J_\alpha = 1/\hat{c}_{\alpha \alpha}$.
\begin{figure*}[h!]
    \centering\includegraphics[height=5.4cm,trim=0 20 0 0]{figs/EvolFiniteK_inset.png}
    \caption{\textbf{Left:} Evolution of $J_\alpha^{(k)}(t)$ for different values of $k$ (different colors) and initial conditions for $J_\alpha^{(k)}(0)$. \textbf{Right:} Convergence for a large training time $t$, where we clearly show that the higher $k$, the faster the convergence; dash-dotted lines result from the numerical integration of Eq.~\eqref{eq:grad_Jak}, full lines show the exponential fit $J_\alpha^{(k)} (t) - J_\alpha^{(k),\infty}\! \sim\! A\exp(-t/\tau)$ with $\tau$ given by Eq.~\eqref{eq:taugauss}. \textbf{Inset:} the asymptotic values $J_\alpha^{(k)}(t\!\to\!\infty)$ as a function of $k$. The gray shaded region corresponds to the non-physical solutions and the vertical line to the limit $k = (\hat{c}_{\alpha \alpha} - \av{x_\alpha^2}_0)/2$.}
    \label{fig:EvolFiniteKInset}
\end{figure*}
\\
However, we are interested here in the case where non-convergent MCMC chains are used to approximate the correlation function $\av{\bm x \bm x^{\text{T}}}_{p_{\bt}}$, and what it changes for Eqs.~\eqref{eq:gradJ}-\eqref{eq:grad_Ja}-\eqref{eq:grad_ua}. We assume for simplicity that the dataset is centered with zero mean, $\av{x_i}_{\pemp} = 0$.
The evolution of the MCMC chain is described by the following Langevin equation
\begin{equation}\label{eq:gradientk}
    \frac{\der x_i(k)}{\der k} \!=\! -\!\textstyle \sum_j J_{ij} x_j (k)\! + \!\eta_i(k)
\end{equation}
with $\bm \eta$ a Gaussian white noise of zero mean and variance $\overline{\eta_i (k) \eta_j(k')} = 2 \delta_{ij} \delta(k-k')$, where $\overline{(\cdot)}$ denotes the average with respect to different independent trajectories. These equations can be decoupled by projecting them on the eigenmodes of $\bJ$, leading upon sampling to
\begin{equation*}
    x_{\alpha}(k)=x_{\alpha}(0)e^{-J_\alpha k}+\textstyle\int_0^k \mathrm{d}\tilde{k} \ e^{J_\alpha (\tilde{k}-k)}\eta_\alpha(\tilde{k}).
\end{equation*}
We see that the convergence properties are dominated by the eigenvalues $J_\alpha$. Hence, for $k \lesssim J_\alpha^{-1}$, the MCMC chain will be far away from its equilibrium value. We can obtain a similar result for the correlation function; let us directly look at the average over realizations at time $k
\begin{equation}
    \av{x_{\alpha}x_{\beta}}_{k,\bm p_0} \!=\!\av{x_\alpha x_\beta}_0 e^{-(J_\alpha+J_\beta)k}\!+\!\frac{\delta_{\alpha\beta}}{J_\alpha}\paren{ 1\!-\! e^{- 2 J_\alpha k}}, \label{eq:corrLangeving}
\end{equation}
where $\av{\cdot}_0$ is the average over the initialization distribution. Thus, the asymptotic value is given by $\av{x_{\alpha}x_{\beta}}_\infty = \delta_{\alpha \beta}/J_\alpha$.
\\
We now analyze how the learning dynamics are affected when the model average $\av{x_{\alpha}x_{\beta}}_{p_\bt}$ is replaced in the gradient by $\av{x_{\alpha}x_{\beta}}_k$, i.e. the averages computed at fixed time $k$ when the chains are initialized at a fixed random initial distribution $p_0(\bx)$. We first consider only the evolution of the eigenmodes given by Eq.~\eqref{eq:grad_Ja}, and keep the eigenvectors $\bm{u}^\alpha$ fixed. Under these conditions, the gradient calculated with non-convergent chains is quite different to the real one
\begin{align*}
    &\left.\left(\frac{\partial \mathcal{L}}{\partial \bm{J}}\right)_{\alpha\beta}\right|_{k,{\bm p}_0} ={-\left[\left\langle x_\alpha x_\beta \right\rangle_{\pemp}-\av{x_{\alpha} x_{\beta}}_{{k,\bm p_0}} \right]}
    \\
    &=\left(\frac{\partial \mathcal{L}}{\partial \bm{J}}\right)_{\alpha\beta} + \paren{\av{x_\alpha x_\beta}_0 - \frac{\delta_{\alpha \beta}}{J_\alpha}}e^{-\paren{J_\alpha+J_\beta} k}.
\end{align*}
Note that the gradient now depends on both the sampling time
$k$ and the initialization ${\bm p}_0$ of the chains.
Then the evolution for the eigenvalues of $\bm J$ in Eq.~\eqref{eq:grad_Ja} in the finite $k$ case transforms into
\begin{equation}\label{eq:grad_Jak}
    \frac{\der J_\alpha^{(k)}}{\der t}\! =\! -\hat c_{\alpha\alpha} + \frac{1}{J_\alpha^{(k)}}  \!+\!  \left(\av{x_\alpha^2}_0 \! -\! \frac{1}{J_\alpha^{(k)} }\right) e^{-2 J^{(k)}_\alpha k},
\end{equation}
whose stationary solution at ${t \to \infty}$, $J_\alpha^{(k),\infty}$, satisfies
\begin{equation}
    \hat{c}_{\alpha \alpha} \!-\! \av{x_\alpha^2}_0 e^{- 2 J_\alpha^{(k),\infty} k}\! =\! \frac{1}{J_\alpha^{(k),\infty} }\! \left( 1\!-\!e^{-2 J_\alpha^{(k),\infty} k} \right). \label{eq:fixedpoint_gauss}
\end{equation}
\begin{theorem}
    \label{thm:conv_gauss}
    If $k > (\hat{c}_{\alpha \alpha} - \av{x_\alpha^2}_0)/2$ holds, then Eq.~\eqref{eq:fixedpoint_gauss} admits a strictly positive solution. Otherwise, $J_\alpha^{(k),{\infty}}$ is negative and the non-convergent learning procedure diverges.
\end{theorem}
This theorem is
proved by analyzing the function $g(J) = J \hat{c}_{\alpha \alpha}-1 + e^{- 2 J k}(1-J\av{x_\alpha^2 }_0)$, see the appendix~\ref{ap:thm_k} for more details. Therefore, we can define a lower bound for the sampling time $k$ below which the learning dynamics leads to an unphysical solution (unbounded Gaussian integrals).
From Eq.~\eqref{eq:fixedpoint_gauss}
we recover straightforwardly
the correct solution $J_\alpha^* = 1/\hat{c}_{\alpha \alpha}$
in the limit $k \to \infty$.

Equation~\eqref{eq:grad_Jak} can be integrated numerically, and the results for various $k$ and initial conditions are shown in
Fig.~\ref{fig:EvolFiniteKInset} (left).
The smaller $k$, the slower the convergence to the fixed point $J_\alpha^{(k),\infty}$.
Indeed, one can
show that the typical training time $\tau$ required to reach the stationary point decreases with $k$ as
\begin{equation}
\label{eq:taugauss}
\frac{1}{\tau} \!\sim\! k\caja{\av{x_\alpha^2}_0 \!-\!\frac{ e^{-2k J_\alpha^{(k),\infty}}}{J_\alpha^{(k),\infty}}} \!+\! \frac{\left(1-e^{-2k J_\alpha^{(k),\infty}} \right)}{2\left(J_\alpha^{(k),\infty}\right)^2},
\end{equation}
by simply assuming an exponential relaxation for a large training time $t$, i.e., $J_\alpha^{(k)} (t)\! \sim\! J_\alpha^{(k),\infty} + A\exp(-t/\tau)$ in Eq.~\eqref{eq:grad_Jak}. Therefore, from a computational point of view, a smaller value of $k$ means by definition faster sampling, but the time required for the learning process to converge increases as $k$ decreases. Thus, the optimal value of $k$ results from a trade-off between these two contributions.
\begin{figure*}[h!]
    \centering
    \includegraphics[height=5.2cm,trim=0 20 0 0]{figs/EvalRot.png}
    \includegraphics[height=5.2cm,trim=0 20 0 0]{figs/GaussCorrError-new.png}
    \caption{\textbf{Left:} Numerical resolution of the learning dynamics in presence of two modes. The dash-lines represent the resolution using a convergent dynamics $k \to \infty$, while the plain ones correspond to $k=1$. \textbf{Right:} \textbf{(inset)} evolution of $J_\alpha^{(k)}(t)$ for $k=1$ at various stages of the learning.  \textbf{(Main)}  The error on the correlation function, $\mathcal{E}^{(2)}=\textstyle(\av{x_\alpha^2}_{k^\prime,\bm p_0}-\hat{c}_{\alpha\alpha})^2$, of the machine at different stages of the learning. We can clearly see the peak of best performance (of several orders of magnitude) at a $k^{\dagger}$ that moves to $k$ as the learning converges. The color lines match the circles on the inset. The horizontal dashed-line represent the asymptotic value of the error for $k\to \infty$.}
    \label{fig:EvolJResampling}
\end{figure*}

One can also study the effect of the non-convergent chains on the alignment of the eigenvectors of $\bm J$,
which we kept fixed until now.
For this, we rely on the simplified situation proposed in~\cite{decelle2021restricted}, where a reduced 2D subspace is considered in which the 2 first eigenvectors of $\bJ$ form an angle $\phi$ with the 2 principal directions of the dataset. In this case, one obtains a set of $3$ dynamical equations for $\phi(t)$, and the 2 eigenvalues $J_1(t)$ and $J_2(t)$
\begin{align*}
    \textstyle\frac{\der\phi}{\der t}  &\!=\! \textstyle\left[\sin(2\phi) (\hat{c}_1\!-\!\hat{c}_2)/2 \!+\! \av{x_1 x_2}_0 e^{- (J_1\!+\!J_2) k}\right]\!/\!(J_2\! -\! J_1),\\
    \textstyle\frac{\der J_1}{\der t} &\!=\! \textstyle J_1^{-\!1\!} \!-\! \hat{c}_1\! \cos^2\!(\phi) \!-\! \hat{c}_2\! \sin^2\!(\phi) \!-\! \left(J_1^{-\!1}\!-\!\av{x_1^2}_0\right) e^{\!-\!2 J_1 k}, \\
    \textstyle\frac{\der J_2}{\der t} &\!=\! \textstyle J_2^{-1\!}\! -\! \hat{c}_1\! \sin^2\!(\phi)\! -\! \hat{c}_2\! \cos^2\!(\phi) \!-\! \left(J_2^{-\!1\!}\!-\!\av{x_2^2}_0\right) e^{\!-\!2 J_2 k}.
\end{align*}
(see Appendix~\ref{ap:rot} for details) whose numerical solution is given in Fig.~\ref{fig:EvolJResampling}-.left.
In the limiting case $k\!\to\!\infty$, we see as expected the eigenvectors align (or anti-align) with the eigenvectors of the covariance matrix of the dataset and  $\phi\!\to\! 0 \text{ or } \pi$.
Yet, for finite $k$ this is not true, and an asymptotic finite $\phi$ is reached which embodies the sampling mismatch.

In Sec.~\ref{sec:master_eq}, we showed that generic EBMs trained with non-convergent MCMC samplings can be used to generate
samples that perfectly match the statistics of the dataset even before the learning has converged.
For that purpose, we only need to mimic the same dynamic process used for training in the generation phase. We can then
quantify this prediction in this solvable case,
given the
parameters of the model $J_\alpha^{(k)}(t)$ after a learning time $t$,
by computing the average correlations
through Eq.~\eqref{eq:corrLangeving}.
\begin{theorem}
     \label{thm:ReSamplingK}
     We consider a Gaussian model trained with non-convergent Markov chains of length $k$ along the modes $\alpha$ (i.e. keeping fixed the eigenvectors $\lbrace \bm{u}^\alpha \rbrace$), and parameterized by $\lbrace J_{\alpha}^{(k),\infty} \rbrace$. It
     exactly reproduces the correlations of the dataset when sampled with the same dynamic process used in the training: same transition rules, initialization distribution, and chain length.
\end{theorem}
The proof is straightforward by inserting the fixed point equation~\eqref{eq:fixedpoint_gauss} in~\eqref{eq:corrLangeving}. In particular, let us generate new data $\bm y$ with the trained model.
When $k'=k$, the averaged correlation along  each mode $\alpha$
\begin{equation*}
    \av{y_\alpha^2}_{k,\bm p_0}=\hat{c}_{\alpha \alpha} + \left( \av{y_\alpha^2 }_0 -  \av{x_\alpha^2 }_0\right)e^{-2 J_\alpha^{(k),\infty} k},
\end{equation*} match those of the dataset but only if the chains are initialized from the same distribution used during the training, i.e, same $\av{x_\alpha^2 }_0$. In Fig.~\ref{fig:EvolJResampling}--right, we show the mean squared error $\mathcal{E}^{\left(2\right)}$ in the correlations committed at different stages of learning, i.e., using the numerically integrated $J_\alpha^{(k)}(t)$. As anticipated in Sect. 3, one can also reproduce the statistics of the dataset before convergence, but at an effective sampling time $k^{\dagger}$ that moves to $k$ as learning progresses.

As a final remark, from Eq.~\eqref{eq:corrLangeving} we see that the cross-correlation terms $\av{x_\alpha x_\beta}_{k}$ for $\alpha \neq \beta$ are nonzero only when the initial conditions are de-correlated along both modes, and thus the eigenvectors $\bm{u}^\alpha$ easily align with the principal directions of the dataset. In purely random initializations and to some extent in the PCD recipe, the cross terms are zero, but in the CD scheme this is not the case in the initial stages of learning, which should affect the rotational dynamics of the eigenvectors of the coupling matrix.

\begin{figure*}[h!]
    \centering\includegraphics[height=5.5cm,trim=0 20 0 0]{figs/figBM.pdf}
    \caption{\textbf{Left:} Error over the covariance matrices $\mathcal{E}^{(2)}=\sum_{i<j}\textstyle(\av{x_i x_j}_{k^\prime,\bm p_0}-\av{x_i x_j}_{p_\mathcal{D}})^2 \slash \binom{N}{2}$ between the training set and data generated at different learning ages vs the generation time $k'$, for a BM trained with data sampled from a $2$D ferromagnetic Ising model with $N=7^2$ at $k=5$.  \textbf{Inset:} Same quantity computed for different $k$ for the oldest machine ($t=10^5$), generated at $k^\prime=k$. Similarly to Fig. \ref{fig:EvolJResampling} (right), the error shows a minimum at $k' = k$. \textbf{Right:} Error over the learned couplings w.r.t. the true ones as a function of the training time $t$ for different values of $k$. A comparison with a BM trained using a persistent chain (PCD) is also shown.
    The learning rate is set to $\gamma=10^{-2}$.
    \label{fig:BM_noneqeffect_2dFerro}}
\end{figure*}

\section{Exact results on the Boltzmann machine}
The previous results can be extended to the case of binary variables, such as the classical Ising inverse problem, or the Boltzmann machine (BM)~\cite{zecchinareviewinverseising}. In this context, we have $\bx = \{\pm 1\}^N$, for which the distribution is defined for all $\bJ$.In general, the learning of the BM cannot be solved exactly, since there is no general closed-form expression for the correlation function. However, we can deal with the so-called mean-field case (MF), where the coupling matrix scales with $1\slash N$. In this regime, BM can be solved analytically using linear response theory in many MF regimes~\cite{10.1162/089976698300017386,Nguyen_2012,Ricci-Tersenghi_2012}.
It is then possible to derive the evolution equation for the magnetization and the correlation function of a given MCMC dynamics .
Following~\cite{suzuki1968dynamics}, we first define the heat-bath dynamics transition matrix as the probability of a given spin $s_j$ when all others are fixed: $W_j(s_j^{t+1} | \{s_i^{t}\}_{i \neq j})$. Then, using the master equation as Eq.~\eqref{eq:mastereqcon}, it is possible to write
the evolution equation of the moments of the spins, which in the case of the correlations gives
\begin{equation*}
    \frac{\der\!\av{s_i s_j}_{k}}{\der k} \!=\! -2 \av{s_i s_j}_{k}+ \av{s_i \text{tanh} (\bm J \bm s )_j}_{k}
    + \av{s_j \text{tanh} \left(\bm J \bm s \right)_i}_{k}.
\end{equation*}
The above equation is exact for any coupling matrix $\bJ$, but cannot be solved in general. In the MF regime, we can expand the $\tanh$ at first order and solve the set of linear differential equations. Since it is very similar to the Gaussian case, we only provide
the results. At leading order, the evolution upon sampling of the two-body correlations, in the eigenbasis of the couplings $\bm J$, can be written \begin{equation*}
    \frac{\der \av{x_\alpha x_\beta}_{k}}{\der k} \!=\! -\! \av{x_\alpha x_\beta}_k \left( 2-  J_\alpha - J_\beta \right) + 2 \delta_{\alpha \beta},
\end{equation*}
whose solution is
\begin{align*}
    \av{x_\alpha x_\beta}_k\! =\!
         \frac{\delta_{\alpha \beta}}{1-J_\alpha} \!+\! \caja{\av{x_\alpha x_\beta}_0\!+\!\frac{\delta_{\alpha \beta}}{1-J_\alpha}}  e^{-(2- J_\alpha- J_\beta) \frac{k}{\tau}}.
\end{align*}
This analytical description is valid only in the high-temperature regime, where ${\rm max}\{J_\alpha\} < 1$. Nevertheless, in this regime we can reproduce the same analysis of the Gaussian model and reach the same conclusion: BMs trained with non-convergent MCMC  to compute the gradient fit models that reproduce the training set statistics, but only when sampled following the dynamics of the training: same initial conditions, transition rules, and number of steps $k$.

\section{Numerical experiments}

This non-monotonic effect in the quality of the generated samples and its relation to the mixing times of the algorithm has been intensively numerically investigated in the case of RBM~\cite{decelle2021equilibrium}, yielding a picture fully compatible with the results obtained here. For example, the evolution of their ``quality estimators" with the sampling time are qualitatively very similar to those shown in Fig.~\ref{fig:EvolJResampling} for discrete sampling times. Here instead we choose to test our analytic predictions on the BM for the inverse Ising problem, where $\bm x \in \{ -1,1 \}^N$, and the energy function is again quadratic in the variables $\bx$: $E_\bt (\bx) = -\sum_{i<j} J_{ij} x_i x_j - \sum_i h_i x_i$. In the following experiment, the dataset is made of configurations drawn from a $2$D ferromagnetic Ising model, with $N=7^2$ spins, $J_{ij} = \beta$, $h_i = 0$, where $(i,j)$ are nearest neighbors on the $2$D lattice. We consider $\beta = 0.44$, i.e a temperature slightly above the critical point.
Learning is done with $k$-steps of nonconvergent heat-bath Markov chains initialized with random conditions. For $k=5$ and for machines obtained at different stages of the learning process, we show in Fig.~\ref{fig:BM_noneqeffect_2dFerro}--left the evolution of the error in the 2-point correlations for generated samples as a function of the number of MCMC steps used to synthesize them.  We clearly observe the same behavior as in the Gaussian case, with a best-quality peak developing at $k' = k$. Note that since $k,k'$ are now discrete, we are not able to see the high quality peaks occurring before convergence. We can also observe that, consistently with our previous predictions, the quality of the best generated samples (those obtained at $k'=k$ at the end of the learning process), does not depend on $k$. It is particularly remarkable that trainings with $k=1$ generate samples as good as those with $k=50$.
Finally, in Fig.~\ref{fig:BM_noneqeffect_2dFerro}--right, we show the error in the inferred couplings $\bJ$ for different values of $k$ as a function of training time. These curves illustrate how short times $k$ can lead to very poor inferred models,
while
models learned with short nonconvergent chains can still catch important information about the data that is useful for interpretive applications.

\section*{Conclusions}

Our results show that EBMs trained with non-persistent and non-convergent MCMCs to estimate the log-likelihood gradient are models with memory of their training scheme. We also show analytically how to exploit such memory effects to generate high-quality samples in a very short time, through a very precise sampling process that must perfectly match that of the training: i.e., same update rules, initializations and sampling time. EBMs trained with the contrastive divergence recipe fall into this category. However, they have the disadvantage that they cannot be used to generate new data because they provide minimal variation in the sample used for initialization. In this sense, we show that EBMs trained with short \emph{randomly} initialized MCMC runs can be effectively used as diffusion models for sample generation. \\On the other hand, we also prove that the equilibrium measure of an EBM trained in this way cannot correctly describe the statistics of the data and should be used with caution for interpretive applications. The lack of convergence of MCMC runs is particularly critical for datasets with an effectively low dimension, i.e., datasets that follow multimodal distributions, as mixing times become excessively long~\cite{bereux2022learning,decelle2021exact} to be thermalized by brute force.

\section*{Acknowledgements}

G.C., A.D. and B.S. acknowledge financial support by the Comunidad de Madrid and the Complutense University of Madrid (Spain) through the Atracción de Talento programs (Refs. 2019-T1/TIC-13298 for A.D. and G.C. and 2019-T1/TIC-12776 for B.S.), the Banco Santander and the UCM (grant PR44/21-29937), and Ministerio de Econom\'{\i}a y Competitividad, Agencia Estatal de Investigaci\'on and Fondo Europeo de Desarrollo Regional (FEDER) (Spain and European Union) through the grant PID2021-125506NA-I00. E.A. acknowledges support from the Swiss National Science Foundation by the SNSF Ambizione Grant No.~PZ00P2{\_}173962, and from the Simons Foundation (Grant No.~348126 to Sidney Nagel).