Deep PDF: Probabilistic Surface Optimization and Density Estimation
Dmitry Kopitkov Vadim Indelman

	
	
	Magellan Adaptive Optics Imaging of PDS 70: Measuring the Mass Accretion Rate of a Young Giant Planet within a Gapped Disk
    Kevin Wagner1,2,3,^⋆, Katherine B. Follette4, Laird M. Close1, Dániel Apai1,3,5,6, Aidan Gibbs1, Miriam Keppler6, André Müller6, Thomas Henning6, Markus Kasper7, Ya-Lin Wu8, Joseph Long1, Jared Males1, Katie Morzinski1, and Melissa McClure9
    April 14, 2020
====================================================================================================================================================================================================================================================


		
		A probability density function (pdf) encodes the entire stochastic knowledge about data distribution, where data may represent stochastic observations in robotics, transition state pairs in reinforcement learning or any other empirically acquired modality. Inferring data pdf is of prime importance, allowing to analyze various model hypotheses and perform smart decision making. However, most density estimation techniques are limited in their representation expressiveness to specific kernel type or predetermined distribution family, and have other restrictions. For example, kernel density estimation (KDE) methods require meticulous parameter search and are extremely slow at querying new points. In this paper we present a novel non-parametric density estimation approach, DeepPDF, that uses a neural network to approximate a target pdf given samples from thereof. Such a representation provides high inference accuracy for a wide range of target pdfs using a relatively simple network structure, making our method highly statistically robust. This is done via a new stochastic optimization algorithm, Probabilistic Surface Optimization (PSO), that turns to advantage the stochastic nature of sample points in order to force network output to be identical to the output of a target pdf. Once trained, query point evaluation can be efficiently done in DeepPDF by a simple network forward pass, with linear complexity in the number of query points. Moreover, the PSO algorithm is capable of inferring the frequency of data samples and may also be used in other statistical tasks such as conditional estimation and distribution transformation. We compare the derived approach with KDE methods showing its superior performance and accuracy.
		
		
		
	
	
	
	

§ INTRODUCTION

	
	
	
	Density estimation is a fundamental statistical problem which is essential for many scientific fields and application domains such as robotics, AI, economics, probabilistic planning and inference. Today's  most statistically robust techniques that work on a wide range of probabilistic distributions are non-parametric where arguably the leading is kernel density estimation (KDE). Though a large amount of research was done on these techniques, they still are computationally expensive at both estimation and query stages. They require searching for hyper-parameters like bandwidth and kernel type, and their expressiveness in representing the estimated probability density function (pdf) is high but still limited to the selected kernel. On the other hand, neural networks (NN) are known for their high flexibility and ability to approximate any function. Moreover, recently strong  frameworks (e.g. <cit.>) were developed that allow fast and sophisticated training of NNs using GPUs. 
	
	With this motivation, 
	in this paper we present a novel method to perform non-parametric density estimation using Deep Learning (DL) where we exploit the approximation power of NN in full. Unlike any existing density estimation techniques, our approach, named DeepPDF, estimates a target pdf from sampled data using our novel stochastic optimization that forces NN's output to be identical to the output of a target pdf function, thereby approximating the latter by a NN.
	
	
	The KDE method is one of the most expressive methods that is flexible enough to approximate a broad diapason of target pdfs. It represents the estimated pdf through a sum of kernel functions, with kernel for each sample point. A naive implementation of such an approach has complexity of O(Nm), where N is number of sample points and m is number of query points. Although faster KDE extensions exist (e.g. <cit.>), they are still relatively slow. Specifically, query stage of state-of-the-art fastKDE method <cit.> has quasi-linear time complexity O(m log m). Additionally, most KDE methods require to specify kernel type and bandwidth for better performance. Optimal values of these parameters are pdf-specific and parameter fine-tuning needs to be done. 
	Moreover, the target pdf may have a complex structure with different parts requiring different optimal KDE parameters. In such case, standard KDE methods will return a statistically sub-optimal solution.
	
	
	In contrast, the presented herein DeepPDF approach (see also Figure <ref>) allows to approximate the target pdf via a NN of arbitrary architecture and takes advantage of approximation robustness embedded within DL methods. Our approach trains a NN using a Probabilistic Surface Optimization (PSO). The PSO, being one of our main contributions, turns to advantage the stochastic nature of each data sample point in order to push surface that is represented by the NN. Likewise, this surface is pushed at points sampled from another known distribution. We use specially discovered stochastic pdf loss to push the neural surface in the opposite directions at data samples and at known distribution samples, where balance (minimum of pdf loss) is achieved only when the NN surface is identical to the target density. Such a balance is a substantial attribute of pdf loss which implicitly forces the estimator function (the surface) to have total integral of 1 and to be non-negative, releasing us from intricate dealing with these constraints explicitly. Moreover, training DeepPDF can be done in batch mode where at each iteration only a small number of sample points is used, lowering the requirement on memory size.
	DeepPDF is highly expressive and able to use facilities and theory developed lately for efficient and accurate DL training, including fast-converging optimizers, regularization techniques and execution on GPU cards. The pdf value at query point can be calculated by just computing output of trained NN with query point as input, which can be done very fast with complexity only depending on the NN's size, in contrast to quasi-linear complexity of fastKDE query stage.
	
	
	
	To summarize, our main contributions in this paper are as follows: 
	(a) we develop a Probabilistic Surface Optimization (PSO) and pdf loss that force any approximator function to converge to the target density;
	(b) we use PSO to approximate a target density using a NN;
	(c) we train a NN via PSO in batch mode using exponential learning rate decay for better convergence, naming the entire algorithm DeepPDF;
	and (d) we analyze different DL aspects of DeepPDF.
	Additionally, for simplicity the derivations are gathered together under Appendix, at the end of the paper.
	
	
	
	
	
	
	
	
	

§ RELATED WORK

	
	
	
	
	Statistical density estimation may be divided into two different branches - parametric and non-parametric. Parametric methods assume data to come from a probability distribution of a specific family, and infer parameters of that family, for example via minimizing the negative log-probability of data samples. Non-parametric approaches are distribution-free in the sense that they do not take any assumption over data population a priori. Instead they infer the distribution density totally from  data.
	
	
	
	The main advantage of parametric approaches is their statistical efficiency. Given the assumption of a specific distribution family is correct, parametric methods will produce more accurate density estimation for the same number of samples compared to non-parametric techniques. 
	However, as was also observed in our previous work <cit.>, in case the assumption is not entirely valid for a given population, the estimation accuracy will be poor, making parametric methods not statistically robust. For example, one of the most flexible distribution families is a Gaussian Mixture Model (GMM) <cit.>. One of its structure parameters is the number of mixtures. Using a high number of mixtures, it can represent multi-modal populations with high accuracy. Yet, in case the real unknown distribution has even higher number of modes, or sometimes even an infinite number (see Section <ref> and Figure <ref>), the performance of a GMM will be low.
	
	
	To handle the problem of unknown number of mixture components in parametric techniques, Bayesian statistics can be applied to model a prior over parameters of the chosen family. Models such as Dirichlet process mixture (DPM) and specifically Dirichlet process Gaussian mixture model (DPGMM) <cit.> can represent uncertainty about the learned distribution parameters and as such can be viewed as infinite mixture models. Although such hierarchical models are more statistically robust (flexible), they still require to manually select a base distribution for DPM, limiting its robustness. Also, Bayesian inference applied in these techniques is more theoretically intricate and computationally expensive <cit.>.
	
	
	On the other hand, non-parametric approaches can infer distributions of arbitrary form. Methods such as data histogram and kernel density estimation (KDE) <cit.> use frequencies of different points within data samples in order to conclude how a population pdf looks like. In general, these methods require more samples but also provide more robust estimation by not taking any prior assumptions. Both histogram and KDE require selection of parameters - bin width for histogram and kernel type/bandwidth for KDE which in many cases require manual parameter search <cit.>. Although an automatic parameter deduction was proposed in several studies <cit.>, it is typically computationally expensive and its performance is not always optimal. Moreover, one of the major weaknesses of the KDE technique is its time complexity during query stage. Even the most efficient KDE methods (e.g. fastKDE <cit.>) require above linear complexity (O(m log m)) in the number of query points m. In contrast, the method presented herein, DeepPDF, is a robust non-parametric approach that calculates the pdf value of a query point by a single forward pass of a NN, making query complexity to be linear in m. Such query complexity supremacy makes DeepPDF a much better alternative when many query points require evaluation. Moreover, as we will see in Section <ref>, DeepPDF provides comparable and sometimes even a much higher approximation accuracy.
	
	
	
	
	
	Using NNs for density estimation was studied for several decades <cit.>. The Mixture Density Network (MDN) presented by <cit.> represents parameters of a GMM through a NN, taking the advantage of its flexibility and expressiveness. Additionally, Williams at el. <cit.> introduced an efficient representation to parametrize a positive-definite covariance matrix of a Gaussian distribution through a NN. In image classification CNN networks produce discrete class probabilities <cit.> for each image. Yet, most of the work was done in the context of parametric approaches, arguably, due to their relative analytical simplicity.
	
	
	Another body of research in DL-based density estimation was explored
	in <cit.>. The described autoregressive methods NADE, RNADE and MADE decompose the joint distribution of a  multivariate data into a product of simple conditional densities where a specific variable ordering needs to be selected for better performance. Although these approaches provide high statistical robustness, NADE and MADE are purposed to infer density of discrete binary multi-dimensional data which limits their applicability to real world data. Further, RNADE was developed to overcome this limitation, inferring each simple conditional density through MDN. It provides a Gaussian Mixture approximation for each simple density, thus introducing assumptions about data. Instead, our approach focuses on continuous multi-dimensional data without taking any a priori data bias.
	
	
	
	A unique work combining DL and non-parametric inference was done by Baird et al. <cit.>. The authors represent a target pdf via Jacobian determinant of a bijective NN that has an implicit property of non-negativity and total integral to be 1. Additionally, their pdf learning algorithm has similarity to our pdf loss described in Section <ref>. Although the authors did not connect their approach to virtual physical forces that are pushing a neural surface, their algorithm can be seen as a simple instance of the more general DeepPDF method that we contribute.
	
	
	Additionally, in <cit.> the authors describe a loss for pdf ratio estimation which is highly similar to ours. Although our pdf loss was developed independently and through different mathematical concepts, Nam et.al. came to the same conclusion by using theory of density ratio estimation <cit.>.
	Thus,  insights from <cit.> are highly valuable and relevant also to our approach. The contributions of both our paper and of <cit.> are primarily different in the focused problem. While in <cit.> the authors estimate a ratio between two pdf densities, we use our loss to infer an entire pdf function which is a more challenging task.
	
	
	Further, in <cit.> the authors present a method, Deep Energy Estimator Network (DEEN), to infer energy function of the data pdf - an unnormalized model that is equal to (negative) log probability density function, up to a constant. Based on previous research of score matching <cit.> and its connection to denoising autoencoders <cit.>, a new approach <cit.> is presented that uses NN as a model for energy approximation. The applied loss can be shown as a force that shapes the surface (the energy approximation) to be parabolic around each data point x, eventually providing final surface that in some sense resembles unnormalized KDE estimation. It is done by optimizing the slope of a surface for each input dimension independently, enforcing the slope around each data point x to be zero at x, to be positive at nearby smaller points x^- < x and to be negative at nearby larger points x^+ > x. Such scheme will have an optimum when surface is parabolic around x. Eventually, since the parabolic shape is enforced around each data point, the final estimated energy function will have a structure very close to the real data density, thus approximating it. Note that since the applied loss deals only with the slope of energy function, the final approximation will be very far from its normalized version - that is, the normalization constant can be arbitrary big, probably depending on the initialization of NN's weights. In contrast, the presented herein method, DeepPDF, provides pdf estimation with an overall integral being very close to the true value 1 since our method forces approximation to be point-wise equal to the real pdf function. Hence, although not explicitly normalized, the pdf model produced by DeepPDF is approximately normalized.
	
	
	Finally, recently Generative Adversarial Networks (GANs) <cit.> became popular to generate new data samples (e.g. photo-realistic images). 
	GAN learns a generative model of data samples, thus implicitly learning also the data distribution. Additionally, this technique is applied for distribution transformation and for feature representation learning of high-dimensional measurements, making it extremely important for image processing problems. In contrast, our DeepPDF method learns data distribution explicitly and in theory can also be used for data generation.
	Nonetheless, there is a strong connection between the herein proposed pdf loss and recently developed novel GAN losses in <cit.>, which we are going to investigate in our future work.
	
	
	
	
	
	
	

§ PROBABILISTIC SURFACE OPTIMIZATION

	
	
	
	
	
	
	

 §.§ Derivation

	
	In this paper we present an approach to learn an arbitrary probability density function (pdf), without making any assumption on the specific distribution family. Concretely, define a multivariate random variable X_U ∈^n whose pdf UX we want to learn and samples from which are available to us. Define a neural network f of arbitrary structure (see more details below) with weights θ that gets X ∈^n as input and returns a scalar. Our goal is for f(X;θ) to approximate UX. To that end, we develop a  probabilistic optimization method, named Probabilistic Surface Optimization (PSO), that uses stochastic pdf loss as described below.
	
	
	
	
	Both f(X;θ) and UX functions can be seen as surfaces in ^n + 1:  n dimensions from X support and one additional probabilistic dimension for probability values. We look for an optimization technique with loss (w.r.t. θ) that will force surface f(X;θ) to be similar to surface UX. Moreover, we look for a loss with the property that given enough data samples from U the optimization convergence point f(X;θ^*) will coincide with UX for any X∈^n. We will derive such an optimization method in several simple steps.
	
	
	Our method is based on the main observation that performing optimization for loss L(θ, X) = f(X;θ) is identical to pushing the virtual surface f(X;θ), parametrized by θ, at point X with some point-wise force. Indeed, consider a gradient descent (GD) update of a single optimization step for data sample X, with θ_i = θ_i-1 - δ·∇ L where δ is the step size and ∇ L = ∂ f(X;θ)/∂θ.
	Then the height change (differential) of the surface  f(X';θ) at another point X' ∈^n can be approximated via a first-order Taylor expansion as (see derivation in Appendix A)
	

	
	
	
    df(X')
    	≜
    	f(X';θ_i) - f(X';θ_i-1)
    	=
    	- δ· g(X', X, θ_i-1)
    	,

	
	
	
    g(X_1, X_2, θ) ≜∂ f(X_1;θ)/∂θ^T
    	·∂ f(X_2;θ)/∂θ
    	,

	
	where the gradient similarity function g(X_1, X_2, θ) expresses gradient correlation at two points. Thus, at the point X where GD optimization was performed, the surface change will be - δ· g(X, X, θ_i-1), i.e. proportional to the L2 norm of θ gradient at this point. In other words, performing such a GD step can be thought of as pushing the surface at point X by df(X). Changing sign L(θ, X) = - f(X;θ) will simply change the direction of the differential/force, pushing the surface . Further, the above first-order Taylor expansion was observed empirically to be very accurate for learning rate δ below 0.01. During the typical training process the learning rate is less than 0.01 for most of the optimization epochs, thus the Eq. (<ref>) provides a good estimation for the real differential.
	
	
	Intuitively, when pushing a physical surface at point X, surface height at another point X' ≠ X also changes, according to elasticity and physical properties of the surface. Typically, such change diminishes with larger distance between X and X'. We argue that the same is true for the NN: according to  Eq. (<ref>), f(X';θ) at point X' ≠ X also changes due to optimization of loss at X, with differential at X' being proportional to g(X', X, ·). Therefore, g(X', X, ·) is responsible for correlation in height change at different points. The gradient similarity function g(X', X, ·) may be seen as a filter, deciding how the push at point X affects surface at point X'.
	
	In our experiments we saw that typically there is an opposite trend between g(X', X, ·) and the Euclidean distance d(X', X). Yet, g(X', X, ·) does not go entirely to zero for large d(X', X), suggesting that there exist (small) side influence even between far away points.
	
	
	
	
	
	
	
	Next, consider an example optimization algorithm where at each iteration we sample X_U ∼U and X_D ∼D, with D being some other known distribution (e.g. Gaussian). In case the loss is
	
	
    L(θ, X_U, X_D) = -f(X_U;θ) + f(X_D;θ),

	
	
	a typical GD (or any other optimization method) update at each iteration will push (along probabilistic dimension) the surface f(X;θ) at point X_U and will push the surface at point X_D. Such loss, also known as the critic loss of Wasserstein GAN <cit.>, is partially reasonable for pdf estimation since also the target pdf function UX on average returns high values for data coming from distribution U and low values for data coming from other distributions. 
	Yet, such simple loss does not converge as we will see below. 
	
	
	Further,
	
	since the considered loss in Eq. (<ref>) is probabilistic with pushed points X_U and X_D being sampled,
	we can calculate what will be (first-order Taylor approximation) the average differential [df(X)] at a  specific point X ∈^n when performing a single GD step (see Appendix B):
	
	
	
    [df(X)]
    	=
    	δ·∫
    	[UX' - DX']
    	·
    	g(X', X, θ)
    	dX'
    	.

	
	The above equation can be thought of as a  convolution of the signal [UX' - DX'] around point X according to the convolution operator g(X', X, θ).
	At each such point X' there are two opposite forces, force F_U(X') ≜UX' and force F_D≜DX', that push the surface with total force F_T(X') ≜ F_U(X') - F_D(X'). Note  these two forces express the impact for each of the two terms in  Eq. (<ref>), the term -f(X_U;θ) and the term + f(X_D;θ).
	
	
	
	
	Since  F_T(X') does not change along the optimization (it does not involve θ), it forces the surface to constantly move which prevents the optimization convergence.
	This can be simply shown in case g(X', X, θ) is a Dirac delta function as follows (for a  general function g(X', X, θ) see proof in Appendix C).
	
	
	The Dirac delta assumption may be seen as taking model f(X;θ) to its expressiveness upper bound, where for all points { X } disjoint subsets of θ are allocated to represent surface around each X. This will make f(X;θ) arbitrary flexible since each surface area can be changed independently of others. Yet, this expressiveness limit is not achievable since it will require θ of infinite size. However, the expressiveness assumption can be useful for convergence analysis since if the most flexible model can not converge then for sure the typical approximation model (e.g. neural network) will not converge also.
	
	
	
	
	
	Given this assumption, from Eq. (<ref>) we will have [df(X)] = δ· [F_U(X) - F_D(X)]
	· g(X, X, θ). It can be interpreted as the opposite forces F_U and F_D are pushing the surface f(X;θ) at each point X ∈^n, without any side influence between different points. Since both forces are point-wise constant along the optimization process, the surface f(X;θ) will be pushed asymptotically to an infinite height for points X with F_U(X) > F_D(X) ⟺UX > DX; at the rest of the points,  f(X;θ) will be pushed to a negative infinite height. Clearly, such an optimization algorithm will never converge to anything meaningful, unless properly regularized (e.g. see <cit.>). 
	
	
	Note that the derived above asymptotic values of f(X;θ)'s outputs, ±∞, typically cannot be really achieved due to limited expressiveness of a given NN architecture. Instead, usually during the optimization of loss in Eq. (<ref>), the f(X;θ) is increased to very large positive heights at areas of U density modes, and is decreased to very large negative heights at areas of D density modes. Such increase/decrease of height is unbounded due to unequal forces as was explained above. Eventually, the optimization fails due to numerical instability that involves too large/small numbers.
	
	
	One way to enforce the convergence of simple loss in Eq. (<ref>) is by adapting/changing D density along the optimization. Indeed, this is the main idea behind the contrastive divergence (CD) method presented in <cit.> and further improved in <cit.>. In CD, a point X_D is approximately sampled from the current pdf estimation f(X;θ), by performing Monte Carlo with Langevin dynamics <cit.> or Stein Variational Gradient Descent (SVGD) <cit.>. As explained above, in case at specific X ∈^n we have F_U(X) > F_D(X) ⟺UX > DX, the surface f(X;θ) will be pushed up. It in its turn will increase D at X since DX≈ f(X;θ). Eventually, such optimization will converge when UX = DX = f(X;θ). A similar idea was also applied in GAN setting in <cit.>, where sample X_D is generated by the generator NN and where generator output's density is forced to approximate current pdf estimation f(X;θ) via KL minimization between two.
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	Alternatively, consider the following pdf loss:
	
	
    L_pdf(θ, X_U, X_D) = 
    
    	-f(X_U;θ) ·DX_U
    	+
    	f(X_D;θ) ·[
    	f(X_D;θ)
    	]^mf,

	
	where DX_U is an analytically calculated pdf value at X_U according to distribution D. We use square brackets [ ·]^mf to specify  magnitude function terms that only produce magnitude coefficients but whose gradient w.r.t. θ is not calculated (stop gradient in Tensorflow <cit.>). Such terms do not change the direction of the gradient, affecting only its magnitude.
	
	The loss in Eq. (<ref>) is still pushing the surface f(X;θ) up at points X_U (via term -f(X_U;θ) ·DX_U) and down at points X_D (via term f(X_D;θ) ·[
	f(X_D;θ)
	]^mf). However, the expected differential [df(X)] has changed due to the new loss terms (see Appendix D):
	
	
    [df(X)]
    	=
    	δ·∫
    	[F_U(X') - F_D(X'; θ)]
    	·
    	g(X', X, θ)
    	dX'
    	,

	
	with forces F_U and F_D being changed to:
	
	
	
	
	
	
	
	
	
    F_U(X') =
    	UX'·DX',

	
	
	
    F_D(X'; θ) = 
    	DX'· f(X';θ).

	
	
	The above differential will become zero only at F_U(X') = F_D(X'; θ), and this is where GD will converge, using an  appropriate learning rate decay. Note, that unlike loss in Eq. (<ref>), the force F_D is now a function of θ and can adapt along optimization to stabilize the opposite force F_U.
	
	
	The PSO optimization is based on pushing the probabilistic surface, represented by a NN f(X;θ), through forces F_U and F_D (see also Figure <ref>). It does so by applying the pdf loss defined in Eq. (<ref>) iteratively, where this loss can also be written in batch mode, with multiple samples from data densities (see Section <ref>). Additionally, pdf loss can be also derived in an alternative way as a sampled approximation of (see Appendix E):
	

	
	
    L(θ) =
    	_X ∼D[
    	[
    	UX - f(X;θ)
    	]^2
    	]
    	,

	
	similarly to the ratio estimation methods in <cit.>. Although such derivation is  faster, it lacks the intuition about the physical forces and the virtual surface. Yet, this intuition is very helpful for the below analysis of learning stability.
	
	
	
	
	
	

 §.§ PSO Analysis

	
	
	
	
	
	
	
	
	
	Here we analyze PSO convergence by assuming, for simplicity, that g(X', X, θ) is a Dirac delta function. In other words, we assume that there is no correlation and side-influence among pushed points. As explained above, this can also be seen as an expressiveness limit assumption for extremely flexible surface f(X;θ) where a push at X affects only the surface height at X. We emphasize this assumption is undertaken only to simplify the converegence analysis, while in our experiments g(X', X, θ) was not constrained in any way. Nevertheless, the provided results (Section <ref>) support the below analysis. A more general analysis is part of future research.
	
	
	Considering Eqs. (<ref>)-(<ref>),
	
	we can analyze the  work done by PSO in two different perspectives. First is the typical optimization analysis of loss in Eq. (<ref>), while second is the analysis of dynamical physical forces on surface in some physical system using Eqs. (<ref>) and (<ref>). Both perspectives are two different sides of the same coin. Below we will use each of these perspectives to analyze the properties of PSO. 
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	Considering an arbitrary point X ∈^n and using a physical perspective, the balance between up and down forces at the point happens when:
	F_U(X) = F_D(X; θ) 
	  ⟹  
	f(X;θ) = UX.
	
	
	
	
	
	
	
	
	Convergence to such balance can be explained as follows. In case f(X;θ) < UX, we can see from Eqs. (<ref>) and (<ref>) that F_U(X) > F_D(X; θ), so the surface at X will be pushed up by the pdf loss. If, on the other hand, f(X;θ) > UX will be true, then it applies F_U(X) < F_D(X; θ) and the opposite will happen, i.e. the surface will be pushed down. Eventually, when balance f(X;θ) = UX is achieved, both up and down forces will be equal and cancel each other; the surface at this point will not be affected anymore. Thus, PSO will indeed optimize f(X;θ) to be UX in a time limit.
	
	
	
	
	
	
	
	
	Intuitively, considering the above analysis and Eq. (<ref>) we would like the gradient similarity g(X', X, θ) to resemble kernel function with bandwidth being close to average distance between training points. In this way at each point X ∈^n the surface f(X;θ) will be modified only by pushes at nearby training points with whom X has similar statistical frequency (pdf value). Yet, in typical NN the gradient similarity doesn't have the kernel shape, although it has rough opposite correlation with distance between two points. This in its turn raises an interesting research direction for NN architecture with more kernel-like function g(X', X, θ). We leave this direction for the future research.
	
	
	
	
	

 §.§ Choosing D

	
	
	For D we can use any candidate distribution that covers the support of a target distribution U. This can be concluded from the following reasons. Consider points {X : UX > 0  and DX = 0 }. From a physical perspective, Eqs. (<ref>) and (<ref>) tell us that both F_U and F_D are zeros at such points. From an optimization perspective, we can conclude the same by noting that the second term in Eq. (<ref>), f(X_D;θ) ·[
	f(X_D;θ)
	]^mf, will never be activated at such points since X cannot be sampled from distribution D. Moreover, the first term, -f(X_U;θ) ·DX_U, has a zero multiplier inside, thus nullifying it. Therefore, at such points there will be no optimization done and f(X;θ) typically will be just interpolation of surface heights at points around X.
	
	Further, consider points {X : UX = 0  and DX > 0 }. The force F_U will not be activated since X cannot be sampled from distribution U. The force F_D will push surface down at such X but using optimization perspective we can see that once surface height becomes negative, f(X;θ) < 0, the force F_D will change its direction and push surface up due to sign canceling in term f(X_D;θ) ·[
	f(X_D;θ)
	]^mf in Eq. (<ref>). Thus, for such points X the f(X;θ) will oscillate around zero and given an appropriate learning rate decay it will eventually converge to zero which is also the value of U at this point.
	
	Concluding from above, the only requirement on distribution D is to wrap support of U. Using such D, the surface f(X;θ) at points {X : UX > 0 } will converge to UX, and at points {X : UX = 0 } it will converge to zero. Empirically we observed that the above "wrap" requirement on D is typically enough for proper learning convergence. However, some other constraints on selecting D may also exist. For instance, if there is an area where {X : UX≫DX} we may have many samples from U but no samples from D. Then the force balance in this area will not be reached and the surface will be pushed to infinity, causing instability in the entire optimization process. Intuitively, for proper convergence the ratio UX/DX should be bounded in some dynamical range. A more detailed investigation is required to clarify this point which we consider as future work.
	
	
	
	
	
	
	
	
	Typically limits of U's support can be easily derived from the training dataset and a proper candidate for D can be picked up. Furthermore, in case we do not know U's support range and still want to learn a  target pdf without experiencing optimization divergence, we can use a "support-safe" version of pdf loss: 
	

	
	
    L_pdf^sprt-safe(θ, X_U, X_D) = 
    
    	-f(X_U;θ) ·DX_U·[
    	sign(P_max - f(X_U;θ))
    	]^mf
    	+
    	f(X_D;θ) ·[
    	f(X_D;θ)
    	]^mf,

	
	where the new sign term constrains the maximal height of surface f(X;θ) to be P_max. Note that such a  constraint will prevent f(X;θ) to accurately approximate target U at points {X : UX > P_max}.
	
	
	
	Finally, the density estimation can also be done via a variant of the pdf loss without the magnitude function terms:
	
	
    L_pdf(θ, X_U, X_D) = 
    	-f(X_U;θ) ·DX_U
    	+ 
    	[
    	f(X_D;θ)
    	]^2,

	
	which has exactly the same optimization gradient as the standard pdf loss from Eq. (<ref>). 
	
	
	
	
	
	
	
	
	

§ DEEP PDF APPROACH

	
	
	
	
	As explained above, the described PSO and pdf loss L_pdf(θ, X_U, X_D) will force a NN f(X;θ) to approximate UX.
	
	We can also use PSO in batch regime by sampling sets of points from distributions U and D, B^U ≜{X_U^i}_i = 1^N and B^D ≜{X_D^i}_i = 1^N, and updating θ appropriately by minimizing loss L_T(θ) = 1/N∑_i = 1^N L_pdf(θ, X_U^i, X_D^i). Using batch of points for a single θ update can be seen as many points on the surface f(X;θ) being pushed /in order to reach a force balance. Such an update should be much more accurate since the number of points in batch will reduce the stochastic nature of the pdf loss, and make real differential df(X) closer to its expected value in Eq. (<ref>). The entire batch learning algorithm, named DeepPDF, is detailed more scrupulously in Appendix F.
	
	
	Below, we provide details on different deep learning aspects of the presented approach.
	
	
	
	
	
	
	
	
	

  
Learning rate and convergence
	
	
	
	
	
	
	In our experiments we saw that convergence to the target pdf happens very fast, requiring us only to apply "good" learning rate decay policy. Empirically observed, during the first several thousand iterations the surface f(X;θ) reaches high similarity with surface UX, and its heights at various points continue to oscillate near the ground truth heights of target surface. Such oscillation around the target surface can be explained by the fact that at each iteration only N points are pushed up/down. These N points are sampled at each iteration from U/D and /forces are applied only at them during this iteration. Thus, considering only one iteration the true forces are very stochastic and may be far from their averaged expressions in Eqs. (<ref>) and (<ref>), but still are not too far since  f(X;θ) becomes similar to the target. 
	Intuitively, we can see that the proper thing to do is to wait till f(X;θ) starts to oscillate around some surface and then slowly decrease the learning rate δ. Such a heuristic policy can be achieved for example by exponential decay learning rate which at first stage is constant letting f(X;θ) to get closer to UX, sequentially reducing learning rate to its minimal value δ_min via 
	δ = a ·
	b^floor(t / s) + δ_min,
	
	
	
	
	
	
	
	where a is the start learning rate, b is decay rate, s is number of steps before each learning rate reduction and t is iteration index. We used Adam optimizer <cit.> with start learning rate a = 10^-3 and decay rate b = 0.5. In our experiments we observed that using s = 200000 and δ_min = 10^-7 provides better convergence and higher density estimation accuracy. Allowing PSO to take more time to converge empirically was seen to provide more accurate approximation. Still, it is up to the network architect to define learning rate decay policy according to time vs accuracy preference.
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	Convergence of pdf loss is stochastic as was explained above. Empirically observed, the pdf loss from Eq. (<ref>) does not provide clear convergence metric since it is too noisy and just oscillates around zero for most part of the optimization. Yet, the value of pdf loss without stop-gradient (Eq. (<ref>)) provides better insight about the current optimization convergence, since it has direct connection to the squared error in Eq. (<ref>). It reduces monotonically during the optimization and may be used to
	
	decide when accuracy stops improving (e.g. for early stopping).
	Alternatively, it is possible to use learning rate with exponential decay and stop training when learning rate get very close to  δ_min
	
	and θ does not change anymore.
	
	Further, value of Eq. (<ref>) is also useful for choosing the most accurate model between several trained candidates.
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	

  
Overfitting
	
	
	Like any other machine learning method, DeepPDF can overfit to training dataset. This can happen for example, if the used neural network f(X;θ) is highly expressive and the amount of available points from U is small. In such case f(X;θ) will converge to zero-height surface with peaks at available training points. In a sense, this would not be a mistake since the distribution of such a small and scarce training dataset has indeed flat pdf with peaks around training points. However, in some cases data is very difficult to generate and we need to do the best with a little data that we have. 
	
	In order to detect  overfitting, the typical approach also applicable here is to have a separate validation dataset of data generated from U, and check its performance (for example, via pdf loss in Eq. (<ref>)). Once validation performance starts to decrease,  overfitting takes place.
	
	
	To deal with the overfitting problem, the network architect can decrease expressiveness of the network by reducing its size or adding L2 regularization on weights θ. Empirically we saw that such methods can indeed help. Reducing expressiveness of the network will of course increase the approximation error between f(X;θ) and UX. Yet, we cannot expect high approximation accuracy when there is only scarce data amount available.
	
	
	Still, there is one additional solution to the problem of a little data. We can perform data augmentation to add some diversity noise into our samples, based on our knowledge about the real distribution U. This technique is very popular in image classification field where it proved to be very useful <cit.>, producing object classifier that is robust to various camera light conditions and which is less prone to the overfitting. In case of the density estimation such method will increase variability of sample points and reduce the probability pike around each such point. Yet, the specific data augmentation process will also introduce bias (prior knowledge) in our pdf inference - the properties of the output pdf surface (e.g. its smoothness) will be affected by it. If such prior knowledge is correct and if data is indeed scarce, the data augmentation can tremendously improve accuracy of the density estimation.
	
	
	
	
	
	
	
	
	

  
Positive constraint
	
	
	
	
	
	
	
	
	
	Since the network f(X;θ) is an approximation of a non-negative function UX, it makes sense to explicitly constrain its output to be positive/non-negative. For example, we can use activation functions Relu or Softplus at the end of the network. Alternatively, we can use an exponential function whose outputs are positive. This will be especially appealing since the input of an exponential function will be the log-probability function logUX. Such log-probability is very useful, for example, in estimation applications.
	
	However, using such activation functions (especially exponential) can damage the numerical stability of  learning. For example, exp (a) will return infinity for a > 88. During training the network signal is not stable and can go above 88 depending on learning rate, causing training process to fail. Such a problem can be dealt with by initializing θ weights in a different way or by reducing learning rate. Yet, such handling is time-consuming and by itself may slow down the learning process.
	
	
	Alternatively, we can use f(X;θ) with an  unconstrained output. As was explained in Section <ref>, the pdf loss will eventually force it to be non-negative at points sampled from both U and D. Of course a little oscillation around zero may occur. Also, at other points, where {X : UX = 0  and DX = 0 }, the f(X;θ) may return negative values since these areas are not affected by PSO optimization (at least not directly). But this can be easily handled by a proxy pdf function:
	
	
    f̅(X;θ) = 
    	
    	0,   if  f(X;θ) < 0  or DX = 0
    
    	f(X;θ),                 otherwise

	
	
	
	
	
	
	

  
Inner network architecture
	
	
	
	Expressiveness of a NN f(X;θ), i.e. how well it can approximate a specific surface UX, depends on its inner structure. Number, type and size of inner layers, all these will eventually affect f(X;θ)'s elasticity and the final approximation accuracy. Yet, while the optimal inner architecture of the network depends on the target pdf UX and its surface properties, in Section <ref> we show that even a relatively simple network can approximate various target pdfs with high accuracy, being competitive to the state-of-the-art KDE methods. Still, finding an inner architecture that works perfectly for all target surfaces is an open research question.
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	Additionally, we tested the usage of batch-normalization (BN) and of dropout within our FC layers. We saw empirically that BN prevents DeepPDF from learning the  target pdf. It can be explained by the fact that BN changes the distribution of each batch while DeepPDF basically relies on this distribution to approximate the target pdf. Similarly, dropout was seen to damage accuracy of density estimation. Thus, BN and dropout harm performance of our method and should be avoided for proper DeepPDF convergence.
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	

  
Batch size vs. estimation quality
	
	
	
	As was explained above the batch size N can be interpreted as amount of points on surface where we push /during a single iteration. The bigger is N and the more points are pushed at once, the closer are the real applied up/down forces to their expected values in Eqs. (<ref>) and (<ref>). Thus, bigger N will produce more stable, fast and accurate approximation results, which was also observed empirically. Up to memory resource constraints, we suggest to use large batches for best performance.
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	

§ EXPERIMENTAL EVALUATION

	
	
	
	
	We evaluate our DeepPDF approach on three target densities, whose analytical expressions can be found in Appendix H. The inferred pdfs conceptually represent complicated multi-modal and discontinuous densities, which need to be inferred in real world applications (e.g. estimating measurement model in robotics domain).
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	In all scenarios, network f(X;θ) contains 3 FC layers of size 1024 with Relu non-linearity, followed by an output FC layer of size 1 with no activation function. Also, we use uniform distribution for DX. During learning the batch size is 1000 and overall size of training samples X_U ∼UX for all evaluated methods is 10^8.
	The effect of smaller dataset size on the inference accuracy is part of future investigation.
	The learning rate was decayed through the same decay policy for all scenarios, without use of early stopping.
	The overall number of iterations was around 6 · 10^6 for each scenario.
	
	Once trained, we use pdf proxy function from
	Eq. (<ref>) as our pdf function approximation.
	DeepPDF was trained with Adam optimizer <cit.> using GeForce GTX 1080 Ti GPU card.
	
	
	
	
	For comparison, we infer the same densities with KDE method, using two KDE implementations: First is KernelDensity provided by sklearn library <cit.> and represents a standard KDE approach. Second is fastKDE <cit.> which is a  state-of-the-art KDE extension that provides better accuracy and computational efficiency. Optimal parameters for Sklearn-KDE's were found via grid search, while fastKDE performs automatic parameter search.
	
	Additionally, we experiment with a simplified version of PSO in Appendix G
	
	
	to demonstrate its power to infer probability of a specific single point.
	
	
	
	
	
	
	
	
	
	
	

 §.§ Learning 2D distributions

	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	Here we use DeepPDF to estimate density of two 2D distributions shown in Figures <ref> (Cosine) and <ref> (Columns). The corresponding density approximations are shown in Figures <ref> and <ref>. As can be seen, both distributions were inferred with high accuracy.
	
	
	
	In Table <ref> we compare performance of all the three approaches: DeepPDF, Sklearn-KDE and fastKDE. The used performance metric is average L2 distance 1/N∑_i = 1^N [ UX_i - f(X_i) ]^2, where f(·) is an  approximation function and points { X_i }_i=1^N form a 257 × 257 grid in UX's support. As can be seen, all methods provide comparable accuracy. For Columns pdf the accuracy of DeepPDF is 36 times higher than other techniques, probably because there are density areas with different optimal bandwidths. KDE methods are using the same bandwidth for all inner kernels and cannot approximate such target pdf optimally. However, there is no such problem in DeepPDF since it is highly flexible and capable of approximating a wide variety of functions without the need to specify bandwidth. For Cosine pdf the accuracy of fastKDE is only 1.9 times lower than DeepPDF. The reason for such relatively comparable accuracy is the smoothness of the Cosine pdf that implies single optimal bandwidth for entire pdf surface.
	
	
	
	
	Also note that query stage of both KDE techniques is very slow - it takes about two days for Sklearn-KDE and about an hour for fastKDE to evaluate 66000 points. In contrast, once trained DeepPDF can evaluate such query number in matter of seconds.
	
	In context of KDE methods, the sample points can be viewed as parameters that represent target pdf approximation function, with such representation being inefficient for the query stage. Thus, the above KDE approaches can be seen as being parametrized by parameter vector of size identical to the size of training dataset (10^8 in our experiments) multiplied by data dimension. On other hand, the size of representation parameter vector θ in our DeepPDF method is only about 2.1 · 10^6 (2103297 for 2D distributions). Moreover, post-training investigation of ∂ f(X;θ)/∂θ on 10^3 test points showed that only small subset of weights has non-zero gradient on at least one test point. That is, only several tens of thousands of parameters are used to represent surface f(X;θ), suggesting that most of θ can be pruned away. Yet, with more compact representation DeepPDF successfully infers the target pdf with even higher accuracy relatively to KDE methods. This supports one of the main arguments of this paper that neural network is more flexible to approximate probability densities than sum of kernels.
	
	
	
	
	
	
	
	

 §.§ Learning 3D range measurement model

	
	
	
	
	
	Further, we evaluate the above density techniques to infer a 3D distribution. We use a RangeMsr distribution that is similar to the  measurement model of a range sensor which is highly adapted in robotics domain. Each sample point has a structure [ x, y, f ] where x and y are uniformly sampled from range [-2, 2] and symbolically represent robot location. Additionally, f is sampled from 𝒩(√(x^2 + y^2), 1) and symbolize distance measurement between robot and some oracle located at origin, plus some noise. Using specific f = 1, the RangeMsr distribution is visualized in 3D space in Figure <ref>. As can be seen, such distribution has infinitely many modes and approximating it with simple methods like GMM would fail.
	
	In Figure <ref> we can see that DeepPDF's approximation of RangeMsr is indeed very accurate. In Table <ref> we compare L2 performance of DeepPDF vs KDE methods. As can be seen, our method performs much more accurate estimation than the standard KDE implementation by sklearn library <cit.>, and even is superior to a state-of-the-art fastKDE approach <cit.>.
	
	
	
	
	
	
	
	
	
	

§ CONCLUSIONS AND FUTURE WORK

	
	
	
	
	
	
	In this paper we presented a novel non-parametric density estimation technique: DeepPDF. It uses a neural network representation of a target pdf which produces high expressiveness, allowing it to approximate an arbitrary multi-dimensional pdf given sample points of thereof. Albeit it requires to design the neural network structure for pdf appropximator, the same simple composition of fully-connected layers with Relu non-linearity was able to approximate all experimented densities with high accuracy. Moreover, once trained, DeepPDF can be efficiently queried for a pdf value at a specific point by a simple network forward pass, making it a much better alternative to the compared KDE methods. We enforce the non-negativity of DeepPDF pdf approximation via proxy function from Eq. (<ref>). The total integral is not explicitly constrained to be 1, yet we verified it via numerical integration and saw it to be almost identical to the correct value 1.
	
	
	Further, we believe the described probabilistic surface optimization (PSO) is a conceptually novel approach to infer data density and can be applied in other statistical frameworks besides density estimation, such as new data generation. As future research directions, we will investigate both DeepPDF theoretical and empirical convergence properties, and its dependence on data dimension and batch size. Additionally, we will explore the relation between the accuracy of DeepPDF density estimation and size of the training dataset.
	
	









10
	
	Antoniak74aos
	C. E. Antoniak.
	Mixtures of Dirichlet processes with applications to Bayesian
	nonparametric problems.
	Annals of Statistics, 2:1152–1174, 1974.
	
	Arjovsky17arxiv
	Martin Arjovsky, Soumith Chintala, and Léon Bottou.
	Wasserstein gan.
	arXiv preprint arXiv:1701.07875, 2017.
	
	Baird05ijcnn
	Leemon Baird, David Smalenberger, and Shawn Ingkiriwang.
	One-step neural network inversion with pdf learning and emulation.
	In 2005 IEEE International Joint Conference on Neural Networks,
		IJCNN'05, volume 2, pages 966–971. IEEE, 2005.
	
	Bengio00nips
	Yoshua Bengio and Samy Bengio.
	Modeling high-dimensional discrete data with multi-layer neural
	networks.
	In Advances in Neural Information Processing Systems (NIPS),
	pages 400–406, 2000.
	
	Bishop94tr
	C.M. Bishop.
	Mixture density networks.
	Technical report, Aston University, Birmingham, 1994.
	
	Caffe_url
	Caffe.
	<caffe.berkeleyvision.org>.
	
	Duong05sjs
	Tarn Duong and Martin L Hazelton.
	Cross-validation bandwidth matrices for multivariate kernel density
	estimation.
	Scandinavian Journal of Statistics, 32(3):485–506, 2005.
	
	Germain15icml
	Mathieu Germain, Karol Gregor, Iain Murray, and Hugo Larochelle.
	Made: Masked autoencoder for distribution estimation.
	In Intl. Conf. on Machine Learning (ICML), pages 881–889,
	2015.
	
	Goodfellow14nips
	Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley,
	Sherjil Ozair, Aaron Courville, and Yoshua Bengio.
	Generative adversarial nets.
	In Advances in Neural Information Processing Systems (NIPS),
	pages 2672–2680, 2014.
	
	Gorur10jcst
	Dilan Görür and Carl Edward Rasmussen.
	Dirichlet process gaussian mixture models: Choice of the base
	distribution.
	Journal of Computer Science and Technology, 25(4):653–664,
	2010.
	
	Gulrajani17nips
	Ishaan Gulrajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron C
	Courville.
	Improved training of wasserstein gans.
	In Advances in Neural Information Processing Systems (NIPS),
	pages 5769–5779, 2017.
	
	Heidenreich13asta
	Nils-Bastian Heidenreich, Anja Schindler, and Stefan Sperlich.
	Bandwidth selection for kernel density estimation: a review of fully
	automatic selectors.
	AStA Advances in Statistical Analysis, 97(4):403–433, 2013.
	
	Hinton06
	G.E. Hinton, S. Osindero, and Y. Teh.
	A fast learning algorithm for deep belief nets.
	Neural Computation, 18:1527–1554, 2006.
	
	Hinton02nc
	Geoffrey E Hinton.
	Training products of experts by minimizing contrastive divergence.
	Neural Computation, 14(8):1771–1800, 2002.
	
	Hyvarinen05mlr
	Aapo Hyvärinen.
	Estimation of non-normalized statistical models by score matching.
	Journal of Machine Learning Research, 6(Apr):695–709, 2005.
	
	Hyvarinen08nc
	Aapo Hyvärinen.
	Optimal approximation of signal priors.
	Neural Computation, 20(12):3087–3110, 2008.
	
	Kim16arxiv
	Taesup Kim and Yoshua Bengio.
	Deep directed generative models with energy-based probability
	estimation.
	arXiv preprint arXiv:1606.03439, 2016.
	
	Kingma14arxiv
	Diederik P Kingma and Jimmy Ba.
	Adam: A method for stochastic optimization.
	arXiv preprint arXiv:1412.6980, 2014.
	
	Kopitkov18iros
	D. Kopitkov and V. Indelman.
	Robot localization through information recovered from cnn
	classificators.
	In IEEE/RSJ Intl. Conf. on Intelligent Robots and Systems
		(IROS). IEEE, October 2018.
	Accepted.
	
	Krizhevsky12nips
	Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.
	Imagenet classification with deep convolutional neural networks.
	In Advances in neural information processing systems, pages
	1097–1105, 2012.
	
	Larochelle11aistats
	Hugo Larochelle and Iain Murray.
	The neural autoregressive distribution estimator.
	In Proceedings of the Fourteenth International Conference on
		Artificial Intelligence and Statistics, pages 29–37, 2011.
	
	Ledig16arxiv
	Christian Ledig, Lucas Theis, Ferenc Huszár, Jose Caballero, Andrew
	Cunningham, Alejandro Acosta, Andrew Aitken, Alykhan Tejani, Johannes Totz,
	Zehan Wang, et al.
	Photo-realistic single image super-resolution using a generative
	adversarial network.
	arXiv preprint, 2016.
	
	Liu16nips
	Qiang Liu and Dilin Wang.
	Stein variational gradient descent: A general purpose bayesian
	inference algorithm.
	In NIPS, pages 2378–2386, 2016.
	
	Liu17arxiv
	Qiang Liu and Dilin Wang.
	Learning deep energy models: Contrastive divergence vs. amortized
	mle.
	arXiv preprint arXiv:1707.00797, 2017.
	
	MacEachern98
	S. N. MacEachern and P. Muller.
	Estimating mixture of dirichlet process models.
	Journal of Computational and Graphical Statistics, 7:223–238,
	1998.
	
	Mao17iccv
	Xudong Mao, Qing Li, Haoran Xie, Raymond YK Lau, Zhen Wang, and Stephen Paul
	Smolley.
	Least squares generative adversarial networks.
	In Intl. Conf. on Computer Vision (ICCV), pages 2813–2821.
	IEEE, 2017.
	
	McLachlan88
	G.J. McLachlan and K.E. Basford.
	Mixture Models: Inference and Applications to Clustering.
	Marcel Dekker, New York, 1988.
	
	Mohamed16arxiv
	Shakir Mohamed and Balaji Lakshminarayanan.
	Learning in implicit generative models.
	arXiv preprint arXiv:1610.03483, 2016.
	
	Mroueh17nips
	Youssef Mroueh and Tom Sercu.
	Fisher gan.
	In Advances in Neural Information Processing Systems (NIPS),
	pages 2510–2520, 2017.
	
	Nam15toias
	Hyunha Nam and Masashi Sugiyama.
	Direct density ratio estimation with convolutional neural networks
	with application in outlier detection.
	IEICE TRANSACTIONS on Information and Systems,
	98(5):1073–1079, 2015.
	
	OBrien16csda
	Travis A O’Brien, Karthik Kashinath, Nicholas R Cavanaugh, William D Collins,
	and John P O’Brien.
	A fast and objective multidimensional kernel density estimation
	method: fastkde.
	Computational Statistics & Data Analysis, 101:148–160, 2016.
	
	Pytorch_url
	PyTorch.
	<pytorch.org>.
	
	Radford15arxiv
	Alec Radford, Luke Metz, and Soumith Chintala.
	Unsupervised representation learning with deep convolutional
	generative adversarial networks.
	arXiv preprint arXiv:1511.06434, 2015.
	
	Raphan11nc
	Martin Raphan and Eero P Simoncelli.
	Least squares estimation without priors or supervision.
	Neural Computation, 23(2):374–420, 2011.
	
	Saremi18arxiv
	Saeed Saremi, Arash Mehrjou, Bernhard Schölkopf, and Aapo Hyvärinen.
	Deep energy estimator networks.
	arXiv preprint arXiv:1805.08306, 2018.
	
	Scott15book
	David W Scott.
	Multivariate density estimation: theory, practice, and
		visualization.
	John Wiley & Sons, 2015.
	
	Sethuraman82sdtrt
	Jayaram Sethuraman and Ram C Tiwari.
	Convergence of dirichlet measures and the interpretation of their
	parameter.
	In Statistical decision theory and related topics III, pages
	305–315. Elsevier, 1982.
	
	Silverman18
	Bernard W Silverman.
	Density estimation for statistics and data analysis.
	Routledge, 2018.
	
	Sklearn_KDE_url
	Sklearn KernelDensity.
	<scikit-learn.org/stable/modules/generated/sklearn.neighbors.KernelDensity.html>.
	
	Smolensky86rep
	Paul Smolensky.
	Information processing in dynamical systems: Foundations of harmony
	theory.
	In D. E. Rumelhart and J. L. McClelland, editors, Parallel
		Distributed Processing, volume 1. The MIT press, Cambridge, MA, 1986.
	
	Sugiyama12book
	Masashi Sugiyama, Taiji Suzuki, and Takafumi Kanamori.
	Density ratio estimation in machine learning.
	Cambridge University Press, 2012.
	
	Tensorflow_url
	TensorFlow.
	<www.tensorflow.org>.
	
	Uria2013nips
	Benigno Uria, Iain Murray, and Hugo Larochelle.
	Rnade: The real-valued neural autoregressive density-estimator.
	In Advances in Neural Information Processing Systems (NIPS),
	pages 2175–2183, 2013.
	
	Williams96nc
	Peter M Williams.
	Using neural networks to model conditional multivariate densities.
	Neural Computation, 8(4):843–854, 1996.
	
	Zhao16arxiv
	Junbo Zhao, Michael Mathieu, and Yann LeCun.
	Energy-based generative adversarial network.
	arXiv preprint arXiv:1609.03126, 2016.
	








§ APPENDIX A: SURFACE CHANGE AFTER GRADIENT DESCENT UPDATE OF L(Θ, X) = F(X;Θ)


Performing a single gradient descent step w.r.t. loss L(θ, X) = f(X;θ) at a specific optimized point X results in a weights update θ_i = θ_i-1 - δ·∇ L, where δ is step size and ∇ L = ∂ f(x;θ)/∂θ|_x = X, θ = θ_i-1. After such update, the surface height at point X, f(X;θ), can be approximated via a Taylor expansion as:


    f(X;θ_i) =
    f(X;θ_i-1)
    +
    (θ_i - θ_i-1)^T
    ·∂ f(x;θ)/∂θ|_x = X, θ = θ_i-1
    +
    
    +
    ·
    (θ_i - θ_i-1)^T
    · H_θ·
    (θ_i - θ_i-1)
    +
    ⋯
    =
    
    =
    f(X;θ_i-1)
    -
    δ·
    (∂ f(x;θ)/∂θ|_x = X, θ = θ_i-1)^T
    ·∂ f(x;θ)/∂θ|_x = X, θ = θ_i-1
    +
    
    +
    δ^2 ·
    (∂ f(x;θ)/∂θ|_x = X, θ = θ_i-1)^T
    · H_θ·
    (∂ f(x;θ)/∂θ|_x = X, θ = θ_i-1)
    +
    ⋯
    ,


where:


    H_θ≜∂^2 f(x;θ)/∂θ^2|_x = X, θ = θ_i-1
    .


Above we can see the first two elements of Taylor expansion where the second element depends on a quadratic step size δ^2. The (learning) step size is typically less than one. Further, assuming it is very small number (δ≪ 1), f(X;θ_i) can be approximated good enough by a first-order Taylor expansion as:


    f(X;θ_i) = f(X;θ_i-1)
    -
    δ·
    (∂ f(x;θ)/∂θ|_x = X, θ = θ_i-1)^T
    ·∂ f(x;θ)/∂θ|_x = X, θ = θ_i-1
    ,


and define


    df(X)
    ≜
    f(X;θ_i) - f(X;θ_i-1)
    =
    - δ· g(X, X, θ_i-1)
    ,




    g(X_1, X_2, θ) ≜∂ f(X_1;θ)/∂θ^T
    ·∂ f(X_2;θ)/∂θ
    ,


where g(·) function calculates the cross-gradient similarity between two different points X_1 and X_2, i.e. the dot product between θ gradients at these two points. Further, g(X, X, θ_i-1) is the L2 norm of θ gradient at point X. 

Similarly, the change of surface height at point X', which is different from the optimized point X (i.e. X' ≠ X), can be approximated as:


    df(X')
    =
    - δ· g(X', X, θ_i-1)
    .


Note that the cross-gradient similarity function is symmetric, i.e. g(X', X, ·) = g(X, X', ·). Thus, from the above we can see that for any two arbitrary points X ∈^n and X' ∈^n, pushing at one point will change height at another point according to the cross-gradient similarity g(X', X, ·). Therefore, g(X', X, ·) is responsible for correlation in height change at different points. Given |g(X', X, ·)| is small (or zero), optimizing loss L (pushing surface) at point X will not affect the surface at point X'. The opposite is also true: if |g(X', X, ·)| is large, then pushing the  surface at point X will have a big impact at point X'.

▪









§ APPENDIX B: SURFACE CHANGE AFTER GRADIENT DESCENT UPDATE OF SIMPLE LOSS



In case of simple loss L(θ, X_U, X_D) = -f(X_U;θ) + f(X_D;θ) where 
X_U and X_D are sampled from U and D respectively, consider weights update θ_i = θ_i-1 - δ·∇ L. Similarly to Eq. (<ref>), the change of surface height at some point X ∈^n can be approximated as a first-order Taylor expansion:


    df(X)
    =
    δ· g(X_U, X, θ_i-1)
    - δ· g(X_D, X, θ_i-1)
    .


Further, considering the stochastic nature of X_U and X_D, the average differential [df(X)] can be calculated as:


    [df(X)]
    =
    δ·∫UX'·
    g(X', X, θ_i-1)
    dX'
    - 
    
    -
    δ·∫DX'·
    g(X', X, θ_i-1)
    dX'
    =
    δ·∫
    [UX' - DX']
    ·
    g(X', X, θ_i-1)
    dX'
    .


Again, we see that g(X', X, ·) has the role of a filter that controls correlations between points. In contrast, the term [UX' - DX'] depends only on point X'.



▪





§ APPENDIX C: PROOF OF DIVERGENCE OF SIMPLE LOSS


In case of simple loss L(θ, X_U, X_D) = -f(X_U;θ) + f(X_D;θ) the average differential [df(X)],


    [df(X)]
    =
    δ·∫
    [UX' - DX']
    ·
    g(X', X, θ_i-1)
    dX'
    ,


contains signal term [UX' - DX'] and filter term g(X', X, θ_i-1). Since the signal term does not depend on θ and is unchanging along the optimization, there is no stable point for the optimization convergence. Thus, the learning process diverges. This can be proved by contradiction as follows.


Assume, in contradiction, that the above simple loss converges at optimization step t to the loss's local minimum. Then, for the following optimization steps,  θ, and so g(X', X, θ), will remain (almost) constant - GD will stay in the minimum area. Hence, the integral in Eq. (<ref>) becomes a constant for each given point X, since θ does not change anymore. This constant is non-zero, unless both and densities are identical UX' = DX'; or the learned f(X;θ) is flat (g(X', X, θ) = 0 for any two points). Yet, such specific cases are degenerate and we exclude them.

Concluding from the above, after optimization step t the [df(X)] became non-zero constant for each X. Therefore, on average the surface height continues to change by a non-zero constant at each point X, going towards ±∞. Yet, such change of surface f(X;θ) is impossible since θ is (almost) constant after t. We got a contradiction, meaning that the above convergence assumption is wrong; that is, the  simple loss  L(θ, X_U, X_D) = -f(X_U;θ) + f(X_D;θ) cannot converge due to non-zero signal [UX' - DX'].





§ APPENDIX D: SURFACE CHANGE AFTER GRADIENT DESCENT UPDATE OF PDF LOSS



In case of pdf loss L_pdf(θ, X_U, X_D) = 

-f(X_U;θ) ·DX_U
+
f(X_D;θ) ·[
f(X_D;θ)
]^mf where 
X_U and X_D are sampled from U and D respectively, consider weights update θ_i = θ_i-1 - δ·∇ L. Similarly to Eq. (<ref>), the change of surface height at some point X ∈^n can be approximated as first-order Taylor expansion:


    df(X)
    =
    δ·DX_U· g(X_U, X, θ_i-1)
    - δ·
    f(X_D;θ)
    · g(X_D, X, θ_i-1)
    .


Further, considering stochastic nature of X_U and X_D, the average differential [df(X)] can be calculated as:


    [df(X)]
    =
    δ·∫UX'·DX'·
    g(X', X, θ_i-1)
    dX'
    - 
    
    -
    δ·∫DX'·
    f(X';θ)
    ·
    g(X', X, θ_i-1)
    dX'
    =
    
    =
    δ·∫DX'·
    [UX' - f(X';θ)]
    ·
    g(X', X, θ_i-1)
    dX'
    .



▪






§ APPENDIX E: PDF LOSS AS EMPIRICAL APPROXIMATION OF CONTINUOUS SQUARED LOSS



Consider a squared loss between NN surface f(X;θ) and a target function UX, weighted by the pdf D of distribution:


    L(θ) =
    _X ∼D[
    [
    UX - f(X;θ)
    ]^2
    ]
    =
    ∫DX·[
    UX - f(X;θ)
    ]^2
    dX
    .


Similarly to the ratio estimation methods in <cit.>, we can show that the above loss can be minimized by minimizing our pdf loss, as follows:


    L(θ) =
    ∫DX·[
    UX]^2
    +
    DX·[
    f(X;θ)
    ]^2
    -
    2 ·DX·UX·
    f(X;θ)
    dX
    .


Since the first term doesn't depend on θ, we can rewrite the loss as:


    L(θ) =
    ∫DX·[
    f(X;θ)
    ]^2
    -
    2 ·DX·UX·
    f(X;θ)
    dX
    =
    
    =
    ∫DX·[
    f(X;θ)
    ]^2
    dX
    -
    2 ·∫UX·DX·
    f(X;θ)
    dX
    =
    
    =
    _X ∼D [f(X;θ)^2]
    -
    2 ·_X ∼U 
    [DX·
    f(X;θ)]
    .


The above can be approximated by N samples { X_U^i }_i = 1^N from distribution U and N samples { X_D^i }_i = 1^N from distribution D as:


    L(θ) =
    1/N∑_i = 1^N
    f(X_D^i;θ)^2
    -
    2 ·1/N∑_i = 1^NDX_U^i·
    f(X_U^i;θ)


where for notation simplicity we can rewrite it as one-sample version:


    L(θ) =
    f(X_D;θ)^2
    -
    2 ·DX_U·
    f(X_U;θ)
    .


Next, note that the above loss has the same gradient w.r.t. θ as the following loss:


    L(θ) =
    2
    ·
    f(X_D;θ)
    ·[
    f(X_D;θ)
    ]^mf
    -
    2 ·DX_U·
    f(X_U;θ)
    ,


where [ ·]^mf specifies  magnitude function terms that only produce magnitude coefficients but whose gradient w.r.t. θ is not calculated (stop gradient in Tensorflow <cit.>).
Thus, both losses are equal in context of the gradient-based optimization.

After changing order of terms in Eq. (<ref>), and deviding by 2 (e.g. moving it into learning rate), we will finally get the pdf loss:


    L_pdf(θ, X_U, X_D) = 
    
    -f(X_U;θ) ·DX_U
    +
    f(X_D;θ) ·[
    f(X_D;θ)
    ]^mf
    .



▪





















































































§ APPENDIX F: DEEPPDF ALGORITHM IN MORE DETAIL











§ APPENDIX G: POINT OPTIMIZATION EXPERIMENT








In this experiment we apply our loss only at one specific point X_z = [0, 0]^T in ^2 space, to demonstrate that the described PSO works in a simple case of single point optimization. The used loss is:


    L_zero(θ) = 
    -f(X_z;θ) ·_U[X = X_z]
    + 
    f(X_z;θ) ·[
    f(X_z;θ)
    ]^mf·_D[X = X_z],


where _U[X = X_z] is indicator random variable that is equal to 1 with probability p_U = 0.9 and is 0 with probability 1 - p_U. Similarly, _D[X = X_z] is 1 with probability p_D = 0.15 and is 0 otherwise. In such a way we simulate a sampling process from  UX and DX where we assume UX_z = p_U and DX_z = p_D. According to the analysis from Section

[3]
, the balance between up/down forces using L_zero(θ) should be at:


    F_U(X_z) = F_D(X_z) 
      ⟹  
    p_U = p_D · f(X_z;θ)
      ⟹  
    f(X_z;θ) = 
    p_U/p_D
    =
    6.


We apply Adam optimizer <cit.> with loss L_zero(θ) during several thousand iterations and show value of f(X_z;θ) along this optimization in Figure <ref>. As can be seen, f(X_z;θ) quickly reaches height 6, following stochastic oscillation around this value. With learning rate decay the oscillations become smaller and the final f(X_z;θ) value ends to be very close to 6. This simple experiment demonstrates the power of PSO to approximate height of a target pdf. Additionally, it provides insight on how PSO can be used to infer density of discrete random variable instead of continuous, which was the main focus of this paper.







§ APPENDIX H: USED DENSITIES - EQUATIONS


Below are the analytical density functions for distributions used in our experiments.


	
	
	
  * 
		Cosine:
	
	
	
	
    Ux_1, x_2
    	=
    	1/μ·[
    	cos
    	(4
    	x_1
    	·
    	x_2
    	)
    	+ 1
    	]
    	,   μ = 17.631302268269998

	
	
	
  * 
		Columns:
	
	
	
	
    Ux_1, x_2
    	=
    	p(x_1)
    	·
    	p(x_2)
    	,

	
	where p(·) is mixture distribution with 5 components:

	Uniform(-2.3, -1.7), 𝒩(-1.0, std = 0.2), 𝒩(0.0, std = 0.2), 𝒩(1.0, std = 0.2), Uniform(1.7, 2.3)

	Each component has weight 0.2.
	
	
	
	
	
  * 
		RangeMsr:
	
	
	
	
    Ux,y,f
    	=
    	p(x) · p(y) · p(f|x,y)
    	,

	
	where
	
	
    p(x)
    	=
    	p(y)
    	=
    	Uniform(-2.0, 2.0)
    	,

	
	and
	
	
    p(f|x,y) = 𝒩(√(x^2 + y^2), std = 1.0)
    	.

	
	
	





