Knowledge-Based Regularization in Generative Modeling

TL;DR

This paper introduces a regularization technique that integrates prior knowledge of feature relations into general-purpose generative models, improving their learning process without altering existing architectures.

Contribution

It proposes a flexible regularizer that encodes feature relation knowledge into various generative models using standard backpropagation, enhancing their ability to incorporate domain expertise.

Findings

Effective across multiple datasets and models

Compatible with VAEs and GANs

Improves generative modeling with prior knowledge

Abstract

Prior domain knowledge can greatly help to learn generative models. However, it is often too costly to hard-code prior knowledge as a specific model architecture, so we often have to use general-purpose models. In this paper, we propose a method to incorporate prior knowledge of feature relations into the learning of general-purpose generative models. To this end, we formulate a regularizer that makes the marginals of a generative model to follow prescribed relative dependence of features. It can be incorporated into off-the-shelf learning methods of many generative models, including variational autoencoders and generative adversarial networks, as its gradients can be computed using standard backpropagation techniques. We show the effectiveness of the proposed method with experiments on multiple types of datasets and generative models.

Tables

Table 1. Table 1: Test mean cross-entropy of VAEs on ccMNIST . Only the case of dim ( MLP ) = 2 10 dimension MLP superscript 2 10 \dim(\text{MLP})=2^{10} is reported, and the other cases are similar.

$dim\left(𝒛\right)$	L2 only	out-layer	proposed	dedicated decoder
25	$332$ ($3.4$)	$331$ ($3.6$)	$\mathrm{𝟑𝟐𝟏}$ ($2.4$)	$\mathrm{𝟑𝟐𝟓}$ ($1.5$)
50	$326$ ($3.3$)	$325$ ($3.1$)	$\mathrm{𝟑𝟏𝟐}$ ($2.6$)	$\mathrm{𝟑𝟎𝟏}$ ($4.3$)
100	$332$ ($2.9$)	$329$ ($2.1$)	$\mathrm{𝟑𝟏𝟕}$ ($2.0$)	$\mathrm{𝟑𝟎𝟓}$ ($5.9$)
200	$333$ ($4.4$)	$333$ ($3.2$)	$\mathrm{𝟑𝟐𝟐}$ ($2.4$)	$\mathrm{𝟑𝟎𝟓}$ ($4.2$)

Table 2. Table 2: Test mean per-feature reconstruction errors by VAEs on the three datasets. Here only the cases of dim ( MLP ) = 2 5 , 2 6 , 2 7 dimension MLP superscript 2 5 superscript 2 6 superscript 2 7 \dim(\text{MLP})=2^{5},2^{6},2^{7} (resp. for the three datasets) are reported. Other cases are similar.

dataset	$dim\left(𝒛\right)$	L2 only	out-layer	proposed
Toy	4	$0.76$ ($.15$)	$0.41$ ($.08$)	$0.41$ ($.07$)
Plant	4	$8.73$ ($.17$)	—	$8.67$ ($.10$)
	7	$8.21$ ($.15$)	—	$8.06$ ($.16$)
	11	$8.11$ ($.23$)	—	$7.95$ ($.17$)
Solar	4	$4.51$ ($.96$)	$4.32$ ($.46$)	$2.57$ ($.33$)
	13	$3.22$ ($.54$)	$3.15$ ($.47$)	$1.97$ ($.16$)
	54	$2.73$ ($.43$)	$2.71$ ($.49$)	$1.90$ ($.16$)
Significant diff. from “L2 only” with  italic:  or  bold: .

Table 3. Table 3: Sets of candidates for hyperparameter λ 𝜆 \lambda , from which λ 𝜆 \lambda was chosen by validation performance.

Model	Dataset	Candidates for $\lambda$
FA	ccMNIST	$5\times {10}^3$, ${10}^4$, or $2\times {10}^4$
	Plant	$2.5\times {10}^4$, $5\times {10}^4$, or ${10}^5$
	Solar	$5\times {10}^4$, ${10}^5$, or $2\times {10}^5$
VAE	Toy	$250$, $500$, or $1000$
	ccMNIST	$5\times {10}^3$, ${10}^4$, or $2\times {10}^4$
	Plant	$2.5\times {10}^4$, $5\times {10}^4$, or ${10}^5$
	Solar	$2.5\times {10}^5$, $5\times {10}^5$, or ${10}^6$

Table 4. Table 7: Test mean per-feature reconstruction errors by FA models. Averages (standard deviations) over 10 random trials are shown.

Setting		Results
Dataset	$dim\left(𝒛\right)$	L2 only	L2 + out-layer	L2 + proposed
ccMNIST	$200$	$1.970\times {10}^{-5}$ ($3\times {10}^{-9}$)	$1.970\times {10}^{-5}$ ($3\times {10}^{-9}$)	$1.958\times {10}^{-5}$ ($1\times {10}^{-8}$) ∗∗∗
	$800$	$8.150\times {10}^{-6}$ ($2\times {10}^{-9}$)	$8.148\times {10}^{-6}$ ($3\times {10}^{-9}$)	$7.828\times {10}^{-6}$ ($2\times {10}^{-8}$) ∗∗∗
Plant	$7$	$6.318$ ($0.021$)	—	$6.323$ ($0.019$) ${}^{}$
	$11$	$4.949$ ($0.039$)	—	$4.936$ ($0.025$) ${}^{}$
Solar	$13$	$1.232$ ($0.006$)	$1.232$ ($0.006$)	$1.227$ ($0.005$) ${}^{\ast }$
	$54$	$0.842$ ($0.004$)	$0.841$ ($0.004$)	$0.831$ ($0.004$) ∗∗∗

Equations

$$\underset{\theta }{\mathrm{minimize}}L(\theta ).$$

$$\text{HSIC}(\mathcal{F},\mathcal{G},p_{xy})=\mathbb{E}_{x,x^{\prime},y,y^{\prime}}[k(x,x^{\prime})l(y,y^{\prime})]\\ +\mathbb{E}_{x,x^{\prime}}[k(x,x^{\prime})]\mathbb{E}_{yy^{\prime}}[l(y,y^{\prime})]\\ -2\mathbb{E}_{x,y}\big{[}\mathbb{E}_{x^{\prime}}[k(x,x^{\prime})]\mathbb{E}_{y^{\prime}}[l(y,y^{\prime})]\big{]}\geq 0.$$

$$H_0:{\rho }_{x,y,z}\le 0\text{and}{H}_1:{\rho }_{x,y,z}>0,$$

$${\hat{\rho }}_{x,y,z}:={\widehat{\text{HSIC}}}_{x,y}(\mathcal{D})-{\widehat{\text{HSIC}}}_{x,z}(\mathcal{D}).$$

$$p\leq 1-\Phi\big{(}\hat{\rho}_{x,y,z}(\sigma^{2}_{xy}+\sigma^{2}_{xz}-2\sigma_{xyxz})^{-1/2}\big{)},$$

$$\mathcal{K}:=\big{\{}(J_{s}^{\text{ref}},J_{s}^{+},J_{s}^{-})\mid s=1,\dots,|\mathcal{K}|\big{\}},$$

$$\underset{\theta }{\mathrm{minimize}}L(\theta )+\lambda R_{\mathcal{K}}(\theta ),$$

$$\begin{array}{c}{R}_{\mathcal{K}}(\theta ):=\frac{1}{|\mathcal{K}|}\sum _{s=1}^{|\mathcal{K}|}{R}_{\mathcal{K},s}(\theta ),\\ R_{\mathcal{K},s}(\theta ):=\max (0,{\nu }_{\alpha }-{\hat{\rho }}_s(\theta )).\end{array}$$

$${\hat{\rho }}_s(\theta ):={\widehat{\text{HSIC}}}_{x_{J_s^{\text{ref}}},x_{J_s^{+}}}(\hat{\mathcal{D}})-{\widehat{\text{HSIC}}}_{x_{J_s^{\text{ref}}},x_{J_s^{-}}}(\hat{\mathcal{D}}),$$

$${\rho }_s:=\text{HSIC}({\mathcal{F}}_{J_s^{\text{ref}}},{\mathcal{F}}_{J_s^{+}})-\text{HSIC}({\mathcal{F}}_{J_s^{\text{ref}}},{\mathcal{F}}_{J_s^{-}}),$$

$$H_{s,0}:{\rho }_s\le 0\text{and}{H}_{s,1}:{\rho }_s>0,$$

$$p_s\le 1-\Phi ({\nu }_{\alpha }/{\tau }_s),\exists {\tau }_s>0.$$

$${\tau }_s=({\sigma }_{J_s^{\text{ref}}{J}_s^{+}}^2+{\sigma }_{J_s^{\text{ref}}{J}_s^{-}}^2-2{\sigma }_{J_s^{\text{ref}}{J}_s^{+}{J}_s^{\text{ref}}{J}_s^{-}}{)}^{1/2}>0.$$

$$\mathcal{K}_{\text{Toy}}=\big{\{}(J^{\text{ref}}_{1}=\{1\},\ J^{+}_{1}=\{2\},\ J^{-}_{1}=\{3\})\big{\}}.$$

$$\displaystyle\mathcal{K}_{\text{ccMNIST}}=\big{\{}($$

$$\displaystyle J^{-}_{1}=\{1569,...,2352\})\big{\}}.$$

$$(J_1^{\text{ref}}=J_{\text{compressor}},J_1^{+}=J_{\text{reactor}},J_1^{-}=J_{\text{stripper}}),\dots$$

$$R_{{\mathcal{K}}_{\text{toy}}}^{\text{out-layer}}:=\max (0,ϵ+\Vert {\mathbf{w}}_1-{\mathbf{w}}_2{\Vert }_2^2-\Vert {\mathbf{w}}_1-{\mathbf{w}}_3{\Vert }_2^2),$$

$$\displaystyle\mathcal{K}_{\text{Solar}}=\big{\{}$$

$$(J_{\text{Separator}},J_{\text{Stripper}},J_{\text{Feed2}}),(J_{\text{Separator}},J_{\text{Compressor}},J_{\text{Feed1}}),(J_{\text{Separator}},J_{\text{Compressor}},J_{\text{Feed2}}),$$

$$(J_{\text{Compressor}},J_{\text{Reactor}},J_{\text{Feed1}}),(J_{\text{Compressor}},J_{\text{Reactor}},J_{\text{Feed2}}),(J_{\text{Compressor}},J_{\text{Reactor}},J_{\text{Stripper}}),$$

$$\displaystyle~{}(J_{\text{Stripper}},\ J_{\text{Reactor}},\ J_{\text{Feed1}}),~{}~{}(J_{\text{Stripper}},\ J_{\text{Reactor}},\ J_{\text{Purge}}),~{}~{}(J_{\text{Stripper}},\ J_{\text{Reactor}},\ J_{\text{Compressor}})~{}\big{\}},$$

$$J_{\text{Feed1}}=\{1,2,3\},J_{\text{Feed2}}=\{4\},J_{\text{Reactor}}=\{5,6,7,8,9,21\},J_{\text{Separator}}=\{11,12,13,14,22\},$$

$$J_{\text{Stripper}}=\{15,16,17,18,19\},J_{\text{Compressor}}=\{20\},\text{and}{J}_{\text{Purge}}=\{10\}.$$

$$\mathbf{x}\approx \mathbf{W}\mathbf{z}+\mathbf{\mu },$$

$$\tilde{\mathbf{z}}={\mathbf{W}}^{\mathsf{T}}(\mathbf{x}-\mathbf{\mu }).$$

$$\Vert \mathbf{x}-\mathbf{W}{\mathbf{W}}^{\mathsf{T}}(\mathbf{x}-\mathbf{\mu })-\mathbf{\mu }{\Vert }_2^2.$$

$$\begin{gathered} \text{ELBO}={\mathbb{E}}_{q_{{\varphi }_{\text{enc}}}(\mathbf{z}\mid \mathbf{x})}[\log p_{{\theta }_{\text{dec}}}(\mathbf{x}\mid \mathbf{z})] \\ -\mathrm{KL}(q_{{\varphi }_{\text{enc}}}(\mathbf{z}\mid \mathbf{x})\Vert p_{{\theta }_{\text{dec}}}(\mathbf{x}\mid \mathbf{z})), \end{gathered}$$

$$\log D_{\phi_{D}}(\bm{x})+\log\big{(}1-D_{\phi_{D}}\big{(}G_{\theta_{G}}(\bm{z})\big{)}\big{)},$$

Full text

1 Introduction

Generative modeling plays a key role in many scientific and engineering applications. It has often been discussed with the notion of graphical models (Koller and Friedman, 2009), which are built upon knowledge of conditional independence of features. We have also seen the advances in neural network techniques for generative modeling, such as variational autoencoders (VAEs) (Kingma and Welling, 2014) and generative adversarial networks (GANs) (Goodfellow et al., 2014).

For efficient generative modeling, a model should be designed following prior knowledge of target phenomena. If one knows the conditional independence of features, it can be encoded as a graphical model. If there are some insights on the physics of the phenomena, a model can be built based on a known form of differential equations. Otherwise, a special neural network architecture may be created in accordance with the knowledge. However, although some tools (e.g., Koller and Pfeffer (1997)) have been suggested, it is often labor-intensive and possibly even infeasible to design a generative model meticulously so that the prior knowledge specific to each problem instance is hard-coded in the model.

On the other hand, we may employ general-purpose generative models, such as kernel-based models and neural networks with common architectures (multi-layer perceptrons, convolutional nets, and recurrent nets, etc.). Meanwhile, a general-purpose model, if used as is, is less efficient than specially-designed ones in terms of sample complexity because of large hypothesis space. Therefore, we want to incorporate as much prior knowledge as possible into a model, somehow avoiding the direct model design.

A promising way to incorporate prior knowledge into general-purpose models is via regularization. For example, structured sparsity regularization (Huang et al., 2011) is known as a method for regularizing linear models with a rigorous theoretical background. In a related context, posterior regularization (Ganchev et al., 2010; Zhu et al., 2014) has been discussed as a methodology to impose expectation constraints on learned distributions. However, the applicability of the existing regularization methods is still limited in terms of target types of prior knowledge and base models.

In generative modeling, we often have prior knowledge of relationship between features. For instance, consider generative modeling of sensor data, which is useful for tasks such as control and fault detection. In many instances, we know to some extent how units and sensors in a plant are connected (Figure 1a), from which the dependence of a part of features (here, sensor readings) can be anticipated. In the example of Figure 1a, features $x_{1}$ and $x_{2}$ would be more or equally dependent than $x_{1}$ and $x_{3}$ because of external disturbances in the processes between the units. Another example is when we know the pairwise relationship of features (Figure 1b), which can be derived from some side information.

The point is that we rarely know the full data structure, which is necessary for building a graphical model. Instead, we know a part of the data structure or pairwise relationship of some features, and there remain unknown parts that should be modeled using general-purpose models. Most existing methods cannot deal with such partial knowledge straightforwardly. In this paper, we propose a regularizer to incorporate such knowledge into general-purpose generative models based on the idea that statistical dependence of features can be (partly) anticipated from the relationship of features. The use of this type of knowledge has been actively discussed in several contexts of machine learning, but there have been surprisingly few studies in the context of generative modeling.

The proposed regularizer is defined using a kernel-based criterion of dependence (Gretton et al., 2005), which is advantageous because its gradients can be computed using standard backpropagation without additional iterative procedures. Consequently, it can be incorporated into off-the-shelf gradient-based learning methods of many general-purpose generative models, such as latent variable models (e.g., factor analysis, topic models, etc.), VAEs, and GANs. We conducted experiments using multiple datasets and generative models, and the results showcase that a model regularized using prior knowledge of feature relations achieves better generalization. The proposed method can provide a trade-off between performance and workload for model design.

2 Background

In this section, we briefly review two technical building blocks: generative modeling and dependence criteria.

2.1 Generative Modeling

We use the term “learning generative models” in the sense that we are to learn $p_{\theta}(\bm{x})$ explicitly or implicitly. Here, $\bm{x}$ denotes the observed variable, and $\theta$ is the set of parameters to be estimated. Many popular generative models are built as Bayesian networks and Markov random fields Koller and Friedman (2009). Also, the advances in deep learning technique include VAEs (Kingma and Welling, 2014), GANs (e.g., (Goodfellow et al., 2014)), autoregressive models (e.g., van den Oord et al. (2016)), and normalizing flows (e.g., Dinh et al. (2018)).

The learning strategies for generative models are usually based on minimization of some loss function $L(\theta)$:

[TABLE]

A typical loss function is the negative log-likelihood or its approximation (e.g., ELBO in variational Bayes) (Koller and Friedman, 2009; Kingma and Welling, 2014). Another class of loss functions is those designed to perform the comparison of model- and data-distributions, which has been studied recently often in the context of learning implicit generative models that have no explicit expression of likelihood. For example, GANs (Goodfellow et al., 2014) are learned via a two-player game between a discriminator and a generator.

2.2 Dependence Criteria

Among several measures of statistical dependence of random variables, we adopt a kernel-based method, namely Hilbert–Schmidt independence criterion (HSIC) (Gretton et al., 2005). The advantage of HSIC is discussed later in Section 3.3. Below we review the basic concepts.

HSIC is defined and computed as follows (Gretton et al., 2005). Let $p_{xy}$ be a joint measure over $(\mathcal{X}\times\mathcal{Y},\Gamma\times\Lambda)$, where $\mathcal{X}$ and $\mathcal{Y}$ are separable spaces, $\Gamma$ and $\Lambda$ are Borel sets on $\mathcal{X}$ and $\mathcal{Y}$, and $(\mathcal{X},\Gamma)$ and $(\mathcal{Y},\Lambda)$ are furnished with probability measure $p_{x}$ and $p_{y}$, respectively. Given reproducing kernel Hilbert spaces (RKHSs) $\mathcal{F}$ and $\mathcal{G}$ on $\mathcal{X}$ and $\mathcal{Y}$, respectively, HSIC is defined as the squared Hilbert–Schmidt norm of a cross-covariance operator $C_{xy}$, i.e., $\text{HSIC}(\mathcal{F},\mathcal{G},p_{xy}):=\|C_{xy}\|_{\mathrm{HS}}^{2}$. When bounded kernels $k$ and $l$ are uniquely associated with the RKHSs, $\mathcal{F}$ and $\mathcal{G}$, respectively,

[TABLE]

HSIC works as a dependence measure because $\text{HSIC}=0$ if and only if two random variables are statistically independent (Theorem 4, Gretton et al. (2005)).

HSIC can be empirically estimated using a dataset $\mathcal{D}=\{(x_{1},y_{1}),\dots,(x_{m},y_{m})\}$ from (2). Gretton et al. (Gretton et al., 2005) presented a biased estimator with $O(m^{-1})$ bias, and Song et al. (Song et al., 2012) suggested an unbiased estimator. In what follows, we denote an empirical estimation of $\text{HSIC}(\mathcal{F},\mathcal{G},p_{xy})$ computed using a dataset $\mathcal{D}$ by $\widehat{\text{HSIC}}_{x,y}(\mathcal{D})$. While the computation of empirical HSICs requires $O(m^{2})$ operations, it may be sped up by methods such as random Fourier features (Zhang et al., 2018).

As later discussed in Section 3.2, we particularly exploit relative dependence of the features. Bounliphone et al. (Bounliphone et al., 2015) discussed the use of HSIC for a test of relative dependence. Here, we introduce $\mathcal{H}$ as a separable RKHS on another separable space $\mathcal{Z}$. Then, the test of relative dependence (i.e., which pair is more dependent, $(x,y)$ or $(x,z)$?) is formulated with the null and alternative hypotheses:

[TABLE]

where $\rho_{x,y,z}:=\text{HSIC}(\mathcal{F},\mathcal{G},p_{xy})-\text{HSIC}(\mathcal{F},\mathcal{H},p_{xz})$. For this test, they define a statistic

[TABLE]

From the asymptotic distribution of empirical HSIC, it is known (Bounliphone et al., 2015) that a conservative estimation of $p$-value of the test with $\hat{\rho}_{x,y,z}$ is obtained as

[TABLE]

where $\Phi$ is the CDF of the standard normal, and $\sigma$’s denote the standard deviations of the asymptotic distribution of HSIC. See (Bounliphone et al., 2015) for the definition of $\sigma^{2}_{xy}$, $\sigma^{2}_{xz}$, and $\sigma_{xyxz}$, which can also be estimated empirically. Here, (4) means that under the null hypothesis $H_{0}$, the probability that $\rho_{x,y,z}$ is greater than or equal to $\hat{\rho}_{x,y,z}$ is bounded by this value. In the proposed method, we exploit this fact to formulate a regularization term.

3 Proposed Method

First, we manifest the type of generative models to which the proposed method applies. Then, we define what we expect to have as prior knowledge. Finally, we present the proposed regularizer and give its interpretation.

3.1 Target Type of Generative Models

Our regularization method is agnostic of the original loss function of generative modeling and the parametrization, as long as the following (informal) conditions are satisfied:

Assumption 1**.**

Samples from $p_{\theta}(\bm{x})$, namely $\hat{\bm{x}}$, can be drawn with an admissible computational cost.

Assumption 2**.**

Gradients \nabla_{\theta}\mathbb{E}_{p_{\theta}(\bm{x})}\big{[}f(\bm{x})\big{]} can be (approximately) computed with an admissible computational cost.

While Assumption 1 is satisfied in most generative models, Assumption 2 is a little less obvious. We note that the efficient computation of \nabla_{\theta}\mathbb{E}_{p_{\theta}(\bm{x})}\big{[}f(\bm{x})\big{]} is often inherently easy or facilitated with techniques such as the reparameterization trick (Kingma and Welling, 2014) or its variants, and thus these assumptions are satisfied by many popular methods such as factor analysis, VAEs, and GANs.

3.2 Target Type of Prior Knowledge

We suppose that we know the (partial) relationship of features that can be encoded as plausible relative dependence of the features. Such knowledge is frequently available in many practices of generative modeling, and in fact, the use of such knowledge has also been considered in other contexts of machine learning (see Section 4.2). Below we introduce several motivating examples and then give a formal definition.

3.2.1 Motivating Examples

The first motivating example is when we know (a part of) the data-generating process. Suppose to learn a generative model of sensor data of an industrial plant. We usually know how units and sensors in the plant are connected, which is an important source of prior knowledge. In Figure 1a, suppose that Unit $i$ has internal state $z_{i}$ for $i=1,\dots,3$. Also, suppose there are relations $z_{2}=h_{12}(z_{1},\omega_{1})$ and $z_{3}=h_{2}(z_{2},\omega_{2})$, where $\omega$’s are random noises and $h$’s are functions of physical processes. Then, $(p(z_{1}),p(z_{2}))$ would be more statistically dependent than $(p(z_{1}),p(z_{3}))$ because of the presence of $\omega_{2}$. If the sensor readings, $\{x_{i}\}$, are determined by $x_{i}=g(z_{i})$ with observation functions $\{g_{i}\}$, then, $(p(x_{1}),p(x_{2}))$ would be more statistically dependent than $(p(x_{1}),p(x_{3}))$ analogously to $z$. Figure 1a is revisited in Example 1. This type of prior knowledge is often available also for physical and biological phenomena.

Another type of example is when we know pairwise similarity or dissimilarity of features (Figure 1b). For example, we may estimate similarities of distributed sensors from their locations. Also, we may anticipate similarities of words using ontology or word embeddings. Moreover, we may know the similarities of molecules from their descriptions in the chemical compound analysis. We can anticipate relative dependence from such information, i.e., directly similar feature pairs are more dependent than dissimilar feature pairs. Figure 1b is revisited in Example 1. Here, we do not have to know every pairwise relation; knowledge on some feature pairs is sufficient to anticipate the relative dependence.

3.2.2 Definition of Prior Knowledge

As the exact degree of feature dependence can hardly be described precisely from prior knowledge, we use the dependence of a pair of features relative to another pair. This idea is formally written as follows:

Definition 1** (Knowledge of feature dependence).**

Suppose that the observed random variable, $\bm{x}$, is a tuple of $d$ random variables, i.e., $\bm{x}=(x_{1},\dots,x_{d})$. Knowledge of feature dependence is described as a set of triples:

[TABLE]

where $J^{\text{ref}},J^{+},J^{-}\subset\{1,\dots,d\}$ are index sets. The semantics are as follows. Let $J=\{i_{1},\dots,i_{|J|}\}$, and let $\bm{x}_{J}$ be the subtuple of $\bm{x}$ by $J$, i.e., $\bm{x}_{J}=(x_{i_{1}},\dots,x_{i_{|J|}})$. Then, triple $(J_{s}^{\text{ref}},J_{s}^{+},J_{s}^{-})$ encodes the following piece of knowledge: “$\bm{x}_{J_{s}^{\text{ref}}}$ is more dependent on $\bm{x}_{J_{s}^{+}}$ than on $\bm{x}_{J_{s}^{-}}$.”

Example 1**.**

$\mathcal{K}=\{(\{1\},\{2\},\{3\})\}$* in Figure 1a and $\mathcal{K}=\{(\{1\},\{2\},\{3,4\})\}$ in Figure 1b.*

More examples of $\mathcal{K}$ are in Section 5.

3.3 Proposed Regularization Method

We define the proposed regularization method and give its interpretation as a probabilistic penalty method.

3.3.1 Definition of Regularizer

We want to force a generative model $p_{\theta}(\bm{x})$ to follow the relations encoded in $\mathcal{K}$. To this end, the order of HSIC between the marginals of $p_{\theta}(\bm{x})$ should be as consistent as possible to the relations in $\mathcal{K}$. More concretely, the HSIC of $(p_{\theta}(\bm{x}_{J^{\text{ref}},s}),p_{\theta}(\bm{x}_{J^{+},s}))$ should be larger than that of $(p_{\theta}(\bm{x}_{J^{\text{ref}},s}),p_{\theta}(\bm{x}_{J^{-},s}))$, for $s=1,\dots,|\mathcal{K}|$.

Because directly imposing constraints on the true HSIC of the marginals is intractable, we resort to a penalty method using empirical HSIC. Concretely, we add a term that penalizes the violation of the knowledge in $\mathcal{K}$ as follows.

Definition 2** (Knowledge-based regularization).**

Let $L(\theta)$ be the original loss function in (1). The learning problem regularized using $\mathcal{K}$ in (5) is posed as

[TABLE]

where $\lambda\geq 0$ is a regularization hyperparameter, and

[TABLE]

Here, $\nu_{\alpha}\geq 0$ is another hyperparameter. $\hat{\rho}_{s}(\theta)$ is the empirical estimation of the relative HSIC corresponding to the $s$-th triple in $\mathcal{K}$, that is,

[TABLE]

where $\hat{\mathcal{D}}$ denotes a set of samples drawn from $p_{\theta}(\bm{x})$.

3.3.2 Interpretation

By imposing the regularizer $R_{\mathcal{K}}$, we expect that the hypothesis space to be explored is made smaller implicitly, which would result in better generalization capability. This is also supported by the following interpretation of the regularizer.

While $R_{\mathcal{K}}(\theta)$ depends on the empirical HSIC values, making (each summand of) $R_{\mathcal{K}}(\theta)$ small can be interpreted as imposing “probabilistic” constraints on the order of the true HSIC values via a penalty method. Now let $\rho_{s}$ denote the true value of $\hat{\rho}_{s}$ (we omit argument $\theta$ for simplicity), i.e.,

[TABLE]

where $\mathcal{F}_{J}$ denotes a separable RKHS on the space of $\bm{x}_{J}$. Moreover, consider a test with hypotheses:

[TABLE]

for $s=1,\dots,|\mathcal{K}|$. On the test with $H_{s,0}$ and $H_{s,1}$, the following proposition holds:

Proposition 1**.**

If $R_{\mathcal{K},s}=0$ is achieved, then the null hypothesis $H_{s,0}$ can be rejected in favor of the alternative $H_{s,1}$ with $p$-value $p_{s}$ upper bounded by

[TABLE]

Proof.

From (4), we have $p_{s}\leq 1-\Phi(\hat{\rho}_{s}/\tau_{s})$, where

[TABLE]

Also, when $R_{\mathcal{K},s}=0$, we have $\hat{\rho}_{s}\geq\nu_{\alpha}$. Consequently, as $\Phi$ is monotonically increasing, we have (11). ∎

3.3.3 Choice of Hyperparameters

The proposed regularizer has two hyperparameters, $\lambda$ and $\nu_{\alpha}$. While $\lambda$ is a standard regularization parameter that balances the loss and the regularizer, $\nu_{\alpha}$ can be interpreted as follows. Now let $\alpha:=1-\Phi(\nu_{\alpha}/\tau_{s})$ (i.e., the right-hand side of (11)), which corresponds to the required significance of our probabilistic constraints. We can determine $0<\alpha<1$ based on the plausibility of the prior knowledge, such as $\alpha=0.05$ or $\alpha=0.1$ (a smaller $\alpha$ requires higher significance). Also, $\tau_{s}$ can be estimated as it is defined with the variances of the asymptotic distribution of HSIC (Bounliphone et al., 2015). Hence, as $\nu_{\alpha}=\tau_{s}\Phi^{-1}(1-\alpha)$, we can roughly determine the value of $\nu_{\alpha}$ that corresponds to a specific significance, $\alpha$. The need to tune $\nu_{\alpha}$ (or $\alpha$) remains, but the above interpretation is useful to determine the search range for $\nu_{\alpha}$.

3.3.4 Optimization Method

From Assumptions 1 and 2, the gradient of $\hat{\rho}_{s}$ with regard to $\theta$ can be computed with the backpropagation via the samples from $p_{\theta}(\bm{x})$, which are used to compute empirical HSIC. Hence, if the solution of the original optimization, (1), is obtained via a gradient-based method, the regularized version, (6), can also be solved using the same gradient-based method. The situation is especially simplified if the dataset is a set of $d$-dimensional vectors, i.e., $\mathcal{D}=\{\bm{x}_{i}\in\mathbb{R}^{d}\mid i=1,\dots,n\}$. In the optimization process, sample $\hat{\mathcal{D}}=\{\hat{\bm{x}}_{i}\in\mathbb{R}^{d}\mid i=1,\dots,m\}$ from $p_{\theta}(\bm{x})$ being learned, and use them for computing the the regularization term in (7) and their gradients. Then, incorporate them into the original gradient-based updates. These procedures are summarized in Algorithm 1.

The computation of $R_{\mathcal{K}}$ requires $O(m^{2}|\mathcal{K}|)$ operations in a naive implementation, where $m$ is the number of samples drawn from $p_{\theta}(\bm{x})$. This will be burdensome for a very large $|\mathcal{K}|$, which can be alleviated by carefully choosing $J^{\text{ref}}$’s so that the computed HSIC values can be reused for many times.

4 Related Work

4.1 Generative Modeling with Prior Knowledge

A perspective related to this work is model design based on prior knowledge. For example, the object-oriented Bayesian network language (Koller and Pfeffer, 1997) helps graphical model designs when the structures behind data can be described in an object-oriented way. Such tools are useful when we can prepare prior knowledge enough for building a full model. However, it is not apparent how to utilize them when knowledge is only partially given, which is often the case in practice. In contrast, our method can exploit prior knowledge even if only a part of the structures is known in advance. Another related perspective is the structure learning of Bayesian networks with constraints from prior knowledge (Cussens et al., 2017; Li and van Beek, 2018).

Posterior regularization (PR) (Ganchev et al., 2010; Zhu et al., 2014; Mei et al., 2014; Hu et al., 2018) is known as a framework for incorporating prior knowledge into probabilistic models. In fact, our method can be included in the most general class of PR, which was briefly mentioned in (Zhu et al., 2014). However, practically, existing work (Ganchev et al., 2010; Zhu et al., 2014; Mei et al., 2014; Hu et al., 2018) only considered limited cases of PR, where constraints were written in terms of a linear monomial of expectation with regard to the target distribution. In contrast, our method tries to fulfill the constraints on statistical dependence, which needs more complex expressions.

The work by Lopez et al. [2018] should be understood as a kind of technical complement of ours (and vice versa). While Lopez et al. [2018] regularize the amortized inference of VAEs to make independence of latent variables, we regularize a general generative model itself to make dependence of observed variables compatible with prior knowledge. In other words, the former regularizes an encoder, whereas the latter regularizes a decoder. Moreover, our method is more flexible than Lopez et al. (2018) in terms of applicable type of knowledge; they considered only independence, but our method can incorporate both dependence and independence. Also, while Lopez et al. [2018] discussed the method only for VAEs, our method applies to any generative models as long as the mild assumptions in Section 3.1 are satisfied.

4.2 Knowledge of Feature Dependence

In fact, the use of knowledge of feature relation has been discussed in different or more specific contexts. In natural language processing, feature similarity is often available as the similarity between words. Xie et al. (Xie et al., 2015) exploited correlations of words for topic modeling, and Liu et al. (Liu et al., 2015) incorporated semantic similarity of words in learning word embeddings. These are somewhat related to generative modeling, but the scope of data and model is limited. Moreover, there are several pieces of research on utilizing feature similarity (Krupka and Tishby, 2007; Li and Li, 2008; Sandler et al., 2009; Li et al., 2017; Mollaysa et al., 2017) for discriminative problems. Despite these interests, there have been surprisingly few studies on the use of such knowledge for generative modeling in general.

5 Experiments

5.1 Datasets and Prior Knowledge

We used the following four datasets, for which we can prepare plausible prior knowledge of feature relations.

Toy data.

We created a cylinder-like toy dataset as exemplified in Figure 3. On this dataset, we know that $x_{1}$ and $x_{2}$ are more statistically dependent than $x_{1}$ and $x_{3}$ are. This knowledge is expressed in the manner of (5) as

[TABLE]

Constrained concatenated MNIST.

We created a new dataset from MNIST, namely constrained concatenated MNIST (ccMNIST, examples in Figure 3) as follows. Each image in this dataset was generated by vertically concatenating three MNIST images. Here, the concatenation is constrained so that the top and middle images have the same label, whereas the bottom image is independent of the top two. We created training and validation sets from MNIST’s training dataset and a test set from MNIST’s test dataset. On this dataset, we can anticipate that the top third and the middle third of each image are more dependent than the top and the bottom are. If each image is vectorized in a row-wise fashion as a $2352$-dim vector, this knowledge is expressed as

[TABLE]

Plant sensor data.

We used a simulated sensor dataset (Plant), which corresponds to the example in Figure 1a. The simulation is based on a real industrial chemical plant called the Tennessee Eastman process (Downs and Vogel, 1993). In the plant, there are four major unit operations: reactor, vapor-liquid separator, recycle compressor, and product stripper. On each unit operation, sensors such as level sensors and thermometers are attached. We used readings of 22 sensors. On this dataset, we can anticipate the relative dependence between sets of sensors based on the process diagram of the Tennessee Eastman process (Downs and Vogel, 1993). The knowledge set for this dataset ($\mathcal{K}_{\text{Plant}}$) contains, e.g.,

[TABLE]

where $J_{\text{compressor}}$ is the set of indices of sensors attached to the compressor, and so on. This is plausible because the reactor is directly connected to the compressor while the stripper is not. We prepared $\mathcal{K}_{\text{Plant}}$ comprising $|\mathcal{K}_{\text{Plant}}|=12$ relations based on the structure of the plant Downs and Vogel (1993).

Solar energy production data.

We used the records of solar power production222www.nrel.gov/grid/solar-power-data.html (Solar) of 137 solar power plants in Alabama in June 2006. The data of the first 20 days were used for training, the next five days were for validation, and the last five days were for test. This dataset corresponds to the example in Figure 1b because we can anticipate the dependence of features from the pairwise distances between the solar plants. We created a knowledge set $\mathcal{K}_{\text{Solar}}$ as follows; for the $i$-th plant ($i=1,\dots,137$), if the nearest plant is within 10 [km] and the distance to the second-nearest one is more than 12 [km], then add $(\{i\},\{j_{i}\},\{k_{i}\})$ to $\mathcal{K}_{\text{Solar}}$, where $j_{i}$ and $k_{i}$ are the indices of the nearest and the second-nearest plants, respectively. This resulted in $|\mathcal{K}_{\text{Solar}}|=31$.

5.2 Settings

Models.

We used factor analysis (FA, see also Zhang (2009)), VAE, and GAN. We trained FA and VAEs on all datasets and GANs on the Toy dataset. For VAEs and GANs, encoders, decoders, and discriminators were modeled using multi-layer perceptrons (MLPs) with three hidden layers. Hereafter we denote the dimensionality of VAE/GAN’s latent variable by $\dim(\bm{z})$ and the dimensionality of MLP’s hidden layer by $\dim(\text{MLP})$. We tried $\dim(\text{MLP})=2^{5},2^{6},\dots,2^{10}$ and used several values of $\dim(\bm{z})$ for the different datasets.

Baselines.

As the simplest baseline, we trained every model only with L2 regularization (i.e., weight decay). We tried another baseline, in which the weight of the last layer of decoder MLPs are regularized using $\mathcal{K}$; e.g., for $\mathcal{K}_{\text{Toy}}=\{\{1\},\{2\},\{3\}\}$, we used a regularization term

[TABLE]

where $\bm{w}_{i}$ denotes the $i$-th row of the weight matrix of the last layer of the decoder MLP, and $\epsilon$ was tuned in the same manner as $\nu_{\alpha}$. We refer to this baseline as “out-layer.” It cannot be applied to the Plant dataset because the feature sets in $\mathcal{K}_{\text{Plant}}$ do not have any one-to-one correspondences.

Hyperparameters.

We computed HSIC with $m=128$ using Gaussian kernels with the bandwidth set by the median heuristics. No improvement was observed with $m>128$. The hyperparameters were chosen based on the performance on the validation sets. The search was not intensive; $\lambda$ was chosen from three candidate values that roughly adjust orders of $L$ and $R_{\mathcal{K}}$, and $\nu_{\alpha}$ was chosen from .01 or .05.

5.3 Results

Below, the quantitative results are reported mainly with cross entropy (for ccMNIST) or reconstruction errors (for the others) as they are a universal performance criterion to examine the generalization capability. For VAEs, the significance of the improvement by the proposed method did not change even when we examined the ELBO values.

Evaluation by test set performance.

The test set performance of VAEs are shown in Tables 1 and 2. We can observe that, while the improvement by the out-layer baseline was quite marginal, the proposed regularization method resulted in significant improvement. The performance of FA slightly improved (details omitted due to space limitations).

Comparison to a dedicated model.

We compared the performance of the proposed method to that of a model designed specifically for ccMNIST dataset (termed “dedicated”). The dedicated model uses an MLP designed in accordance with the prior knowledge that the top and the middle parts are from the same digit. The performance of the dedicated model is shown in the right-most column of Table 1. We can observe that the proposed method performed intermediately between the most general case (L2 only) and the dedicated model. As designing models meticulously is time-consuming and often even infeasible, the proposed method will be useful to give a trade-off between performance and a user’s workload.

Effect of knowledge set size.

Another interest lies in how the amount of provided prior knowledge affects the performance of the regularizer. We investigated this by changing the number of tuples in $\mathcal{K}_{\text{Solar}}$ used in learning VAEs. We prepared knowledge subsets by extracting some tuples from the original set $\mathcal{K}_{\text{Solar}}$. When creating the subsets, tuples were chosen so that a larger subset contained all elements of the smaller ones. Figure 5 shows the test set performances along with the subset size. The performance is improved with a larger knowledge set. For example, the difference between $|\mathcal{K}|=6$ and $|\mathcal{K}|=31$ is significant with $p\approx.003$.

Inspection of generated samples.

We inspected the samples drawn from GANs trained on the Toy dataset. When the proposed regularizer was not applied, we often observed a phenomenon similar to “mode collapse” (see, e.g., Metz et al. (2017)) occurred as in the left-side plots of Figure 5. In contrast, the whole geometry of the Toy dataset was always captured successfully when the proposed regularizer was applied, as in the right-side plots of Figure 5.

Running time.

For example, the time needed for training a VAE for 50 epochs on the ccMNIST dataset was 150 or 117 seconds, with or without the proposed method, respectively. As the complexity is linear in $|\mathcal{K}|$, this would not be prohibitive, while any speed-up techniques will be useful.

6 Conclusion

In this work, we developed a method for regularizing generative models using prior knowledge of feature dependence, which is frequently encountered in practice and has been studied in contexts other than generative modeling. The proposed regularizer can be incorporated in off-the-shelf learning methods of many generative models. The direct extension of the current method includes the use of higher-order dependence between multiple sets of features (Pfister et al., 2017).

Acknowledgments

This work was supported by JSPS KAKENHI Grant Numbers JP19K21550, JP18H06487, and JP18H03287.

Appendix A Detailed settings of experiments

A.1 Datasets and knowledge sets

A.1.1 Toy data

In creating the Toy data, we generated 1,000 samples for a training set, 2,000 samples for a validation set, and another 2,000 samples for a test set.

A.1.2 Concatenated MNIST

In creating the ccMNIST data, we generated a training set of size 20,000 and a validation set of size 10,000 from MNIST’s training dataset and a test set of size 10,000 from MNIST’s test dataset.

A.1.3 Plant sensor data

Dataset

We used the Tennessee Eastman (TE) process data available online333github.com/camaramm/tennessee-eastman-profBraatz [retrieved 30 December 2018]., which were generated using teprob.f444The code and the related materials are also archived online at depts.washington.edu/control/LARRY/TE/download.html. originally provided by the authors of (Downs and Vogel, 1993). We downloaded d00.dat and d00_te.dat, which contained the data generated with the normal operating condition, and created three sets from them. We used the whole d00.dat for a training set (size 500) and divided d00_te.dat into a validation set (size 400) and a test set (size 560). The original data consist of 12 manipulated variables and 41 process measurement variables. Within the process measurement variables, we used the measurements of the level sensors, flow rate sensors, thermometers, and manometers, which finally composed the 22-dimensional data. As preprocessing, we normalized the data using the values of average and standard deviation of each variable of the training set. Moreover, we added noise following $\mathcal{N}(0,10^{-4})$.

Knowledge set

Based on the sensor locations shown in the process diagram of the TE process (Downs and Vogel, 1993), we classified each of the 22 variables into one of the following seven groups: Feed 1, Feed 2, Reactor, Vapor-liquid separator, Recycle compressor, Product stripper, and Purge. This grouping is built upon the unit structure of the TE process. The units are interconnected in the process as shown in Figure 6, according to which we designed the knowledge set, $\mathcal{K}_{\text{Plant}}$, as shown in the below of Figure 6.

A.1.4 Solar

Dataset

The Solar dataset was created from the data provided at the website of National Renewable Energy Laboratory555www.nrel.gov/grid/solar-power-data.html [retrieved 9 January 2019].. We downloaded the set of csv files (al-pv-2006.zip) containing power production records of solar plants in Alabama state and concatenated the contents of the 137 files as the 137-dimensional dataset. We used the data of June 2006 because the seasonal variation seemed moderate in that period. We subsampled the original data (sampled every five minutes) by 1/3 and divided them into a training set (days 1–20), a validation set (days 21–25), and a test set (days 26–30). Finally the size of each set is 1,920, 480, and 480, respectively. Because the data contain many zero values, we added $10^{-3}$ to the whole datasets and took natural logarithm, which resulted in the data valued approximately from $-6$ to $4$. As preprocessing, we normalized the data using the values of average and standard deviation of each variable of the training set. Moreover, we added noise following $\mathcal{N}(0,10^{-4})$.

Knowledge set

Knowledge set was created following scheme presented in the main manuscript.

A.2 Generative models

A.2.1 Factor analysis model

In the factor analysis (FA) model we used, the observed variable, $\bm{x}\in\mathbb{R}^{d}$, is modeled as

[TABLE]

where $\bm{z}\in\mathbb{R}^{p}$ ($p\leq d$) is a latent variable, which is inferred by

[TABLE]

The parameters, $\bm{W}$ and $\bm{\mu}$, are learned by minimizing the following square loss:

[TABLE]

This model can be regarded as a linear autoencoder with the constrained encoder and encoder. Other detailed settings are described in Section A.3.

A.2.2 Variational autoencoder

A variational autoencoder (VAE) (Kingma and Welling, 2014) is learned via maximization of the evidence lower bound (ELBO):

[TABLE]

where $p_{\theta_{\text{dec}}}(\bm{x}\mid\bm{z})$ and $q_{\phi_{\text{enc}}}(\bm{z}\mid\bm{x})$ are a decoder and an encoder having sets of parameters, $\theta_{\text{dec}}$ and $\phi_{\text{enc}}$, respectively. We modeled the encoder and the decoder using multi-layer perceptrons (MLPs) in every experiment. Other detailed settings are described in Section A.3.

A.2.3 Generative adversarial network

A generative adversarial network (GAN) (Goodfellow et al., 2014) is learned via minimization of the following cross-entropy loss:

[TABLE]

where $G_{\theta_{G}}$ and $D_{\phi_{D}}$ are a generator and a discriminator with having sets of parameters, $\theta_{G}$ and $\phi_{D}$, respectively. In GANs, the generator generates fake samples while the discriminator tries to distinguish the fake and the real data. We modeled the generator and the discriminator using MLPs. Other detailed settings are described in Section A.3.

A.3 Experimental settings

For settings that are not described below, we used the default values of PyTorch 0.4.1.

A.3.1 Hyperparameter tuning

The hyperparameters of the proposed method were chosen in accordance with the validation set performance. The set of candidates for $\nu_{\alpha}$ was just $\{.01,\ .05\}$ for every experiment. The sets of candidates for $\lambda$ are shown in Table 6. The candidate values of $\lambda$ were set so that the orders of the original loss function and the regularizer were roughly adjusted.

A.3.2 Architectures

The encoder / decoder of VAEs and the generator / discriminator of GANs were modeled using MLPs. Each MLP has fully-connected three hidden layers, each of which has the same number of units. The dimensionality of $\bm{z}$ and the number of units of MLP’s hidden layers are summarized in Table 6. The dimensionality of $\bm{z}$ was determined basically in accordance with the cumulative contributing rates (i.e., explained variance) of PCA on the Plant and Solar datasets. As for the Plant dataset, 4 principal components (PCs) explain 99% variance, 7 PCs explain 99.9%, and 11 PCs explain 99.99%; analogously for the Solar dataset.

For the MLPs, we used ReLU as the activation function except for the VAE on the Solar dataset. For the VAE on the Solar dataset, we used $\operatorname{tanh}$.

In the experiment of VAE on the ccMNIST dataset, for a reference of performance, we also tried a dedicated decoder architecture. In the dedicated decoder, the latent variable, $\bm{z}$, was divided into four parts as $\bm{z}=[\bm{z}_{1}^{\mathsf{T}}~{}\bm{z}_{2}^{\mathsf{T}}~{}\bm{z}_{3}^{\mathsf{T}}~{}\bm{z}_{4}^{\mathsf{T}}]^{\mathsf{T}}$. Then, the top third of each image is modeled only with $\bm{z}_{1}$ and $\bm{z}_{2}$, the middle third is with $\bm{z}_{1}$ and $\bm{z}_{3}$, and the bottom third is with $\bm{z}_{4}$. The detail of the dedicated decoder is shown in Figure 7b. This dedicated structure reflects the fact that the top and the middle parts of the images are from the same digit while the bottom is independent from the top two.

A.3.3 Optimization

Every optimization was done using Adam optimizer. $\alpha$, one of the parameters of Adam, was basically 0.001 with a few exception; we used $\alpha=0.0005$ for learning the FA model on the Plant dataset and for the GAN on the Toy dataset. Also, we used $\beta_{1}=0.5$ for GANs and $\beta_{1}=0.9$ otherwise.

Appendix B Detailed experimental results

B.1 Computational speed

We examined the computational time of training VAEs (using a Tesla V100 and the implementation based on PyTorch) with and without the proposed regularization. We show the runtime of training VAE on the ccMNIST dataset in Table 6 and the Solar dataset in Table 6. The settings are respectively from Sections 5.3.1 and 5.3.3. In both tables, the averages over 10 trials are shown. In summary, training with the proposed regularizer is not prohibitively slow compared to the case without the regularizer ($q=0$ in Table 6 and $|\mathcal{K}|=0$ in Table 6). Also, the runtime is not so sensitive to $q>0$ and admissible for a middle-sized knowledge set.

B.2 Evaluation by test set performance

In Tables 9, 9, and 9, we present the full version of the experimental results introduced in the main manuscript.