title: Low-rank Characteristic Tensor Density Estimation Part I: Foundations
authors: Magda Amiridi Nikos Kargas Nicholas D. Sidiropoulos
text: 
  
  
  
  
  
  Low-rank Characteristic Tensor Density Estimation Part I: Foundations
      Magda Amiridi,
          Nikos Kargas,
          and Nicholas D. Sidiropoulos, Fellow, IEEE
  
  
  M. Amiridi and N.D. Sidiropoulos are with the Department of ECE, University of Virginia, Charlottesville, VA 22904. Author e-mails: (ma7bx,nikos)@virginia.edu. N. Kargas was with the Department of ECE, University of Minnesota; he is now with Amazon, Cambridge, U.K. Author e-mail: karga005@umn.edu
      December 30, 2023
  =======================================================================================================================================================================================================================================================================================================================================================================================================
  
  
  
  
  Effective non-parametric density estimation is a key challenge in high-dimensional multivariate data analysis. In this paper, we propose a novel approach that builds upon tensor factorization tools. Any multivariate density can be represented by its characteristic function, via the Fourier transform. If the sought density is compactly supported, then its characteristic function can be approximated, within controllable error, by a finite tensor of leading Fourier coefficients, whose size depends on the smoothness of the underlying density. This tensor can be naturally estimated from observed and possibly incomplete realizations of the random vector of interest, via sample averaging. In order to circumvent the curse of dimensionality, we introduce a low-rank model of this characteristic tensor, which significantly improves the density estimate especially for high-dimensional data and/or in the sample-starved regime. By virtue of uniqueness of low-rank tensor decomposition, under certain conditions, our method enables learning the true data-generating distribution. We demonstrate the very promising performance of the proposed method using several toy, measured, and image datasets. 
  
  
  
  Statistical learning, Probability Density Function (PDF), Characteristic Function (CF), Tensors,   Rank,  Canonical  Polyadic  Decomposition  (CPD).
  
  
  
  
  
  
  § INTRODUCTION
  
  
  Density estimation is a fundamental yet challenging problem in data analysis and machine learning. Density estimation is the task of learning the joint Probability Density Function (PDF) from a set of observed data points, sampled from an unknown underlying data-generating distribution. A model of the density function of a continuous random vector provides a complete description of the joint statistical properties of the data and can be used to perform tasks such as computing the most likely value of a subset of elements (“features”) conditioned on others, computing any marginal or conditional distribution, and deriving optimal estimators, such as the minimum mean squared error (conditional mean) estimator. Density estimation has a wide range of applications including classification <cit.>, clustering <cit.>, data synthesis <cit.>, data completion <cit.> and reconstruction related applications <cit.>, as well as learning statistical regularities such as skewness, tail behavior, multi-modality or other structures present in the data <cit.>. 
  
  Existing work on density estimation can be mainly categorized into parametric approaches such as Gaussian Mixture Models (GMM) <cit.>, and non-parametric approaches such as Histogram <cit.> and Kernel Density Estimation (KDE) <cit.>. A density model must be expressive – flexible enough to represent a wide class of distributions, and tractable and scalable (computationally and memory-wise) at the same time (expressivity-tractability trade-off). Over the last several years, explicit feed-forward neural network based density estimation methods <cit.> have gained increasing attention as they provide a tractable way to evaluate high-dimensional densities point-wise. On the other hand implicit generative models such as generative adversarial networks <cit.> and variational autoencoders <cit.> can be used to obtain models which allow effective and efficient sampling. 
  
  In this paper, we develop a novel non-parametric method for multivariate PDF estimation based on tensor rank decomposition – known as Canonical Polyadic Decomposition (CPD) <cit.>. CPD is a powerful model that can parsimoniously represent high-order data tensors exactly or approximately, and its distinguishing feature is that under certain reasonable conditions it is unique – see <cit.> for a recent tutorial overview. We show that any compactly supported continuous density can be approximated, within controllable error, by a finite characteristic tensor of leading complex Fourier coefficients, whose size depends on the smoothness of the density. This  characteristic tensor can be naturally estimated via sample averaging from realizations of the random vector of interest. 
  
  The main challenge, however, lies in the fact that the size of this tensor (the number of model parameters in the Fourier domain) grows exponentially with the number of random variables – the length of the random vector of interest. In order to circumvent this “curse of dimensionality” (CoD) and further denoise the naive sample averaging estimates, we introduce a low-rank model of the characteristic tensor, whose degrees of freedom (for fixed rank) grow linearly with the random vector dimension. Low-rank modeling significantly improves the density estimate especially for high-dimensional data and/or in the sample-starved regime. By virtue of uniqueness of low-rank tensor decomposition, under certain conditions, our method enables learning the true data-generating distribution. 
  
  
  In order to handle incomplete data (vector realizations with missing entries) as well as scaling up to high-dimensional vectors, we further introduce coupled low-rank decomposition of lower-order characteristic tensors corresponding to smaller subsets of variables that share `anchor' variables, and show that this still enables recovery of the global density, under certain conditions. As an added benefit, our approach yields a generative model of the sought density, from which it is very easy to sample from. This is because our low-rank model of the characteristic tensor admits a latent variable naive Bayes interpretation. A corresponding result for finite-alphabet random vectors was first pointed out in <cit.>. In contrast to <cit.>, our approach applies to continuous random vectors possessing a compactly supported multivariate density function. From an algorithmic standpoint, we formulate a constrained coupled tensor factorization problem and develop a Block Coordinate Descent (BCD) algorithm. 
  
  The main results and  contributions of this paper can be summarized as follows:
  
  
    * We show that any smooth compactly supported multivariate PDF can be approximated by a finite tensor model, without using any prior or data-driven discretization process. We also show that truncating the sampled multivariate characteristic function of a random vector is equivalent to using a finite separable mixture model for the underlying distribution. Under these relatively mild assumptions, the proposed model can approximate any high dimensional PDF with approximation guarantees.By virtue of  uniqueness of CPD, assuming low-rank in the Fourier domain, the underlying multivariate density is identifiable. 
  
  
    * We show that high dimensional joint PDF recovery is possible under low tensor-rank conditions, even if we only observe subsets (triples) of variables. This is a key point that enables one to handle incomplete realizations of the random vector of interest. To the best of our knowledge, no other generic density estimation approach allows this. To tackle this more challenging version of the problem, we propose an optimization framework based on coupled tensor factorization. Our approach jointly learns lower-order (3-dimensional) characteristic functions, and then assembles tensor factors to synthesize the full characteristic function model.
  
  
    * The proposed model allows efficient and low-complexity inference, sampling, and density evaluation. In that sense, it a more comprehensive solution that neural density evaluation or neural generative models. We provide convincing experimental results on toy, image, and measured datasets that corroborate the effectiveness of the proposed method.
  
  
  This is the first of a two-part paper. The second part builds on this foundation to develop a joint compression (nonlinear dimensionality reduction) and compressed density estimation framework that offers additional flexibility and scalability, but does not provide a density estimate in the original space as the “baseline” method in this first part does. Each approach has its own advantages, but the second builds upon the first. It is therefore natural to present them as Part I and Part II.   
  
  
  
  § BACKGROUND
  
  
  
   §.§ Related work
  
  
  Density estimation has been the subject of extensive research in statistics and the machine learning community. Methods for density estimation can broadly be classified as either parametric or non-parametric. Parametric density estimation assumes that the data are drawn from a known parametric family of distributions, parametrized by a fixed number of tunable parameters. Parameter estimation is usually performed by maximizing the likelihood of the observed data. One of the most widely used parametric models is the Gaussian Mixture Model (GMM). GMMs can approximate any density function if the number of components is large enough <cit.>. 
  However, a very large number of components may be required for good approximation of the unknown density, especially in high dimensions. Increasing the number of components introduces computational challenges and may require a large amount of data <cit.>. Misspecification and inconsistent estimation is less likely to occur with nonparametric density estimation. 
  
  Nonparametric density estimation is more unassuming and in that sense “universal”, but the flip-side is that it does not scale beyond a small number of variables (dimensions). The most widely-used approach for nonparametric density estimation is Kernel Density Estimation (KDE) <cit.>. The key idea of KDE is to estimate the density by means of a sum of kernel functions centered at the given observations. However, worst-case theoretical results show that its performance worsens exponentially with the dimension of the data vector <cit.>. 
  
  Our approach falls under nonparametric methods, and is motivated by Orthogonal Series Density Estimation (OSDE) <cit.>, a powerful non-parametric estimation methodology. OSDE approximates a probability density function using a truncated sum of orthonormal basis functions, which may be trigonometric, polynomial, wavelet etc. However, OSDE becomes computationally intractable even for as few as 10 dimensions, since the number of parameters grows exponentially with the number of dimensions. Unlike OSDE, our approach is able to scale to much higher dimensions.
  
  Recently, several density evaluation and modeling methods that rely on neural networks have been proposed. The Real-valued Neural Autoregressive Distribution Estimator (RNADE) <cit.> is among the best-performing neural density estimators and has shown great potential in scaling to high-dimensional distribution settings. These so-called autoregressive models (not to be confused with classical AR models for time-series) decompose the joint density as a product of one-dimensional conditionals of increasing conditioning order, and model each conditional density with a parametric model. Normalizing Flows (NF) <cit.> models start with a base density e.g., standard Gaussian, and stack a series of invertible transformations with tractable Jacobian to approximate the target density. Masked Autoregressive Flow (MAF) <cit.> is a type of NF model, where the transformation layer is built as an autoregressive neural network. These methods do not construct a joint PDF model but rather serve for point-wise density evaluation. That is, for any given input vector (realization), they output an estimate of the density evaluated at that particular input vector (point). For small vector dimensions, e.g., two or three, it is possible to evaluate all inputs on a dense grid, thereby obtaining a histogram-like density estimate; but the curse of dimensionality kicks in for high vector dimensions, where this is no longer an option. Additionally, these methods cannot impute more than very few missing elements in the input, for the same reason (grid search becomes combinatorial).
  
  
  
   §.§ Notation
  
  In this section we briefly present notation conventions and some tensor algebra preliminaries. We use the symbols 𝐱, 𝐗, 𝐗 for vectors, matrices and tensors respectively. We use the notation 𝐱(n), 𝐗(:,n), 𝐗(:,:,n) to refer to a particular element of a vector, a column of a matrix and a slab of a tensor. Symbols 𝐱_2, 𝐗_F, and 𝐱_∞ correspond to L_2 norm, Frobenius norm, and infinity norm.
  Symbols ∘, ⊛, ⊙  denote the outer, Hadamard and Khatri-Rao product respectively. The vectorization operator is denoted as vec(𝐗), vec(𝐗) for a matrix and tensor respectively. Additionally, diag(𝐱)∈ℂ^K × K denotes the diagonal matrix with the elements of vector 𝐱∈ℂ^K on its diagonal. The set of integers {1,…, K} is denoted as [K]. 
  
  
  
   §.§ Relevant tensor algebra
  
  
   An N-way tensor Φ∈ℂ^K_1 × K_2 ×⋯× K_N is a multidimensional array whose elements are indexed by N indices. Any tensor can be decomposed as a sum of F rank-1 tensors
   
      Φ = ∑_h=1^Fλ(h) 𝐀_1(:,h) ∘𝐀_2(:,h) ∘⋯∘𝐀_N(:,h),
  
  where 𝐀_n ∈ℂ^K_n× F and constraining the columns 𝐀_n(:,h) to have unit norm, the real scalar λ(h) absorbs the h-th rank-one tensor’s scaling. A visualization is shown in Figure <ref> for the case of N=3. 
  
  
  We use Φ = [[λ, 𝐀_1,…,𝐀_N]] to denote the decomposition. When F is minimal, it is called the rank of tensor Φ, and the decomposition is called Canonical Polyadic Decomposition (CPD). A particular element of the tensor is given by Φ(k_1,k_2,…,k_N) = ∑_h=1^F λ(h) ∏_n=1^N 𝐀_n(k_n,h). The vectorized form of the tensor can be expressed as vec(Φ)= (⊙_n=1^N 𝐀_n)λ. We can express the mode-n matrix unfolding which is a concatenation of all mode-n `fibers' of the tensor as  Φ^(n) = (⊙_k ≠ n𝐀_k)  diag(λ)𝐀_n^T,
  where (⊙_k ≠ n𝐀_k)= 𝐀_N ⊙⋯⊙𝐀_n+1⊙𝐀_n-1⊙⋯⊙𝐀_1. 
  
  A key property of the CPD  is that the rank-1 components are unique under mild conditions. For learning probabilistic latent variable models and latent representations, the uniqueness of tensor decomposition can be interpreted as model identifiability. A model is identifiable, if and only iff there is a unique set of parameters that is consistent with what we have observed. 
  
   <cit.>: 
  Let k_𝐀 be the Kruskal rank of 𝐀, defined as the largest integer k such that every k columns of 𝐀 are linearly independent. Given  Φ=[[λ,𝐀_1,…,𝐀_N]], if ∑_n=1^N k_𝐀_n≥2F+N-1, then the rank of Φ is F and the decomposition of Φ in rank-one terms is unique. 
  
  Better results allowing for higher tensor rank are available for generic tensors of given rank.
   <cit.>: Given Φ =[[λ,𝐀_1,𝐀_2,𝐀_3]], assume, without loss of generality, that I_1≤ I_2≤ I_3. Let α,β be the largest integers such that 2^α≤ I_1 and 2^β≤ I_2. If F≤2^α+β-2 the decomposition of Φ in rank-one terms is unique almost surely.
  
  
  
  
  § A CHARACTERISTIC FUNCTION APPROACH
  
  The characteristic function of a random variable X is the Fourier transform of its PDF, and it conveys all information about X. The characteristic function can be interpreted as an expectation: the Fourier transform at frequency ν∈ℝ is E[e^j ν X]. Similarly, the multivariate characteristic function is the multidimensional Fourier transform of the density of a random vector X, which can again be interpreted as the expectation E[e^j ν^T  X], where ν is a vector of frequency variables. The expectation interpretation is crucial, because ensemble averages can be estimated via sample averages; and whereas direct nonparametric density estimation at point x requires samples around x, estimating the characteristic function enables reusing all samples globally, thus enabling better sample averaging and generalization. This point is the first key to our approach. The difficulty, however, is that pinning down the characteristic function seemingly requires estimating an uncountable set of parameters. We need to reduce this to a finite parameterization with controllable error, and ultimately distill a parsimonious model that can learn from limited data and still generalize well. In order to construct an accurate joint distribution estimate that is scalable to high dimensions without making explicit and restrictive prior assumptions (such as a GMM model) on the nature of the density, and without requiring huge amounts of data, we encode the following key ingredients into our model.
  
  
    * Compactness of support. In most cases, the random variables of interest are bounded, and these bounds are known or can be relatively easily estimated.  
  
    * Continuity of the underlying density and its derivatives. The joint distribution is assumed to be sufficiently smooth in some sense, which enables the use of explicit or implicit interpolation.
  
    * Low-rank tensor modeling.
  We show that joint characteristic functions can be represented as higher order tensors. In practice these tensor data are not unstructured. Low-rank tensor modeling provides a concise  representation that captures the salient characteristics (the principal components) of the data distribution in the Fourier domain. 
  
  
  
  
   §.§ The Univariate Case
  
  
  
  Before we delve into the multivariate setting, it is instructive to examine the univariate case. Given a real-valued random variable X with compact support S_X, the Probability Density Function (PDF) f_X and its corresponding Characteristic Function (CF) Φ_X form a Fourier transform pair:
  
      Φ_X(ν)    :=∫_S_Xf_X(x)e^jν xdx  = E[e^jν X], 
       
      f_X(x)   :=1/2π∫_-∞^∞Φ_X(ν)e^-jν xdν.
  
  Note that Φ_X(0)=∫_-∞^∞f_X(x)d x=1. Without loss of generality, we can apply range normalization and mean shifting so that sX+c ∈ [0,1] – the transformation is invertible. We may therefore assume that S_X = [0,1]. Every PDF supported in [0,1] can be uniquely represented over its support by an infinite Fourier series,
  
      f_X(x) = ∑_k=-∞^∞Φ_X[k]e^-j2π k x,
  
  where Φ_X[k]=Φ_X(ν)|_ν=2π k,    k∈ℤ. This shows that countable parameterization through samples of the characteristic function suffices for compactly supported densities. But this is still not enough - we need a finite parametrization.  Thankfully, if f_X is sufficiently differentiable in the sense that f_X∈ C^p i.e.,  all its derivatives ∂ f_X/∂ x, ∂^2 f_X/∂ x^2, ⋯, ∂^p f_X/∂ x^p exist and are continuous we have that
  
   (e.g., see <cit.>): If f_X∈ C^p, then 
  
      |Φ_X[k]|= 𝒪(1/1+|k|^p).
   
  
  
  It is therefore possible to use a truncated series 
  
      f_X(x) = ∑_k=-K^KΦ_X[k]e^-j2π k x,
   
  with proper choice of K that will not incur significant error. Invoking Parseval's Theorem 
  
      f_X - f_X _2^2=∑_|k|>K|Φ_X[k]|^2,
   
  which is controllable by the smoothing parameter K. The k-th Fourier coefficient
  
      Φ_X[k] = ∫_0^1f_X(x)e^j2π k xd x =E[e^j2π k X]
   
  can be conveniently estimated via the sample mean
  
      Φ_X[k]=1/M∑_m=1^Me^j2π k x_m
   
  Here M is the number of available realizations of the random variable X.
  
  A toy example to illustrate the idea is shown in Figure <ref>. For this example, we are given M=500 realizations of a random variable X, which is a mixture of two Gaussian distributions with means μ_1=0.35,  μ_2=0.7 and standard deviations σ_1=0.1, σ_2=0.08. The recovered PDF is very close to the true PDF using only 9 coefficients of the CF.
  
  
   §.§ The Multivariate Case
  
  
  In the multivariate case, we are interested in obtaining an estimate f_ X of the true density f_ X of a random vector X:= [X_1,…, X_N]^T. The joint or multivariate characteristic function of X is a function Φ_ X:ℝ^N→ℂ defined as 
  
      Φ_ X(ν):=
      E[e^j ν^T  X],
  
  where ν := [ν_1,…,ν_N]^T.
  For any given ν, given a set of realizations {𝐱_m }_m=1^M, we can estimate the empirical characteristic function of the sequence as 
  
      Φ_ X(ν) =
      
      1/M∑_m=1^M e^j ν^T𝐱_m.
  
  
  Under mixing conditions such that sample averages converge to ensemble averages, the corresponding PDF can be uniquely recovered via the multidimensional inverse Fourier transform
  
      f_ X(x)=1/(2π)^N∫_ℝ^NΦ_ X(ν)e^-jν^T  x dν.
  
  If the support of the joint PDF f_ X(𝐱) is contained within the hypercube S_ X=[0,1]^N, then similar to the univariate case, it can be represented by a multivariate Fourier series
  
      f_ X(x)= ∑_k_1=-∞^∞⋯∑_k_N=-∞^∞Φ_ X[k]e^-j2π𝐤^T𝐱,
      where Φ_ X[𝐤] = Φ_ X(ν)|_ν=2π𝐤, 𝐤 = [k_1,…,k_N]^T.
  
   (see e.g., <cit.>): For any p ∈ℕ, if the partial derivatives ∂^θ_1/∂ x_1^θ_1⋯∂^θ_N/∂ x_N^θ_N f_ X(𝐱) exist and are absolutely integrable for all θ_1,…,θ_N with ∑_n=1^N θ_n ≤ p then the rate of decay of the magnitude of the 𝐤-th Fourier coefficient |Φ_ X[𝐤]| obeys |Φ_ X[𝐤]|= 𝒪(1/1+𝐤_2^p).
  
  The smoother the underlying PDF, the faster its Fourier coefficients and the approximation error tend to zero. Thus we can view the joint PDF through the lens of functions with only low frequency harmonics. Specifically,
  it is known <cit.>, <cit.> that the approximation error of the truncated series with absolute cutoffs {K_n}_n=1^N is upper bounded by
  
      f_ X- f_ X_∞≤ C  ∑_n=1^N ω_n (∂^θ_n/∂ x_n^θ_n f_ X,1/1+K_n)/(1+K_n)^θ_n,
  
  where ω_n(f_ X,δ):=
  
      sup|x_j-x_j'|≤δ|f_ X(x_1,…,x_j,…,x_N)-f_ X(x_1,…,x_j',…,x_N)|,
  
  and 
  
      C = C_2 (1+ C_1  ∏_n=1^N log K_n).
   
  C_1,C_2 are constants independent of K_n. The  smoother  the  underlying  PDF,  the  smaller  the  obtained  finite parametrization error. It follows that we can approximate f_ X
  
      f_ X(𝐱)=∑_k=-K_1^K_1⋯∑_k_N=-K_N^K_NΦ_ X[𝐤]e^-j2π𝐤^T𝐱.
  
  The truncated Fourier coefficients can be naturally represented using an N-way tensor Φ where
  
      Φ(k_1, …, k_N) = Φ_X[𝐤].
  
  
  
  
  § PROPOSED APPROACH: BREAKING THE CURSE OF DIMENSIONALITY
  
  
  We have obtained a finite parameterization with controllable and bounded error, but the number of parameters (2K_1 + 1) ×⋯× (2K_N+1)  obtained by truncating Φ_ X as above grows exponentially with N. This curse of dimensionality can be circumvented by focusing on the principal components of the resulting tensor, i.e., introducing a low-rank parametrization of the Characteristic Tensor obtained by truncating the multidimensional Fourier series. Keeping the first F principal components, the number of parameter reduces from order of K_1 ×⋯× K_N to order of (K_1 + ⋯ + K_N)F. Introducing the rank-F CPD in Equation (<ref>), one obtains the approximate model
  
      f̃_ X(𝐱)=
      ∑_k_1=-K^K⋯∑_k_N=-K^K∑_h=1^F p_H(h)∏_n=1^NΦ_X_n|H=h[k_n]
      e^-j2π k_n x_n,
  
  where H can be interpreted as a latent (H for `hidden') random variable, Φ_X_n|H=h[k_n] is the characteristic function of X_n conditioned on H=h
  
      Φ_X_n|H=h[k_n]   :=Φ_X_n|H=h(ν|h)|_ν=2π k_n
         =𝔼_X_n|H=h[e^j2π k_n X_n],
  
  and we stress that for high-enough F, this representation is without loss of generality – see, e.g., <cit.>. For the rest of the paper, we consider K=K_1=⋯=K_N for brevity. 
  
  By linearity and separability of the multidimensional Fourier transformation it follows that 
  
      f̃_ X(𝐱)    =∑_h=1^F p_H(h)∏_n=1^N∑_k_n=-K^KΦ_X_n|H=h[k_n]e^-j2π k_n x_n
         =∑_h=1^F p_H(h)∏_n=1^N f_X_n|H(x_n|h).
  
  This generative model can be interpreted as mixture of product distributions <cit.>. The joint PDF f_ X is a mixture of F separable component PDFs, i.e., there exists a `hidden' random variable H taking values in {1,…,F} that selects the operational component of the mixture, and given H the random variables X_1,…,X_N become independent (See Figure <ref> for visualization of this model). We have thus shown the following result: 
  
   Truncating the multidimensional Fourier series (sampled multivariate characteristic function) of any compactly supported random vector is equivalent to approximating the corresponding multivariate density by a finite mixture of separable densities.  
   
   Thus, by choosing appropriate K and F, it is possible to represent and approximate any compactly supported density that it is sufficiently smooth by the proposed model. See Figures <ref>, <ref> where we showcase how each parameter affects the modeling of complex structures in 2D synthetic datasets. 
  
  
  
  Conversely, if one assumes that the sought multivariate density is a finite mixture of separable densities, then it is easy to show that the corresponding characteristic function is likewise a mixture of separable characteristic functions:
  
  
  Φ_X(ν)    = E[e^j ν^T X]
  
     =E_H[E_X|H[e^j ν_1 X_1 ⋯e^j ν_N X_N]]
  
     =E_H[ Φ_X_1|H(ν_1|H) ⋯Φ_X_n|H(ν_N|H)]
  
     =∑_h=1^F p_H(h) ∏_n=1^N Φ_X_n|H(ν_n|h).
  
  
  If we sample the above on any finite N-dimensional grid, we obtain an N-way tensor and its polyadic decomposition. Such decomposition is unique, under mild conditions <cit.>. It follows that: 
   A compactly supported multivariate (N ≥ 3) mixture of separable densities is identifiable from (samples of) its characteristic function, under mild conditions.
  
  
  
  The above analysis motivates the following course of action. Given a set of realizations {𝐱_m  }_m=1^M,
  
  
    * estimate
  
      Φ[𝐤] = 1/M∑_m=1^M e^j 2 π𝐤^T𝐱_m,
  
  
    * fit a low-rank model
  
      Φ[𝐤] ≈∑_h=1^F p_H(h) ∏_n=1^NΦ_X_n|H=h[k_n],
  
  
    * and invert using
  
      f_ X(𝐱)= ∑_h=1^F p_H(h)∏_n=1^N f_X_n|H(x_n|h),  where 
       
       f_X_n|H(x_n|h)=∑_k_n=-K^KΦ_X_n|H=h[k_n]e^-j2π k_n x_n.
  
  
  
  When building any statistical model, identifiability is a fundamental question. 
  A statistical model is said to be identifiable when, given a sufficient number of observed data, it is possible to uniquely recover the data-generating distribution. When applying a non-identifiable model, different structures or interpretations may arise from distinct parametrizations that explain the data equally well. Most deep generative models do not address the question of identifiability, and thus may fail to deliver the true latent representations that generate the observations. Our approach is fundamentally different, because it builds on rigorous and controllable Fourier approximation and identifiability of the characteristic tensor.
  
  In the Appendix (Section <ref>), we provide additional statistical insights regarding the proposed methodology, including the asymptotic behavior of the empirical characteristic function and the mean squared error reduction afforded by low-rank tensor modeling in the characteristic function domain. 
  
  Two issues remain. First, uniqueness of CPD only implies that each rank-one factor is unique, but leaves scaling/counter-scaling freedom in p_H and the conditional characteristic functions. To resolve this, we can use the fact that each conditional characteristic function must be equal to 1 at the origin (zero frequency). Likewise, p_H must be a valid probability mass function. These constraints fix the scaling indeterminacy. 
  
  We note here that, under certain rank conditions (see Section <ref>) on the Fourier series coefficient tensor, the proposed method ensures that the reconstructed density is positive and integrates to one, as it should. This is due to the uniqueness properties of the Fourier series representation and the CPD: if there exists a density that generates a low-rank characteristic tensor, and that tensor can be uniquely decomposed, the sum of Fourier inverses of its components is unique, and therefore equal to the generating density. Under ideal low-rank conditions, this is true even if we ignore the constraints implied by positivity when we decompose the characteristic tensor in the Fourier domain. This is convenient because strictly enforcing those in the Fourier domain would entail cumbersome spectral factorization-type (positive semidefinite) constraints.
  
  We therefore propose the following formulation:
  
      min     Φ - [[λ,𝐀_1,…,𝐀_N  ]] _F^2 
      subject to     λ≥0, 1^Tλ = 1, 
           𝐀_n(K+1,:)=1^T,   n=1… N,
  
  where 𝐀_n(K+1+k_n,h) holds Φ_X_n|H=h[k_n], and λ(h) holds p_H(h). 
  
  The second issue is more important. When N is large, instantiating or even allocating memory for the truncated characteristic tensor is a challenge, because its size grows exponentially with N. Fortunately, there is a way around this problem. The main idea is that instead of estimating the characteristic tensor of all N variables, we may instead estimate the characteristic tensors of subsets of variables, such as triples, which partially share variables with other triples. The key observation that enables this approach is that the marginal characteristic function of any subset of random variables is also a constrained complex CPD model that inherits parameters from the grand characteristic tensor. Marginalizing with respect to the n'-th random variable, we have that  
  
      Φ (k_1,…,k_n'=0,…,k_N)    =∑_h=1^F ∏_n=1
      n≠ n'^NΦ_X_n|H[k_n]Φ_X_n'|H[0]_=1
         =∑_h=1^F ∏_n=1
      n≠ n'^NΦ_X_n|H[k_n].
  
  Thus, a characteristic function of any subset  of  three random  variables X_i,X_j,X_ℓ (triples) can  be written as a third-order tensor, Φ_ijℓ, of rank F. These sub-tensors can be jointly decomposed in a coupled fashion to obtain the sought factors that allow synthesizing the big characteristic tensor, or decomposed independently and `stitched' later under more stringent conditions. Either way, we beat the curse of dimensionality for low-enough model ranks. In addition to affording significant computational and memory reduction, unlike neural network based methods, the above approach allows us to work with fewer and even missing data during the training phase, i.e., only having access to incomplete realizations of the random vector of interest. We estimate lower-order characteristic function values from only those realizations that all three random variables in a given triple appear together. 
  
  In earlier work, we proposed a similar approach for the categorical case where every random variable is finite-alphabet and the task is to estimate the joint probability mass function (PMF) <cit.>. There we  showed  that every joint PMF of a finite-alphabet random vector can be represented by a naïve Bayes model with a finite number of latent states (rank). If the rank is low, the high dimensional joint PMF is almost surely identifiable from lower-order marginals – which is reminiscent of Kolmogorov extension. 
  
  In case of continuous random variables, however, the joint PDF can no longer be directly represented by a tensor. One possible solution could be discretization, but this unavoidably leads to discretization error. In this work, what we show is that we can approximately represent any smooth joint PDF (and evaluate it at any point) using a low-rank tensor in the characteristic function domain, thereby avoiding discretization loss altogether. 
  
  Our joint PDF model enables easy computation of any marginal or conditional density of subsets of variables of X. Using the conditional expectation, the response variable, taken without loss of generality to be the last variable X_N, can be estimated in the following way (see detailed derivation in Section  <ref>.).
  
      E [  X_N|X_1, …, X_N-1]= 1/c_1∑_h=1^F λ(h)∏_n=1^N-1∑_k_n=-K^K
             𝐀_n(k_n,h)e^-j2π k_n x_n∑_k_N=-K^Kc_2,k_N𝐀_N(k_N,h) 
         where  c_1=∑_h=1^F λ(h)∏_n=1^N-1∑_k_n=-K^K𝐀_n(k_n,h)e^-j2π k_n x_n,
         and  c_2,k_N=e^-j2π k_N/-j2π k_N+1-e^-j2π k_N/[-j2π k_N]^2.
  
  
  One of the very appealing properties of the proposed approach is that it is a generative model that affords easy sampling. According to Equation (<ref>), a sample of the multivariate distribution can be generated by first drawing H according to p_H and then independently drawing samples for each variable X_n from the conditional PDF f_X_n|H. The resulting generative model can be visualized in Figure <ref>.
  
  
  
   §.§ Algorithm: Coupled Tensor Factorization
  
  
  
  We formulate the problem as a coupled complex tensor factorization problem and propose a Block Coordinate Descent algorithm for recovering the latent factors
  of the CPD model representing the joint CF. Then, we only need to invert each conditional CF and synthesize the joint PDF. We refer to this approach as Low-Rank Characteristic Function based Density Estimation (LRCF-DE). 
  
  We begin by defining the following coupled tensor factorization problem
  
      min_λ, 𝐀_1,…,𝐀_N   ∑_i∑_j>i∑_ℓ>jΦ_ijℓ- [[λ, 𝐀_i,𝐀_j, 𝐀_ℓ ]]  _F^2               
      subject to     λ≥0, 1^Tλ = 1, 
         𝐀_n(K+1,:)=1^T,   n=1,…, N.
  
  Each lower-dimensional joint CF of triples,   Φ_ijℓ, can be computed directly from the observed data via sample averaging according to equation (<ref>).
  The formulated optimization problem (<ref>) is non-convex and NP-hard. However it becomes convex with respect to each variable if we fix the remaining ones and can be handled using alternating optimization. By using the mode-1 matrix unfolding of each tensor Φ_ijℓ, the optimization problem with respect to 𝐀_i becomes
  
      min_𝐀_i  ∑_j ≠ i∑_ℓ≠ i, ℓ > jΦ_ijℓ^(1) - (𝐀_ℓ⊙𝐀_j) diag (λ) 𝐀_i^T
      _F^2  
         subject to  𝐀_i(K+1,:)=1^T.
  
  The exact update for each factor 𝐀_i can be computed as
  
      𝐀_i←𝐆_i^-1𝐕_i,
  
  where
  
      𝐆_i=(λλ^T) ⊛∑_j ≠ i∑_ℓ≠ i, ℓ > j𝐐_ℓ j^H𝐐_ℓ j, 
         𝐕_i= diag(λ) ∑_j ≠ i∑_ℓ≠ i, ℓ > j𝐐_ℓ j^HΦ_ijℓ^(1),  
         𝐐_ℓ j = 𝐀_ℓ⊙𝐀_j.
  
  For each update, the row of 𝐀_i that corresponds to zero frequency is removed and updating 𝐀_i becomes an unconstrained complex least squares problem. A vector of ones is appended at the same row index after each update 𝐀_i. Due to role symmetry the same form holds for each factor 𝐀_n. 
  
  Now, for the λ-update we solve the following optimization problem
  
  
  
      min_λ  ∑_i∑_j>i∑_ℓ>jvec(Φ_ijℓ)- (𝐀_ℓ⊙𝐀_j⊙𝐀_i) λ_F^2 
         subject to   λ≥0, 1^Tλ = 1.
  
  The optimization problem (<ref>) is a least squares problem with a probability simplex constraint. We use an ADMM algorithm to tackle it. Towards this end, we reformulate the optimization problem by introducing an auxiliary variable λ̂ and rewrite the problem equivalently as
  
      min_λ, λ̂   f  (λ̂ )+r(λ) 
         subject to λ̂^T=λ,
  
  where, f(λ̂)= ∑_i∑_j>i∑_ℓ>jvec(Φ_ij ℓ)- (𝐀_ℓ⊙𝐀_j⊙𝐀_i) λ̂_F^2 
  and 
  r(λ) 
  is the indicator function for the probability simplex. C={λ|λ≥0, 1^Tλ= 1},
  
      r(λ) =
      
      0,   λ∈ C
      ∞,   λ∉C.
  
  At each iteration τ, we perform the following updates
  
      λ̂^τ+1←(𝐆+ρ𝐈)^-1(𝐕+ρ(λ^τ+𝐮^τ)) 
         λ^τ+1←𝒫_C(λ^τ-λ̂^τ+1 + 𝐮^τ) 
         𝐮^τ+1←𝐮^τ+λ^τ+1-λ̂^τ+1,
      where
         𝐆= ∑_i∑_j>i∑_ℓ>j𝐐_ℓ j i^H𝐐_ℓ ji,
         𝐕= ∑_i∑_j>i∑_ℓ>j𝐐_ℓ j i^H vec(Φ),
         𝐐_ℓ ji = 𝐀_ℓ⊙𝐀_j⊙𝐀_i.
  
  𝒫_C(𝐲) denotes the projection operator onto the convex set C – it computes the Euclidean projection of the real part of a point 𝐲= [y_1,…, y_F]^T∈ℂ^F onto the probability simplex
  
      min_𝐱∈ℝ^F   1/2𝐱-(y)  _F^2 
      subject to      𝐱≥0, 1^T𝐱 = 1,
   
  using the method described in <cit.>. The overall procedure is described in Algorithm
  <ref>.
  
  As the final step, the factors are assembled from the triples and the joint CF over all variables is synthesized as Φ = [[ λ,𝐀_1,…,𝐀_N]]. Given, the model of the joint CF, the corresponding joint PDF model can be recovered at any point as
  
      f_ X(𝐱)=∑_h=1^F λ(h)∏_n=1^N∑_k_n=-K^K𝐀_n(k_n,h)e^-j2π k_n x_n.
  
  
  
  
  
  
  
  
  
  § EXPERIMENTS
  
  
  
   §.§ Low Dimensional Toy Data
  
  
  We first show motivating results from modeling low dimensional datasets and showcase the expressivity of the proposed model as well as the significance of each parameter. We begin by modeling  M=2000 samples from the Weight-Height dataset. In Figure <ref>,  we present M'=1500 synthetic samples obtained from the proposed model for different smoothing parameters K∈[1,2,3,4,5] and ranks F ∈[2,4,6,8,10]. By judiciously selecting the parameter search space, our approach yields an expressive generative model that can well-represent the data.
  
  Following the same procedure, we now visualize M=2000 samples from the 2-D Moons and Circles datasets. We fix the number of smoothing coefficients K, K=11, and visualize synthetic M'=1500 samples obtained from the proposed model for different approximation ranks F ∈[1,2,3,4,6] in Figure <ref>. The results show that our model is able to capture complex structures and properties of the data for an appropriate choice of rank F.  
  
  
  
  
  
  
   §.§ Real Data
  
  
  We test the proposed approach on datasets (see a brief description of the datasets in Table <ref>) obtained from the UCI machine learning repository <cit.>.
  
  
  
  
  
  For each dataset we randomly hide 20% of data (testing set) and consider the remaining entries as observed information (training set). The parameters, which include the tensor rank F and the smoothing parameter K, are chosen using cross-validation. We use  20% of the training data as validation data, where we seek to find the optimal parameter values maximizing the average log-likelihood of the validation samples. Once the hyperparameters are chosen, we train the model using all the training data (including the validation data) and measure its performance
  on the testing set. We compare our approach against standard baselines described in section <ref>.
  
  Evaluating the quality of density models is an open and difficult problem <cit.>. Following the approach in <cit.>, we calculate and report the average log-likelihood of unseen data samples (testing set), further averaged over 5 random data splits. The results are shown in Table <ref>.
  LRCF-DE has a higher average test sample log likelihood on almost all datasets. Overall, we observe that our method outperforms the baselines in 4 datasets and is comparable to the winning method in the remaining ones.
  
   Following the derivation in  Equation (<ref>), we test the proposed model in several regression tasks. We evaluate and report the Mean Absolute Error (MAE) in estimating X_N for the unseen data samples in Table <ref> and additional results for multi-output regression are presented in Table <ref>. Overall, we observe that LRCF-DE outperforms the baselines on almost all datasets, and performs comparable to the winning method in the remaining ones. 
   
   We have to stress again the fact that neural network based density estimation methods evaluate multivariate densities point-wise. These methods cannot impute more than a few missing elements in the input as grid search becomes combinatorial. Due to the interpretation of the approximation of the sought density as a finite mixture of separable densities and the coupled tensor factorization approach, our method allows us to easily work with missing data during both training and testing. Here, we showcase the results of LRCF-DE against MAF for simultaneously predicting the last two random variables of each dataset given the remaining ones.
  
  
  
  
  As our last experiment, we train LRCF-DE to learn the joint distribution of grayscale images from the USPS dataset <cit.>, which contains 9298 images of handwritten digits of size 16 × 16 → N=256. The number of examples for each digit is shown in Table <ref>. We sample from the resulting 256-dimensional model, and provide visualization of the generated data. We fix the tensor rank to F=8 and the smoothing parameter to K=15, and draw 8 random samples of each digit (class). The resulting samples are shown in Figure <ref>, and they are very pleasing – in light of the fact that our model is “agnostic”: designed for general-purpose density estimation, not specifically for realistic-looking image synthesis. It is possible to incorporate image modeling domain knowledge in the design of LRCF-DE (such as correlation between adjacent pixel values), but this is beyond the scope of this paper. 
  
  
  
  § CONCLUSIONS
  
  In this work, we have revisited the classic problem of non-parametric density estimation from a fresh perspective – through the lens of complex Fourier series approximation and tensor modeling, leading to a low-rank characteristic function approach. We showed that any compactly supported density can be well-approximated by a finite characteristic tensor of leading complex Fourier coefficients as long as the coefficients decay sufficiently fast. We posed density  estimation as a constrained (coupled) tensor factorization problem and proposed a Block Coordinate Descent algorithm, which under certain conditions enables learning the true data-generating distribution. Results on real data have demonstrated the utility and promise of this novel approach compared to both standard and recent density estimation techniques.
  
  
  
  
  
  
  
  
  § APPENDIX
  
  
  
  
   §.§ Bias, variance, consistency of the empirical characteristic function
  
  
  In this appendix we summarize important properties of empirical characteristic functions as statistical estimators of the corresponding characteristic functions. We refer the reader to <cit.> for proofs and additional results.  
  
  By linearity of expectation, it is easy to see that the empirical characteristic function is an unbiased estimator of the corresponding characteristic function, i.e., 
      E[ Φ_ X(ν) ] = Φ_ X(ν),
  
  for all ν and M ≥ 1. For the remainder of this section, we assume that {𝐱_m }_m=1^M is i.i.d. in m. The variance of the  empirical characteristic function estimate can be shown <cit.> to be    
  
      Var[Φ_ X(ν)]   =  E[|Φ_ X(ν) -E[ Φ_ X(ν) ]|^2 ]
         =  E[|Φ_ X(ν)-Φ_ X(ν)|^2 ]
         = 1/M(1-|Φ_ X(ν)|^2).
  
  Note that 0≤|Φ_ X(ν)|≤1, and therefore Var[Φ_ X(ν)]≤1/M. It follows that 
  
      lim_M →∞ E[|Φ_ X(ν)-Φ_ X(ν)|^2 ] = 0,
  
  i.e., for any fixed ν, Φ_ X(ν) converges to Φ_ X(ν) in the mean-squared sense. By the strong law of large numbers, it also converges almost surely for any fixed ν. Furthermore, for any fixed positive T<∞
  
      lim_M →∞sup _ν_2 ≤T|Φ_ X(ν)-Φ_ X(ν)|=0,
  
  almost surely. It can also be shown <cit.> that for any increasing sequence T_M such that lim_M →∞log(T_M)/M=0, it holds 
  
      lim_M →∞sup _ν_2 ≤T_M|Φ_ X(ν)-Φ_ X(ν)| = 0,
  
  almost surely. In our context, we only use a sampled and truncated version of the characteristic function (corresponding to a truncated multivariate Fourier series), hence T is always finite – we do not need the latter result. 
  
  It is also worth noting that the covariance of different samples of the empirical characteristic function (corresponding to different values of ν) goes to zero ∼1/M, and so does the covariance of its real and imaginary parts. As a result, for large M, the errors in the different elements of the characteristic tensor are approximately uncorrelated, with uncorrelated real and imaginary parts. This suggests that when we fit a model to the empirical characteristic function, it makes sense to use a least squares approach. Another motivation for this comes from Parseval's theorem: minimizing integrated squared error in the Fourier domain corresponds to minimizing integrated squared error between the corresponding multivariate distributions. This is true in particular when we limit the support of the distribution to a hypercube and use the samples of the characteristic function that correspond to the multivariate Fourier series, thereby replacing the multivariate integral in the Fourier domain by a multivariate sum. 
  
  
  
   §.§ Low-rank denoising: reduction of the mean squared error
  
  
  Our use of a low-rank model in the characteristic tensor domain is primarily motivated by the need to avoid the “curse of dimensionality”: using a rank-F model with 2K+1 Fourier coefficients per mode parametrizes the whole N-dimensional multivariate density using just FN(2K+1) coefficients, and avoids instantiating and storing a tensor of size (2K+1)^N, which is close to impossible even for moderate N. However, there is also a variance benefit that comes from this low-rank parametrization. We know from <cit.> that for a square L × L matrix of rank F observed in zero-mean white noise of variance σ^2, low-rank denoising attains mean squared error c LF σ^2 asymptotically in L, where c is a small constant. In practice the asymptotics kick in even for relatively small L <cit.>. Contrast this to the raw L^2 σ^2 if one does not use the low-rank property. 
  
  For an N-way tensor of rank F, assume N is even, K_n=K, ∀ N (for simplicity of exposition), and “unfold” the characteristic tensor into a K^N/2× K^N/2 matrix. In practice we use F far less than K^N/2, and thus the resulting matrix will be very low rank. Invoking <cit.>, low-rank tensor modeling will yield a reduction in mean squared error by a factor of at least F/K^N/2. We say at least, because this low-rank matrix structure is implied but does not imply low-rank tensor structure, which is much stronger. Note also that mean squared error in the characteristic tensor domain translates to mean squared error between the corresponding distributions, by virtue of Parseval's theorem.  
  
  
  
   §.§ Derivation of (<ref>)
  
  In this appendix we present the derivation of (<ref>), which is used to solve  regression  tasks.
  
      E   [  X_N|X_1, …, X_N-1] 
         = ∫_0^1 x_N f_X_N|X_1 ,…, X_N-1(x_N|x_1, …, x_N-1) dx_N 
         = ∫_0^1 x_N f_X_1,…,X_N(x_1,…, x_N)/f_X_1,…,X_N-1(x_1,…,x_N-1) dx_N
         = 1/c_1∫_0^1 x_N f_X_1,…,X_N(x_1,…,x_N) dx_N  
         = 1/c_1∫_0^1 x_N ∑_h=1^F λ(h)∏_n=1^N∑_k_n=-K^K𝐀_n(k_n,h)e^-j2π k_n x_ndx_N  
         = 1/c_1∑_h=1^F λ(h)∏_n=1^N-1∑_k_n=-K^K𝐀_n(k_n,h)e^-j2π k_n x_n
                         ∑_k_N=-K^K𝐀_N(k_N,h)∫_0^1 x_N e^-j2π k_N x_Ndx_N  
         = 1/c_1∑_h=1^F λ(h)∏_n=1^N-1∑_k_n=-K^K𝐀_n(k_n,h)e^-j2π k_n x_n
                                 ∑_k_N=-K^Kc_2,k_N𝐀_N(k_N,h),
  
  where 
  
      c_1   =f_X_1,…,X_N-1(x_1,…,x_N-1)
         =∑_h=1^F λ(h)∏_n=1^N-1∑_k_n=-K^K𝐀_n(k_n,h)e^-j2π k_n x_n,
  
  and 
  
      c_2,k_N   =e^-j2π k_N/-j2π k_N+1-e^-j2π k_N/[-j2π k_N]^2.
  
  
  
  
  IEEEtran
  
  
  
  
  
  
  