Computer Science Mathematics
Computer Science-2021.gz:143170
Computer Science
678f613cb914a392faa260844b86a7a588ff4615
Supervising the Decoder of Variational Autoencoders to Improve Scientific Utility

Probabilistic generative models are attractive for scientific modeling because their inferred parameters can be used to generate hypotheses and design experiments. This requires that the learned model provide an accurate representation of the input data and yield a latent space that effectively predicts outcomes relevant to the scientific question. Supervised Variational Autoencoders (SVAEs) have previously been used for this purpose, where a carefully designed decoder can be used as an interpretable generative model while the supervised objective ensures a predictive latent representation. Unfortunately, the supervised objective forces the encoder to learn a biased approximation to the generative posterior distribution, which renders the generative parameters unreliable when used in scientific models. This issue has remained undetected as reconstruction losses commonly used to evaluate model performance do not detect bias in the encoder. We address this previously-unreported issue by developing a second order supervision framework (SOS-VAE) that influences the decoder to induce a predictive latent representation. This ensures that the associated encoder maintains a reliable generative interpretation. We extend this technique to allow the user to trade-off some bias in the generative parameters for improved predictive performance, acting as an intermediate option between SVAEs and our new SOS-VAE. We also use this methodology to address missing data issues that often arise when combining recordings from multiple scientific experiments. We demonstrate the effectiveness of these developments using synthetic data and electrophysiological recordings with an emphasis on how our learned representations can be used to design scientific experiments.

An interpretable relationship between the observed covariates and model representation, while necessary, is often not sufficient in scientific settings. Many scientific applications require that the learned representation also be predictive of an auxiliary variable [9]. In the neuroscience research motivating this work, this is often a behavioral outcome [10], a genetic phenotype [11], or the presence of some disorder or disability [12]. Obtaining a predictive latent representation allows for researchers to identify patterns relevant to the auxiliary variables and establish causality through experimental modification [13]. Unfortunately, exclusively generative models often fail to represent the desired auxiliary variable, and are typically dominated by other irrelevant sources of variation. Returning to our motivating work as an example, neural dynamics associated with motion [14] and even blinking [15] are often substantially stronger than the auxiliary variables of interest.
One class of generative models called Variational Autoencoders (VAEs) can be encouraged to yield a predictive latent representation by including a supervision loss during training. This yields a standard Supervised Variational Autoencoder (SVAE), which has been commonly used in the machine learning community [16,17,18]. In those applications, the reconstruction loss of the generative model has been motivated as an effective method of improving predictive models, as the reconstructive loss is theoretically and empirically an effective regularization technique [19]. Based on this work it might seem that SVAEs satisfy both the generative and predictive criteria for scientific utility.
Unfortunately, the supervision loss in a standard SVAE (referred to as SVAE in the remaining sections) biases the encoder away from approximating the true posterior distribution of latent representations. We demonstrate this bias by analyzing the fixed points of the SVAE objective. The bias caused by the inclusion of the supervision loss has not been noted in previous work, in part because those applications primarily focused on obtaining a predictive model. When the generative aspects of SVAEs are viewed merely as a convenient and effective regularization technique, the scientific utility of the decoder parameters is irrelevant. However, this bias in the variational objective can have a profoundly negative impact on causal scientific experiments designed using the decoder parameters. In these applications, the model will give misleading conclusions on how modification of the observed covariates will influence the auxiliary variable via the latent representation. As our motivating work uses these generative decoder models as a means for informing scientific manipulations, addressing this issue is critical. Fig. 1: Relationship among the proposed models (SOS-VAE and SOS-DVAE, detailed in Sections III and IV, respectively), a generative model (VAE) and a predictive model (SVAE). Blue indicates a model pursues a lower reconstruction error (better generative performance), and red denotes that a model pursues a higher predictive performance. Both VAE and SOS-VAE try to find latent representations that only use information in the proper generative model family, with SOS-VAE trying to find a representation that is good at prediction within the model family. We note that the learned latent representation from SVAEs can be far from the generative model family, as we empirically show in Section V, and SOS-DVAE allows us to vary between optimal predictive performance and faithfulness to the generative model.
We remove this previously unobserved bias in SVAEs by developing a novel optimization framework using second-order gradient techniques [20,21]. This Second-Order Supervision (SOS) framework maintains the interpretation of the variational encoder as an unbiased approximation to the posterior while inducing the latent representation to be predictive of an auxiliary variable (Fig. 1). This approach yields a model that possesses both properties required for scientific utility, a predictive latent space and an unbiased posterior approximation. We empirically demonstrate these properties on a real dataset which uses Local Field Potentials (LFPs) recordings [4,5] to predict a behavioral trait. We show that our framework yields more accurate predictions when compared to an exclusively generative model [22] while maintaining a more accurate variational approximation of the posterior when compared to an SVAE.
Previous work has shown that interpretable generative models can handicap a predictive objective as compared to a flexible neural network [23] due to the restrictive constraints imposed by the model family. Because of this, our SOS framework can cause degraded predictive performance relative to an SVAE when the variational approximation is unconstrained by the generative model. While sometimes this degradation is scientifically necessary and unavoidable, often a faint amount of bias is acceptable to obtain substantial gains in predictive ability. Our first extension provides a means to relax the constraint on the variational encoder, allowing it to maintain high predictive ability with minimal bias in the approximation. Using Kullback-Liebler (KL) divergences [24], we demonstrate on both synthetic and real data examples that this regularization approach improves predictive performance for a given level of posterior discrepancy compared to an SVAE-based model. We denote this method as Second-Order Supervision Double Variational Autoencoder (SOS-DVAE), and conceptually show the relationship to the other models in Fig. 1. We additionally show that this framework can be helpful in the context of missing data.
Thus, our work provides three major contributions to the field. First, we identify bias in SVAEs using a novel fixed point analysis, and show this bias has a detrimental impact on generative reconstruction. Second, we develop inference techniques to reduce or remove this bias in a computationallyfeasible manner resulting in a more realistic generative model. Finally, we demonstrate how the proposed models (SOS-VAE, SOS-DVAE) can be applied in scientific settings by isolating brain networks of predictive of traits of interest.
II. RELATED WORK Joint Factor Modeling. Joint modeling, such as probabilistic supervised principal component analysis [25] or supervised Gaussian processes [26], assumes that the observations and the outcome are independent given a latent representation. Once the prior distribution of the latent factors and the conditional distributions have been defined, statistical estimation is straightforward by maximum likelihood or Bayesian methods [27]. However, it has been demonstrated that joint models suffer under model misspecification, particularly when the number of learned factors is less than the true latent dimensionality [23,28]. Furthermore, the variance of the outcomes is often small relative to the variance of the observations, which leads to the outcome being poorly characterized [9].
Supervised Autoencoders. To address some of these limitations, SAEs can be used to effectively approximate standard factor models while including a predictive term. SAEs have been shown to ameliorate some concerns about model misspecficiation [28]. SAEs have been shown to give gains in many predictive applications [18,16]. Using an SAE can be viewed as a form of regularization and prevents the latent representation from over-fitting, which has been shown theoretically to enhance generalizability [19]. It has also been shown in deep learning that adding auxiliary tasks can act as a form of regularization [29,30,31]. Recent work has also developed strategies for unsupervised generation of the tasks [21]. However, none of these works evaluated how well the inferred latent space fit the generative model and instead focused solely on prediction.
Semi-Supervised Learning. The standard VAE [24] (referred to as VAE in the remaining contents) has been widely extended in semi-supervised learning tasks [32,33,34,35] using deep neural networks. These VAE-based semi-supervised models pursue a shared latent space for both inference and prediction, and have shown impressive predictive performance in some applications. However, they have the same limitation that the inferred latent space does not necessarily fit the generative model. In addition, most of existing methods were evaluated only on image data [32,33,34] or sequence data [35], which might not suitable for our motivating applications (predicting a characteristic, trait or behavior from electrophysiology). Importantly, these models are less interpretable due to the use of complicated generative network structures and thus not suitable for scientific hypothesis generation.
Second-Order Optimization. Finally, our novel learning technique will require the use of second derivatives to indirectly induce the learned latent variable model to yield a predictive posterior. This is difficult to implement directly via backpropagation. However, by incorporating computational tricks used in some meta-learning and self-supervision techniques [20,36,21], this can be efficiently implemented in modern learning platforms.

III. METHODS
We first introduce notation. Let {x i } i=1,...,N ∈ R p be N independent samples with associated outcomes {y i } i=1,...,N ∈ Y drawn from the true joint probability distribution p d (x, y). Given these data, we fit a joint model defined as x, s ∼ p θ (x, s) and y|s ∼ p ψ (y|s), where s ∈ R L is the L-dimensional latent space. To simplify theoretical derivations we assume that p θ and p ψ are infinitely differentiable, which holds with exponential family models that are commonly used in our applications. We view θ as parameterizing our generative model while ψ parameterizes the predictive model (the classifier) of the outcome given the latent space. We assume that a decoder p θ (x|s) can be derived directly from the likelihood of the true generative model. We approximate our true posterior p θ (s|x) by a variational encoder q φ (s|x) to enable variational inference [24]. In our model, as in SVAEs or any standard supervised model, the variational approximation is not conditioned on y as this variable is unknown when the model is used for prediction.

A. Supervised Variational Autoencoders
If our only objective were to obtain a generative model of x, we could simply parameterize q φ with a flexible neural network and maximize the Evidence Lower Bound Objective (ELBO) used in VAEs. However, focusing exclusively on the generative model often yields poor predictions of y [28]. This motivates the inclusion of a supervised loss, forming the SVAE objective where λ functions as a tuning parameter controlling the relative weight of the supervised loss. This objective can be approximated by empirical risk minimization of the observed data and estimated via gradient descent, as is commonly done with VAEs. For simplicity, we omit E p(x) in further derivations so our losses are defined for a single data sample. We note that in practice we would use an empirical risk minimization formulation rather than an expectation. If λ = 0, then (1) reduces to VAE objective. In practice, λ is usually set to a fairly high value to emphasize prediction, as the variance of y is often substantially outweighed by the variance of x [28].
The SVAE approach may appear ideal: the encoder q φ simultaneously provides an accurate reconstruction of x and a latent representation predictive of y. However, the inclusion of the predictive loss induces a systematic deviation of the encoder from the posterior distribution of the generative model, which we refer to as the deviance. This can be seen by analyzing the stationary point associated with φ in Proposition 1. Proposition 1. The stationary points of (1) can be found using the reparameterization trick used by Kingma et al. [24]. This reparameterization expresses the random variable s ∼ q φ (s|x) as a transformation of a random variable dependant on the observed data x, denoted g φ ( , x). Under this transformation, the stationary point such that ∇ φ L φ,θ,ψ = 0 is The proof and the stationary points for θ and ψ, given φ, match a VAE and are provided in the Supplemental Section A.
In the case where λ = 0 the left hand side of (2) must also be 0. This represents the standard stationary point of a VAE learned without the supervised loss. In a VAE, the values φ minimize the divergence of the variational approximation and the true posterior defined by θ. Thus, the right hand side is a systematic deviation from the posterior induced in the latent representation by the predictive objective. This deviation ensures increased relevance of the latent space to the outcome y, as larger values of λ correspond to an increased emphasis on the predictive objective. As θ is dependent only on the generative loss, the strength of this deviation indicates the amount of information relevant to y that is not included in the purely generative model.

B. A Metric of Scientific Utility for Experimental Design
While scientific goals are highly variable and cannot be completely captured by a single metric, we now explain why this deviation reduces scientific utility in our application. The ultimate goal of our work is to experimentally validate the relationship between the covariates x and the outcome y. The generative model defines this relationship as which posits an hypothesis as to how changing x will impact y based on the parameters θ and ψ. Using such a model, one can design an experiment to test whether the relationship defined by the model is causal. For example, this approach has been used successfully to design neurostimulation protocols [37,6,7]. Certain applications may require real time estimation of the latent space and outcomes [37]. As described previously, a common fast approximation is to use a variational approximation to obtain a predictive distribution as Unfortunately, θ provides no direct insight on p φ (y|x). However, interpreting θ on how to influence p φ (y|x) is reasonable when p θ (y|x) is "close" to p φ (y|x), with this approximation degrading in quality as these distributions diverge. This divergence can be quantified via a variety of ways, and we say a model has high scientific utility if this divergence is small.
Unfortunately, this measure typically cannot be computed exactly without expensive MCMC methods. However, if we modify our definition of scientific utility slightly to use a KLdivergence on the latent space p(s|x) [24], we can obtain a reasonable metric that is straightforward to compute. This is reasonable as the chosen relationship between s and y is a Lipschitz function, ensuring that small posterior divergences will correspond to small predictive differences.
This measure can be estimated as follows. We replace the true posterior p θ (s|x) with the encoder learned by a VAE q φ * (s|x) holding θ fixed. By definition, the VAE encoder minimizes the KL from the true posterior out of all variational approximations. By extension, any approximation close to the variational posterior will also have a small KL from the true posterior. As both q φ * (s|x) and q φ (s|x) are analytic (usually Gaussian), this divergence is straightforward to compute using an analytic formula.
In summary, we measure scientific utility as the KL between the supervised encoder and a "refit" VAE encoder holding the decoder and classifier parameters fixed. Small divergences suggest that the parameters θ provide reasonable insight on how to modify x to induce changes in y. Furthermore, the systematic discrepancy that SVAEs introduce demonstrated in the previous section inherently reduce this utility.

C. A Bilevel Objective to Maintain a Posterior Approximation
As described previously, an SVAE learns a system that predicts accurately but can induce significant bias, decreasing utility in our applications. We mitigate this issue by "recoupling" the encoder of an SVAE to match a VAE encoder via a constraint while still including the predictive loss in the objective function. This combined objective is This objective corresponds to a bilevel optimization problem [38]. To match the terminology of the bilevel optimization literature, the higher-level objective, which we denote F (θ, ψ, φ), matches an SVAE while the lower-level objective, which we denote f (φ, θ) constrains the encoder to approximate the posterior. For notational convenience, we define F x = E q φ (s|x) log p θ (x, s) − log q φ (s|x)] and F y = E q φ (s|x) λ log p ψ (y|s)] as they represent the generative and supervised aspects of the higher-level objective respectively.
We assume that this is an opportunistic problem, signifying that the lower-level objective is not actively opposing the primary objective [39]. Furthermore, we assume that the lowerlevel objective is unconstrained, as careful reparameterization of neural network parameters can remove most constraints. Even with these assumptions, bilevel optimization problems are notoriously difficult to analyze [38]. One common approach is to reformulate the the bilevel optimization as a standard optimization problem. We can achieve this by replacing the lower-level objective with a standard constraint that the gradient ∇ φ f = 0, corresponding to replacing the lower-level objective with its Karush-Kuhn-Tucker (KKT) conditions [40]. These two objectives are equivalent provided that f is convex [41]. This assumption holds over φ in the simplified case of a supervised autoencoder with an affine encoder, a model that has been used previously in multiple applications [42,28]. However, this is unlikely to hold once this is reformulated as a variational autoencoder, and certainly will not hold once φ corresponds to a deep neural network. When convexity fails to hold, not only can suboptimal stationary points be discovered, the global optimum may fail to be a stationary point [43].
We will continue with this substitution even when it violates convexity, which we believe is reasonable approach in practice to build intuition. First, recent work has shown that the function g φ is convex around the optimum in the function space in the infinite width limit [44,45]. The functional form is our primary interest, and the specific parameters of the encoder are irrelevant in our applications. Remarkably, in the finite width regime overparameterized neural networks also possess similar properties to convex problems [46]. Second, even local optima in neural networks yield state-of-the-art performance in many tasks. Work has shown that high-energy stationary points are often saddle points [47], which can be largely avoided with modern gradient descent techniques [48]. Thus, the stationary point associated with this substitution is likely to be reasonably close to the global optimum with models that have been used with great success in a variety of practical tasks.
Ultimately, we rely on a commonly used procedure in machine learning. We first develop an approach using strong assumptions and then provide intuitive justifications for why the approach is reasonable even when the assumptions are not strictly true, and finally we validate these claims empirically with practical applications.
Having replaced the lower level optimization with a gradient constraint, our objective now is We now have a single-level optimization problem with a nonconvex equality constraint.

D. Derivation of the Learning Algorithm
We now derive an algorithm for maximizing (6). We will use stochastic gradient descent to perform inference, as such methods provide computational tractability under the largedata regimes we encounter. However, this descent must be performed moving along the manifold defined by the constraint. The step on ψ is identical to an SVAE, as ψ is not involved in the constraint. The step on φ is straightforward as well; gradient descent is performed with respect to the constraint until convergence. However, the step for θ is more challenging.
One intuitive view of the constraint is that it makes φ a function of θ; the generative parameters determine the best variational approximation to the posterior. If this were a standard function, we could solve for φ in terms of θ in the primal objective and bypass the constraint entirely. We know that this is not possible; the non-convexity of the lower-level objective makes the feasible set a set-valued function. However, we can use the Implicit Function Theorem [49] to make this substitution valid within a neighborhood of the current values of θ and ψ. Let and q is the number of parameters in φ. Our manifold can then be defined as the level set such that G(θ, φ) = 0. The Implict Function Theorem states that if the second derivative of G with respect to θ has a bounded inverse at a point (φ 0 , θ 0 ), then there exists an open neighborhood Θ 0 × Φ 0 wherein an unique differentiable function h : Θ → Φ exists such that G(θ, h(θ)) = G(θ, φ). This unconstrained objective is We can now obtain the gradients of θ with respect to this unconstrained objective. The total derivative is This expression can be simplified by noting that as φ is learned to convergence with respect to the constraint, we can make the simplifying assumption ∂F x /∂h = 0. Thus, our update on θ consists of a standard partial derivative with respect to the generative likelihood as well as a second-order term induced by the constraint. The gradients for ∇ φ F , ∇ ψ F , ∂F x /∂θ and ∂F y /∂h are straightforward, however ∂h/∂θ is not. While the Implicit Function Theorem provides an explicit means to evaluate this gradient, it involves computing the Hessian matrix. While this is possible in low dimensions, it is clearly infeasible when modern neural networks containing millions of parameters are used for the encoder.
Instead, we must find a means to approximate this secondorder term. This can be done using an optimization trick used in meta-learning and self-supervised learning [20]. While our work differs dramatically in objective, we share a common inference technique. This approximation replaces the Hessian ∇ θ ∇ φ with a discrete approximation ∇ θ δ φ . This can be efficiently implemented in modern packages by maintaining the computational graph but estimating the gradients of θ using the updated parameters φ estimated after learning to convergence.
In practice, we only use a single update of φ. As the steps on θ are small, this step φ maintains the constraint within a small tolerance. The parameters ψ are also updated with a single step, and θ is updated with a standard gradient descent step along with a second-order step to account for the dependence of φ on θ. The procedure is summarized in Algorithm 1.

E. Stationary Point Analysis
A natural question at this point is how the constraint on φ is incorporated, as each step in Algorithm 1 is unconstrained. We rely in the strong regularity properties of neural networks for this justification. At each step, the update on θ moves the current parameter values to a point in Θ × Φ × Ψ that violates the constraint. However, our update on φ to convergence returns the parameter values to the manifold defined by the constraint. While in practice this single update approximation has worked well, if necessary φ can be updated many times to more closely satisfy the constraint.
We now sketch a justification for convergence of our algorithm. We first note that the discrete approximation trick used in our algorithm has been shown to converge under mild regularity conditions in meta-learning [50]. Based on Initialize: Network parameters: φ, θ, ψ; learning rate: α, β; weights: λ.
Step 1: decoder update Step 2: encoder update Step 3: classifier update Step 4: second-order decoder update Step 5: Model parameters update the smoothness assumptions of the objective, we can ensure that the primal objective increases at each iteration by selecting a sufficiently small step size. As our objective is a constrained variational autoencoder, it automatically inherits the property of being a lower bound on the marginal log likelihood. This likelihood often corresponds to a simple generative model and consequently we can assume this likelihood is finite. As a bounded monotonically increasing sequence, the primal objective is guaranteed to converge. Furthermore, local convexity in the function space ensures that the model parameters, θ and φ, will be a local optimum. In summary, our algorithm will converge given a set of assumptions that have held empirically in a variety of applications, namely convexity in the function space and smoothness of the networks.
Assuming our given algorithm converges, we can now analyze properties of the resulting stationary point. In particular, we now provide a necessary condition the resulting stationary points. The procedure is largely identical to the stationary points for an SVAE, and hence are deferred to the supplementary material. The only difference is the second-order term stemming from the constraint must be considered with respect to θ. Proposition 2. The stationary points of (6) can be found using the same reparameterization trick in Proposition 1. The stationary point of φ matches a VAE and is However, θ is modified to induce a predictive posterior and has a stationary point of The stationary point form for ψ given φ matches an SVAE.
The stationary points of the SOS-VAE provide an interesting contrast with an SVAE. Our modified objective induces a predictive latent space by influencing the generative parameters, x and y denote an arbitrary sample. In (b), we mark the 4 learning steps from Algorithm 1. The major difference is that a second-order step is used on the supervised loss to update the generative model parameters rather than the encoder. In (c), the 5 update steps from Supplemental Algorithm 2. The major change from SOS-VAE is that a second encoder is used to relax the exact inference strategy and measure divergence. hence our colloquial description of supervising the decoder. As the stationary point of the encoder is solely with respect to the generative model, refitting an encoder with the VAE loss and a fixed θ will not meaningfully change the encoder, making sure we maintain utility for our applications. We demonstrate this claim empirically in Section V-A.

IV. RELAXING THE SOS-VAE FRAMEWORK
SVAEs can be viewed as a regularization technique on a flexible classifier [19]. In contrast, the SOS-VAE explicitly uses the generative model to induce a predictive latent space. If a need for predictive accuracy matches or exceeds the need for informative decoder parameters, it may be preferable to allow the encoder to deviate slightly from the posterior to obtain increased flexibility, as a slight relaxation can yield substantial gains in predictive accuracy with only a minimal decrease in scientific utility. This improved performance can be incredibly valuable when the model is used to track the auxiliary variable, such as in real-time stimulation experiments [51]. We thus extend SOS-VAE by relaxing the constraint on the encoder, which we refer to as the Second-Order Supervision Double VAE (SOS-DVAE). We show that this approach can also be used to model missing data in Appendix C. Such missingness commonly arises when datasets obtained from multiple experiments are combined to broaden conclusions and increase statistical power or do to common electrode failure.
The SOS-DVAE approach defines two distinct encoders, a generative encoder q φ1 (s|x) that exclusively approximates the generative posterior and a predictive encoder q φ2 (s|x) that includes the supervision task. The generative encoder is learned as in Section III-C to approximate the true generative posterior. However, the predictive encoder is learned to maintain high predictive ability but with regularization towards the generative encoder. One objective that allows for this tradeoff is Here, µ is a tuning parameter controlling the strength of regularizing φ 2 to yield proper posterior inference (as µ → ∞, φ 2 is forced to exactly match φ 1 ). The parameters of the generative encoder, φ 1 , are learned in the same way as φ in (6), meaning that it approximates the posterior of the generative model while inducing the generative model parameters to be useful for prediction. The parameters of the predictive encoder φ 2 are free to deviate in approximating the generative posterior in order to yield improved predictions, provided that q φ2 is close to the generative encoder as measured by a KL regularization term. The KL term also allows straightforward estimation of how much deviation there is from the generative model. We note that alternative distribution distances could be used if desired, and KL is one approach that provides straightforward gradients. This learning approach is visualized in Fig. 2c and follows the same logic as SOS-VAE (Fig. 2b).
Once again, we proceed using by replacing the constraint with its KKT conditions and rely on the previously-provided justifications for this substitution. Pseudo-code is provided in Supplemental Algorithm 2. We note that removing the second-order optimization step is empirically detrimental to the predictive accuracy and increases the deviation of the variational approximation. This importance of the second-order step contrasts with other work that uses this technique, which found that such a change was largely inconsequential [20].

V. EXPERIMENTAL RESULTS
We present empirical evaluations on the SOS-VAE and SOS-DVAE inference strategies compared to several baseline approaches: 1) VAE-refit -a sequential fitting strategy with a VAE approximating a generative model and then using a supervised network on the frozen latent representation (this matches the "cutting-the-feedback" strategy in statistics [22]), 2) SVAE -the standard SVAE technique, 3) SVAE-refit -an SVAE where the encoder is refit to optimize only the generative objective while the supervised component remains fixed after initial training, and 4) SDVAE -the SOS-DVAE approach without the second-order step (Supplemental Algorithm 3). These strategies are chosen to compare different aspects of scientific utility. The VAE-refit enforces that the encoder approximate the generative posterior distribution as well as possible, at the cost of a potentially less predictive latent space. The SVAE will yield good predictive and reconstructive performance but has a biased approximation of the posterior distribution. Comparing the performance of SVAE to that of SVAE-refit reveals this inherent bias, as predictive improvements often do not induce changes in the generative parameters. Finally, SDVAE illustrates the need for our second-order strategies.
We first demonstrate our methods on the MNIST dataset as a proof-of-concept to illustrate the issues associated with a biased approximation. We then demonstrate scientific applications on two neural recording datasets, one consisting of LFPs recorded in mice and the other consisting of Electroencephalography (EEG) measurements in humans.
In each application, we use KL to evaluate the discrepancy between the posteriors from the two encoders in the SOS-DVAE caused by the supervision. There are other commonly-used methods for comparing distribution similarity [52,53,54]. However, KL is chosen for two important reasons. First, the likelihood is nearly ubiquitous in statistics to evaluate model fit, and as this work is motivated by generative modeling the KL is a natural choice. More importantly, these models are learned to minimize this KL, and evaluating similarity using the KL makes it clear that the bias stems from the supervised objective rather than the similarity metric used for training. Fig. 3: T-SNE visualization of the learned latent representation space for MNIST dataset using each model type. Each class is represented by a different color. λ = 1e4, feature size of 28x28, learning rate of 1e-3, and a batch size of 128 were used in the SVAE model.
A Multilayer Perceptron (MLP) with a single hidden layer is used as the encoder for all experiments. We evaluate two types of generative decoders. First, we use a Non-negative Matrix Factorization (NMF), which has been frequently used as a model capable of uncovering latent networks from neural electrophysiological data [4,28,55,56]. Second, we use an MLP decoder with a single hidden layer to demonstrate the broad applicability of the methods. All models have a latent space dimension of 20. These parameter choices are made to match common settings where supervised generative models have been used [10,8]. Additional details on implementation and training are given in the Supplemental Section D. The code to reproduce these result is available on Github 1 .

A. Demonstrating the Effect of Posterior Bias
To demonstrate the consequences of bias in the variational approximation of a standard SVAE and the effectiveness of our proposed solution, we perform baseline comparisons in the familiar MNIST database of handwritten digits [57]. The model is trained and evaluated using 10-fold cross-validation on the 70,000 images. Full training parameters are in Supplemental  Table III.
We first visualize the effect of the bias in SVAE models on the distribution of the latent representations in Fig. 3. We find that the latent space of the SVAE model is changed to the point that there is no resemblance to the latent space of the fully generative VAE model. In contrast, the latent space for the SOS-VAE and SVAE-refit models closely resemble that of the VAE. The massive change in latent space from SVAE to SVAE-refit clearly explains why predictive performance is lost in the refit model Fig. 4 shows the effect of tuning λ for an SVAE. We use Peak Signal-to-Noise Ratio (PSNR) and Structural Similarity Index (SSIM) to measure the quality of the generated images relative to the original. When λ goes higher, the generated images are of lower quality when compared to the original. Along with visualizations of the latent space, this demonstrates that the generative structure in the latent space is lost in exchange for good predictive accuracy with high λ.
Next, we evaluate predictive performance of our methods compared to the baselines previously mentioned. Predictive performance is quantified via accuracy (ACC) and area under the ROC curve (AUC) ( Table I). Because MNIST is a multi-class classification problem, the AUC is averaged over all individual classification tasks. Unsurprisingly, among the models that use an MLP decoder there is little variation between our novel inference method and the predictive baselines. The "cutting-the-feedback" method had lower accuracy due to the latent space not focusing on prediction, but the classes are still largely distinguishable in this latent representation. However, refitting the encoder of the SVAE (the refit encoder approximates the true posterior of the SVAE generative model) resulted in a dramatic drop in predictive performance. In other words, the bias introduced when estimating the encoder is absolutely critical for maintaining predictive ability and this predictive information is not incorporated into the unbiased generative model.
The results of the models using an NMF decoder are equally interesting, yielding patterns that are dramatically different from the MLP decoder. Here, the SOS-VAE results in a dramatic drop in predictive ability relative to the SVAE. This drop is unsurprising due to the restrictive constraint imposed by shallow generative models that was noted previously [23]. As an NMF decoder is substantially less complex than an MLP, we would expect this constraint to become more apparent. However, the SOS-DVAE, by relaxing this stringent condition, is able to maintain a high predictive accuracy marginally superior to an SVAE. The novel method without the second-order step (SDVAE) is also able to maintain high predictive accuracy.
We now analyze the other aspect required for scientific utility, that the latent space estimated by the encoder approximates the posterior of the generative model. While it is difficult to inspect this posterior, we can visualize the reconstruction of samples to determine the quality of the variational approximation to the posterior defined by the decoder. An example of this is shown in Fig. 5 using an MLP decoder. The top left represents the original image, with the remaining images depicting the reconstructions provided by various models. The SOS-VAE strongly resembles the original image, which is unsurprising as it is an unbiased approximation to the true posterior. Below, we show the reconstruction of the SVAE, along with the reconstruction after the encoder is refit. Both methods provide reasonable reconstructions of the image, but the refit SVAE has substantially worse predictive performance. This is precisely why bias in the encoder has remained undetected; the SVAE provides reasonable reconstructions and excellent predictions. It is only when the generative parameters are used, either for causal manipulations or scientific interpretation, that the bias manifests itself. Fig. 5 visualizes that the results of the SOS-DVAE as we increasingly relax the requirement that the variational approximation remains unbiased (smaller values of µ). Each pair represents models learned for a particular strength, with the left representing the reconstruction of the generative encoder f φ1 while the right represents the reconstruction using the predictive encoder f φ2 . We can see that a strong emphasis on an unbiased encoder result in minor gains in predictive accuracy, small divergences between the two encoders, and similar reconstructions from both encoders. However, the predictive encoder f φ2 increasingly diverges from the generative encoder as this regularization shrinks, yielding increasingly poor reconstructions from the predictive encoder. However, this relaxation in alignment between the two encoders corresponds to increasingly improved predictions, emphasizing the utility of µ as a tuning parameter. In particular, it is worth noting that a substantial improvement in predictive ability can be obtained with minimal increase in bias as shown when µ = 20.
Finally, we can analyze the impact of the second-order update by performing a similar analysis without this step (SDVAE), with results shown in the bottom row. We can see that the predictive accuracy is worse almost uniformly over the entire range. Furthermore, while the bias in the encoder is smaller relative to the SOS-DVAE initially, it quickly becomes substantially larger. This results in almost unrecognizable reconstructions using the predictive encoder at the weakest regularization strengths. From this, we conclude that the novel second-order optimization step developed in this work substantially improves both predictive and generative performance.

B. Comparison with Other Semi-supervised Generative Models
We are also interested in how realistic the reconstructions provided by our unbiased and relaxed approximations are as compared to alternative deep learning techniques with much more complex generative networks. This is evaluated on the MNIST dataset for the same reason as stated above; unlike LFPs, MNIST is easy to visualize and intuitively evaluate. While there are metrics for evaluating the quality of images, such as the Frechet Inception Distance (FID) [58], these methods have been developed for a substantially different task and are not reliable for MNIST. Instead, we visually evaluate quality, an approach commonly used to compare generative adversarial networks [59].
We provide a visualization of several MNIST reconstructions in Supplemental Fig. 11. We compare the reconstructions using a standard VAE [24], its semi-supervised extension [32], variational flows [33], auxiliary deep generative models [60], semi-supervised deep generative models with smooth-ELBO [61], and our SOS-VAE and SOS-DVAE. Our methods clearly provide more realistic reconstructions. The other four models yield blurry and diffuse reconstructions to the point that some reconstructions are illegible. Our methods provided reconstructions with distinct lines and no visibly missing sections. It is particularly interesting that our methods provided sharper-looking reconstructions as compared to the VAE with an exclusively generative focus. However, this development is not necessarily surprising, given the extensive work showing that auxiliary tasks can improve performance as relevant regularization [62].
Together, the above analyses strongly support multiple conclusions. First, the encoder bias identified in Section III is a serious concern that commonly affects supervised generative models. Second, our developed inference methodology addresses this bias, either completely removing it with the SOS-VAE or a controlled reduction via the SOS-DVAE. And third, the reconstructions provided by our methods are likely more realistic and a better model of the generative dynamics.

1) Decoding Emotional Affect from Electroencephalography
Recordings: We then apply the inference methods described above to the publicly available SEED electroencephalography (EEG) dataset [63,64]. This dataset includes 15 subjects recorded while watching movie clips designated with a negative/neutral/positive emotion label. We analyze signals from a curated subset of 19 of the original 62 electrodes (see Supplemental Fig. 7). We split the signals into non-overlapping 1 s time windows. For each window, we calculate the spectral power and coherence for 5 different frequency bands: 1-4 Hz, 5-8 Hz, 9-12 Hz, 13-30 Hz, and 31-50 Hz, corresponding to the Delta, Theta, Alpha, Beta, and Low Gamma bands typically used in EEG analysis, respectively. Models are trained to reconstruct the power and coherence values and to classify the emotion label associated with each window. We use a leave-one-participant-out cross-validation scheme to select hyperparameters [65]. Table I reports the performance of each inferred model. We see that with the MLP decoder, the SOS-VAE gives much better predictive performance than the SVAE-refit model, indicating that the predictive performance of the SVAE model is driven by features in the latent space that are not relevant to the generative model. The SOS-DVAE sacrifices relatively little performance compared to the SVAE while we again find that SOS-DVAE can get higher predictive performance at smaller KLs compared to SDVAE (see Supplemental Figs. 8c and 8d).
A major goal of this work is to improve the utility of generative models for drawing scientific conclusions. We have shown that our proposed modifications to the SVAE produce generative models that more accurately explain variables of interest through a supervised task. For neural datasets such as this one, we can use the decoder parameters to draw scientific conclusions. By using a non-negative matrix factorization model as our decoder, we can interpret the learned factors of the decoder as network factors of the electrical functional connectome (electome) [4,28,6]. To demonstrate the practical usefulness of this approach, we visualize the latent electome factor with the largest weight in the logistic regression classifier in Fig. 6a. This factor represents a network of brain regions defined by the power and coherence signatures we expect the network to produce. This network is defined by nearly fullbrain synchrony in the alpha band, whereas the other bands have localized coherence between PZ, P4, P8, O1, and O2.
2) Decoding Behavioral Context from Local Field Potential Recordings: We next apply the proposed inference methods to a dataset of Local Field Potentials (LFPs) recorded from 11  data with all and missing channels, respectively (see Appendix C). The outermost set of labels are abbreviated names for each of the brain regions present in the recordings. The inner labels iterate over the modeled frequencies. Colored segments along the outer "wheel" of the image indicate that the factor represents signal power within that region and frequency range. Colored "spokes" between regions indicate that the factor represents coherence between those two regions at the associated frequency. different brain regions across 26 mice [4,5]. Full brain region names are given in Supplemental Table IV, with abbreviated names referenced in Fig. 6. Each mouse is recorded in three different behavioral contexts, which are thought to induce low, medium, and high levels of stress respectively. As with the EEG dataset, recordings are divided into 1 second non-overlapping time windows. Spectral power within each brain region and coherence between brain regions are calculated at frequencies from 1 Hz to 56 Hz in 1 Hz increments for each time window [28]. We use a multinomial logistic regression for supervised classification of the behavioral context and evaluate the model performance via a 5-fold cross-validation over mice. We report results for two different types of decoders, MLP and NMF, as in the previous sections.
For each inferred model we evaluate performance on the supervised task as well as the KL between posterior distributions if possible (see Table II). As expected, SVAE and SOS-DVAE display comparable classification performance while SVAE-refit has near random performance, indicating that the performance of the SVAE model is not associated with the generative aspects of the model. Supplemental Figs. 8e and 8f visualize the tradeoff between predictive performance and KL for SDVAE and SOS-DVAE. The SOS-DVAE is more predictive at the same KL, demonstrating that it is guiding the generative model towards the classification goal. Fig. 6b shows the electome network with the largest weight in the classifier. The network is positively associated with the medium stress behavioral context, and demonstrates significant increases in connectivity between several brain regions in the low Gamma band of 30-50 Hz and synchrony between an overlapping group of regions around 12Hz.

VI. DISCUSSION AND CONCLUSION
Generative latent variable models have great utility in scientific and clinical trial analysis to improve hypothesis generation and experimental design [66,67]. This scientific and clinical utility often depends on two goals, obtaining an accurate representation of the data and ensuring that the latent representation is predictive of an auxiliary task. A commonly used inference method, the SVAE, has been previously used to achieve both of these objectives. However, this objective function induces a previously-undetected bias in the encoder that hinders scientific utility. We have shown, both on synthetic and real data, that the theoretical bias in SVAEs corresponds to a detrimental impact on learned representations in practice. To address this, we developed a novel inference technique that allows for supervision of an auxiliary task while maintaining a generative representation. We demonstrate empirically that our novel inference technique achieves both scientific objectives without bias, and demonstrated the efficacy of the proposed methodology in relevant neuroscience applications. Further-more, we have provided two relevant extensions to our methods that address critical needs in the neuroscience community.
We see two possible limitations of the proposed method. First, our models are designed for interpretable generative models commonly used in scientific analysis. This objective was not typically present in related work, which largely used the generative model as a regularization technique. When pursuing purely a predictive problem, a user can choose to use a much more complex yet difficult-to-interpret deep generative model [32,60,68,20,21], making comparisons with the simpler, interpretable generative model learned via an SOS-VAE irrelevant. Second, our contributions are applicable to models that include both a generative and a supervised component in the class of SAEs, whereas many supervised models have alternative inference strategies, such as [4].
In conclusion, our developed inference techniques are highly relevant to scientific fields such as neuroscience that use latent variable models to design experiments and discover novel relationships in high-dimensional data. These techniques improve the scientific utility of these latent variable models by incorporating predictive information while maintaining a clear understanding of how manipulations of observed covariates will result in changes in the latent space. In the future, we will continue applying this to real-world scientific problems and work to build greater integration with standard black-box variational inference tools [69,70,71].
The gradients of θ and ψ are straightforward. However, the gradients of φ are more difficult, as the ELBO expectation is taken with respect to a random sample from q φ (s|x). However, if we express the random variable s ∼ q φ (s|x) as a transformation of random variable given x and φ, g φ ( , x), the distribution of will be independent of x and φ. The gradients for θ are while the gradients for ψ are (14) Finally, the gradients for φ are (15) Thus, the fixed points for θ and ψ are and the fixed point for φ is 2) Second-Order Supervision: The second-order supervision objective recouples the encoder and decoder by forcing the encoder to approximate solely the generative posterior. Mathematically this objective is The gradients for ψ follow almost identically to the SVAE as The gradients taken with respect to φ are done wrt the constraint as Finally, the gradients with respect to θ are (21) This implies that the fixed points of the SOS-VAE for φ and ψ are while the fixed point for θ is 3) Gradient Update Derivations: As mentioned previously, ∇ θ log p ψ (y|g φ ( , x)) is strange given that θ and φ are both variables in common implementations. However, the constraint induces dependence on φ. To evaluate the gradient on θ, note that the total derivative estimate for θ is (24) The first term is trivial to implement in modern software packages. When the constraint on φ is approximately satisfied ∂φ/∂θ is small and can be ignored. While this term is not close to zero when the networks are initialized, it is still small comparitively and can be reasonably ignored. The third term can be approximately evaluated using the second-order trick described in [20]. obtained via each encoder reconstructs the the observed data while predicting the unobserved covariates. Dependence between the encoders is obtained through a shared θ between all experiments. Given a test sample x * with a new set of observed regions, the posterior can be approximated by learning a new encoder with parameters φ * using the objective (26) If the training encoders approximate exclusively a generative model (a VAE), this objective will be consistent. However, in an SVAE formulation this approach will approximate the "refit" encoder defined by the generative model. If the predictive capability of the SVAE is largely obtained via the deviation, this approximation will lose most of the predictive ability. Our methodology will be resilient to such problems, as the supervision is induced in the decoder parameters θ, which are used to train the new encoder.

D. Implementation Details
For all models, the encoder uses a single hidden layer MLP with 512 nodes. Two decoders are investigated: (i) a single layer MLP with 512 nodes, and (ii) a Non-negative Matrix Factorization (NMF) decoder. In the NMF decoder, the latent space is mapped to non-negative values through a softplus nonlinearity, and then a non-negative linear mapping is learned to project to the outputs.
All models are implemented using Pytorch. The experiments are run on a cluster with a Red Hat Enterprise Linux 7 operating system and a range of Nvidia GPUs, including TitanXPs and RTX2080Tis. Training parameters used for the experiments for each dataset are summarized in Table III. All hyperparameters were selected in the same way. For a standard SVAE and SOS-VAE, the λ, learning rate, and batch size are important hyperparameters. In the SDVAE and SOS-DVAE models an extra hyperparameter is a scaling factor µ balancing the Kullback-Liebler (KL) divergence with the other terms of the objective function. These parameters are tuned in 3 steps: 1) a sparse search with 13 values in a bigger range to locate smaller ranges that provide better prediction performance with 10-fold cross-validation, 2) shrink the searching range and add new values in the new searching range that obtains better performance, 3) run multiple times and chose the parameters that offers the best performance. For λ, the initial search range is [1e -8, 1e8]. For µ, the initial search ranging [1e-5,5] to investigate the trade off between generation and inference for each dataset using MLP decoder (see Figure 8). Note that we used a single hidden layer MLP to demonstrate the broad applicability of the methods, and that the modification to a more complex architecture is straightforward to implement.

1) MNIST:
The MNIST contains 60,000 training images and 10,000 test images, which we concatenate together into one dataset for a 10-fold cross-validation with random splits. No data augmentation is applied to the dataset, since our goal is to compare the performance among different models, not to pursue the best prediction.
2) Local Field Potential (LFP): The Local Filed Potentials (LFPs) are recorded from 11 different brain regions (see Table  IV) for each mouse [5]. Recordings are split into 1s intervals, each with an associated genotype and condition label. We estimate the spectral power features for each time interval using Welch's method [74] and mean squared coherence between pairs of brain regions [75] as measures of frequency-resolved synchrony within and between regions, respectively. These features are calculated at 1 Hz intervals from 1 Hz to 56 Hz, yielding a 3696 dimensional observation space.
The model for the outcome given the latent factors naturally lends itself to multinomial regression for the prediction of low-, medium-, and high-stress contexts corresponding to experimental conditions of home cage, open field, and tail suspension respectively. The NMF decoder naturally lends itself as a biologically interpretable model for neural electrophysiology [28]. It views each observation as a positive sum of nonnegative features, which matches the biological assumption that no network of neural activity can be negatively activated or be associated with with negative power or coherence. We also use an MLP decoder to demonstrate the broad applicability of these methods.

3) Electroencephalography(EEG):
For the electroencephalography (EEG) recordings of the SEED dataset [63,64], we use a subset of 19 electrodes (as shown in Figure 7) that approximate the whole head to simplify visualizations, but it would be a straightforward extension to use the full set of 62 electrodes. Power and coherence features are calculated for the 19 electrodes in 1 second time windows at the following frequency bands: 1-4, 5-8, 9-12, 13-30, and 31-50 Hz, corresponding to the Delta, Theta, Alpha, Beta, and Low Gamma, respectively. The resulting power and coherence features are flattened and concatenated to yield a 950 dimensional observation space. We use the MNE package (an open-source tool available at: https://mne.tools/stable/index.html) to extract the power and coherent features for this EEG dataset.

E. SOS-DVAE Improves Robustness to Missing Data
One of the appealing novelties of our method is that it addresses missing data issues that often arise when combining recordings from multiple scientific experiments. To demonstrate the concept of such application, we first generate a synthetic dataset by randomly removing 3 brain regions in each mouse in the LFP dataset (described in Section V-C2), and then we apply the extension of our proposed method to handling such missing data. Finally, we compare the trained models with that trained on full channels as shown in Figures 6b and 6c.
In practice, data may not be collected consistently, resulting in missing portions of the data. For example, LFP recordings such as the ones used above often suffer from overly noisy channels that prevent signal from being observed in one or more of the brain regions in the study. We simulate this scenario by randomly removing 3 brain regions in each mouse in the LFP dataset from Section V-C2, and we apply the missing data methodology from Section C.
The missing data scenario is much less harmful to the model learned via the SOS-DVAE as compared to an SDVAE. We show in Table V that the SOS-DVAE outperforms the SDVAE at decoding behavioral context in all scenarios. The single latent factor with the largest weight in the logistic regression classifier is shown in Fig. 6c next to the equivalent factor from a model with all channels present (Fig. 6b). We see that these factors share many features in common.

F. Additional results and comparison with other methods
As can be seen from Algorithm 2, µ weights the KLdivergence between the inference encoder q φ1 (·) and the classification q φ2 (·) encoder, with higher values for µ emphasizing that KL-divergence term over the classification loss. We tune µ to investigate the trade-off between prediction performance and the fidelity of generative model. We report some of these results in the main paper, and included additional supplemental results here. Figures 8b, 8d, and 8f show the prediction performance against the KL-divergence on all three datasets using MLP decoders. We can see that the proposed SOS-DVAE obtains higher predictive performance (ACC or AUC) for the same level of KL-divergence. As visualized in Figures 9 for the LFP dataset and Figure 10 for the SEED dataset, a smaller KLdivergence (greater µ) in SOD-DVAE and SDVAE means the φ 2 encoders are more influenced by the generative model and thus are able to reconstruct the input sample much better. In both figures, the trained models with a Multilayer Perceptron (MLP) decoder are used to reconstruct an arbitrary input sample from the hold-out data. Note that the SVAE and SVAE-refit only has a single encoder thus they do not have a KL-divergence value. The SVAE, SOS-DVAE, and SDVAE perform similarly in prediction. Similar to Figure 5 for the MNIST dataset, the SVAE-refit presents a huge decrease in prediction. Despite predicting less well, the refitted model is able to generate the input features better than that of SVAE.
As described in Section V-B in the main manuscript, Figure  9 compares the generative ability to reconstruct an input image using VAE [24], its semi-supervised extension [32], variational flows [33], auxiliary deep generative models [60], semi-supervised deep generative models with smooth-ELBO [61], and the proposed SOS-VAE and SOS-DVAE. The methods in comparison are trained on the train set of MNIST using the default parameters in the publicly released code. The examples shown in Figure 9 are arbitrary unseen images from the test set of MNIST.   9: Visualization of reconstructed power features on LFP dataset using MLP decoder. In the bottom left, we show the results for an SVAE model before and after (SVAE-refit) updating the encoder to prioritize the generative model. The results of SOS-DVAE (above yellow arrow) and SDVAE (above blue arrow) models are shown with 3 different KL-divergences (controlled by µ as detailed in Section IV and in Algorithm 2). For each µ, the left and right pictures are reconstructed from the inference encoder q φ 1 (·) and the classification q φ 2 (·) encoder, respectively. Bigger µ (smaller KL-divergence) yields similar reconstruction from the two encoders, vice versa. The numbers under the picture are KL-divergence in nats and prediction AUC. The models used in this figure is trained on a single split of training set, and the reconstruction is for a random input in the hold-out test set. Fig. 10: Visualization of reconstructed power features on SEED dataset using an MLP decoder. In the bottom left, we show the results for an SVAE model before and after (SVAE-refit) updating the encoder to prioritize the generative model. The results of SOS-DVAE (grouped by the yellow arrow) and SDVAE (grouped by the blue arrow) models are shown with 3 different KL-divergences (controlled by µ in Algorithm 2). For each µ, the left and right pictures are reconstructed from the inference encoder q φ 1 (·) and the classification q φ 2 (·) encoder, respectively. Bigger µ (smaller KL-divergence) yields similar reconstruction from the two encoders, vice versa. The numbers under the picture are KL-divergence in nats and prediction accuracy in percentage. The models used in this figure is trained on a single split of training set, and the reconstruction is for a random input in the hold-out test set. Fig. 11: Demonstration of the promising generative ability of our models. The last 2 columns are from our proposed SOS-DVAE model, with the left column from the inference encoder, and the right column from the classification encoder. The KL-divergence between the 2 encoders is 0.5.