PSD Representations for Effective Probability Models
Year: 2022

Abstract:
Finding a good way to model probability densities is key to probabilistic inference. An ideal model should be able to concisely approximate any probability while being also compatible with two main operations: multiplications of two models (product rule) and marginalization with respect to a subset of the random variables (sum rule). In this work, we show that a recently proposed class of positive semidefinite (PSD) models for non-negative functions is particularly suited to this end. In particular, we characterize both approximation and generalization capabilities of PSD models, showing that they enjoy strong theoretical guarantees. Moreover, we show that we can perform efficiently both sum and product rule in closed form via matrix operations, enjoying the same versatility of mixture models. Our results open the way to applications of PSD models to density estimation, decision theory and inference.

heading: 1 Introduction
text: Modeling probability distributions is a key task for many applications in machine learning [17, 5]. To this end, several strategies have been proposed in the literature, such as adopting mixture models (e.g. Gaussian mixtures) [5], exponential models [32], implicit generative models [11, 13] or kernel (conditional) mean embeddings [16]. An ideal probabilistic model should have two main features: i) efficiently perform key operations for probabilistic inference, such as sum rule (i.e. marginalization) and product rule [5] and, ii) concisely approximate a large class of probabilities. Finding models that satisfy these two conditions is challenging and current methods tend to tackle only one of the two. Exponential and implicit generative models have typically strong approximation properties (see e.g. [32, 29]) but cannot easily perform operations such as marginalization. On the contrary, mixture models are designed to efficient integrate and multiply probabilities, but tend to require a large number of components to approximate complicated distributions.
In principle, mixture models would offer an appealing strategy to model probability densities since they allow for efficient computations when performing key operations such as sum and product rule. However, these advantages in terms of computations, come as a disadvantage in terms of expressiveness: even though mixture models are universal approximators (namely they can approximate arbitrarily well any probability density), they require a significant number n of observations and of components to do that. Indeed it is known that models that are non-negative mixtures of non-negative components lead to learning rates that suffer from the curse of dimensionality. For example, when approximating probabilities on Rd, kernel density estimation (KDE) [21] with non-negative components has rates as slow as n 2/(4+d) (see, e.g. [38, page 100]), and cannot improve with the regularity (smoothness) of the target probability (see e.g. [9, Thm. 16.1] for impossibility results in the one dimensional case).
In the past decades, this limitation has been overcome by removing either the non-negativity of the weights of the mixture (leading to RBF networks [27][35]), or the non-negativity of the components
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
used in the mixture (leading to KDE with oscillating kernel [35]), or both. On the positive side this allows to achieve a learning rate of n /(2 +d), for -times differentiable densities on Rd. Note that such rate is minimax optimal [10] and overcomes the curse of dimensionality when d. Additionally the resulting model is very concise, since only m = O(nd/ ) centers are necessary to achieve the optimal rate [36]. However, on the negative side, the resulting model is not a probability, since it may attain negative values, and so it cannot be used where a proper probability is required.
In this paper, we show that the positive semidefinite (PSD) models, recently proposed in [14], offer a way out of this dichotomy. By construction PSD models generalize mixture models by allowing also for negative weights, while still guaranteeing that the resulting function is non-negative everywhere. Here, we prove that they get the best of both worlds: expressivity (optimal learning and approximation rates with concise models) and flexibility (exact and efficient sum and product rule).
Contributions. The main contributions of this paper are:
i) Showing that PSD models can perform exact sum and product rules in terms of efficient matrix operations (Sec. 2.1).
ii) Characterizing the approximation and learning properties of PSD models with respect to a large family of probability densities (Sec. 3).
iii) Providing a “compression” method to control the number of components of a PSD model (introduced in Sec. 2.2 and analyzed in Sec. 3.3).
iv) Discussing a number of possible applications of PSD models (Sec. 4).
Summary of the results. Table 1 summarizes all the desirable properties that a probability model should have and which of them are satisfied by state-of-the-art estimators. Among these we have: non-negativity, the model should be point-wise non-negative; sum and product rule, the model should allow for efficient computation of such operations between probabilities; concise approximation, a model with the optimal number of components m = O(" d/ ) is enough to approximate with error " a non-negative function in a wide family of smooth -times differentiable functions [36]; optimal learning, the model achieves optimal learning rates of order n /(2 +d) when learning the unknown target probability from i.i.d. samples [35]; efficient sampling the model allows to extract i.i.d. samples without incurring in the curse of dimensionality in terms of computational complexity. We note that for PSD models, this last property has been recently investigated in [15]. We consider under the umbrella of linear models all the models which corresponds to a mixture P m
j=1 wjfj(x) of functions fj : X ! R, with either wj 0, 8j = 1, . . . ,m or fj 0, 8j = 1, . . . ,m or both. We denote by mixture models the mixture models defined as above, where wj 0, fj 0, 8j = 1, . . . ,m. By mean embeddings, we denote the kernel mean embedding estimator framework from [30] see also [16]). With exponential models we refer to models of the form exp( P m
j=1 wjfj(x)) with wj 2 R, fj : X ! R [32]. With PSD models we refer to the framework introduced in [14] and studied in this paper. We note that this work contributes in showing that the latter are the only probability models to satisfy all requirements (we recall that non-negativity and efficient sampling have been shown respectively in [14] and [15]).
Notation. We denote by Rd++ the space vectors in Rd with positive entries, Rn⇥d the space of n ⇥ d matrices, Sn+ = S+(Rn) the space of positive semidefinite n ⇥ n matrices. Given a vector ⌘ 2 Rd, we denote diag(⌘) 2 Rd⇥d the diagonal matrix associated to ⌘. We denote by A B
and A ⌦ B respectively the entry-wise and Kronecker product between two matrices A and B. We denote by kAkF , kAk, det(A), vec(A) and A> respectively the Frobenius norm, the operator norm (i.e. maximum singular value), the determinant, the (column-wise) vectorization of a matrix and the (conjugate) transpose of A. With some abuse of notation, where clear from context we write element-wise products and division of vectors u, v 2 Rd as uv, u/v. Given two matrices X 2 Rn⇥d1 , Y 2 Rn⇥d2 with same number of rows, we denote by [X,Y ] 2 Rn⇥(d1+d2) their concatenation row-wise. The term 1n 2 Rn denotes the vector with all entries equal to 1.

heading: 2 PSD Models
text: Following [14], in this work we consider the family of positive semi-definite (PSD) models, namely non-negative functions parametrized by a feature map : X ! H from an input space X to a suitable feature space H (a separable Hilbert space e.g. Rq) and a linear operator M 2 S+(H), of the form
f(x ; M, ) = (x)>M (x). (1)
PSD models offer a general way to parametrize non-negative functions (since M is positive semidefinite, f(x ; M, ) 0 for any x 2 X ) and enjoy several additional appealing properties discussed in the following. In this work, we will focus on a special family of models of the form (1) to parametrize probability densities over X = Rd. In particular, we will consider the case where: i) = ⌘ : Rd ! H⌘ is a feature map associated to the Gaussian kernel [26] k⌘(x, x0) = ⌘(x)> ⌘(x) = e (x x
0)>diag(⌘)(x x0), with ⌘ 2 Rd++ and, ii) the operator M lives in the span of (x1), . . . , (xn) for a given set of points (xi)ni=1, namely there exists A 2 Sn+ such that M = P ij Aij (xi) (xj)>. We define a Gaussian PSD model by specializing the definition in (1) as
f(x ; A,X, ⌘) = nX
i,j=1
Ai,jk⌘(xi, x)k⌘(xj , x), 8x 2 Rd (2)
in terms of the coefficient matrix A 2 Sn+, the base points matrix X 2 Rn⇥d, whose i-th row corresponds to the point xi for each i = 1, . . . , n and kernel parameter ⌘. In the following, given two base point matrices X 2 Rn⇥d and X 0 2 Rm⇥d, we denote by KX,X0,⌘ 2 Rn⇥m the kernel matrix with entries (KX,X0,⌘)ij = k⌘(xi, x0j) where xi, x0j are the i-th and j-th rows of X,X 0 respectively. When clear from context, in the following we will refer to Gaussian PSD models as PSD models. Remark 1 (PSD models generalize Mixture models). Mixture models (a mixture of Gaussian distributions) are a special case of PSD models. Let A = diag(a) be a diagonal matrix of n positive weights a 2 Rn++. We have f(x ; A,X, ⌘/2) = P n i=1 Aiik⌘/2(xi, x) 2 = P n
i=1 aik⌘(xi, x). Remark 2 (PSD models allow negative weights). From (2), we immediately see that PSD models generalize mixture models by allowing also for negative weights: e.g., f(· ; A,X, ⌘) with A = (1, 12 ; 1 2 , 1 4 ) 2 S 2 +, ⌘ = 1, X = (x1;x2) and x1 = 0, x2 = 2, corresponding to f(x ; A,X, ⌘) = e 2x 2
+ 14e 2(x 2)2 1 e e 2(x 1) 2 , i.e. a mixture of Gaussians with also negative weights.

heading: 2.1 Operations with PSD models
text: In Sec. 3 we will show that PSD models can approximate a wide class of probability densities, significantly outperforming mixture models. Here we show that this improvement does not come at the expenses of computations. In particular, we show that PSD models enjoy the same flexibility of mixture models: i) they are closed with respect to key operations such as marginalization and multiplication and ii) these operations can be performed efficiently in terms of matrix sums/products. The derivation of the results reported in the following is provided in Appendix F. They follow from well-known properties of the Gaussian function.
Evaluation. Evaluating a PSD model in a point x0 2 X corresponds to f(x = x0 ; A,X, ⌘) = K>
X,x0,⌘ AKX,x0,⌘ . Moreover, partially evaluating a PSD in one variable yields
f(x, y = y0 ; A, [X,Y ], (⌘1, ⌘2)) = f(x ; B,X, ⌘1) with B = A (KY,y0,⌘2K > Y,y0,⌘2 ). (3)
Note that f(x ; B,X, ⌘1) is still a PSD model since B is positive semidefinite.
Sum Rule (Marginalization and Integration). The integral of a PSD model can be computed as Z
f(x ; A,X, ⌘) dx = c2⌘ Tr(A KX,X, ⌘2 ) with c⌘ = Z k⌘(0, x) dx, (4)
where c⌘ = ⇡d/2 det(diag(⌘)) 1/2. This is particularly useful to model probabiliy densities with PSD models. Let Z = R f(x ; A,X, ⌘)dx, then the function f(x ; A/Z,X, ⌘) = 1
Z f(x ; A,X, ⌘)
is a probability density. Integrating only one variable of a PSD model we obtain the sum rule.
Proposition 1 (Sum Rule – Marginalization). Let X 2 Rn⇥d, Y 2 Rn⇥d0 , A 2 S+(Rn) and ⌘ 2 Rd++, ⌘0 2 Rd 0
++. Then, the following integral is a PSD model Z
f(x, y ; A, [X,Y ], (⌘, ⌘0)) dx = f(y ; B, Y, ⌘0), with B = c2⌘ A KX,X, ⌘2 , (5)
The result above shows that we can efficiently marginalize a PSD model with respect to one variable by means of an entry-wise multiplication between two matrices. Remark 3 (Integration and marginalization on the hypercube). The integrals in (4) and (5) can be performed also on when X is a hypercube H = Q d
t=1[at, bt] rather than the entire space Rd. This leads to a closed form, where the matrix KX,X, ⌘2 is replaced by a suitable KX,X, ⌘2 ,H that can be computed with same number of operations (the full form of such matrix is reported in Appendix F).
Product Rule (Multiplication). Multiplying two probabilities is key to several applications in probabilistic inference [5]. The family of PSD models is closed with respect to this operation. Proposition 2 (Multiplication). Let X 2 Rn⇥d1 , Y 2 Rn⇥d2 , Y 0 2 Rm⇥d2 , Z 2 Rm⇥d3 , A 2 Sn+, B 2 Sm+ and ⌘1 2 Rd1++, ⌘2, ⌘02 2 Rd2++, ⌘3 2 Rd3++. Then
f(x, y ; A, [X,Y ], (⌘1, ⌘2))f(y, z ; B, [Y 0, Z], (⌘02, ⌘3)) = f(x, y, z ; C,W, ⌘), (6)
is a PSD model, where C = (A ⌦ B) vec(KY,Y 0,e⌘2)vec(KY,Y 0,e⌘2)> , with ⌘̃2 =
⌘2⌘ 0 2
⌘2+⌘02 , base
matrix W = [X ⌦ 1m, ⌘2⌘2+⌘02Y ⌦ 1m + ⌘ 0 2 ⌘2+⌘02 1n ⌦ Y 0, 1n ⌦ Z] and ⌘ = ⌘1, ⌘2 + ⌘02, ⌘3 .
We note that, despite the heavy notation, multiplying two PSD is performed via simple operations such as tensor and entry-wise product between matrices. In particular, we note that X ⌦ 1m 2 Rnm⇥d1 and 1n ⌦ Z 2 Rnm⇥d3 correspond respectively to the nm⇥ d1-matrix containing m copies of each row of X and the nm⇥ d3-matrix containing n copies of Z. Finally, ⌘Y ⌦ 1m + ⌘01n ⌦ Y is the nm⇥ d2 base point matrix containing all possible sums (⌘yi + ⌘0y0j) n,m i,j=1 of points from Y and Y 0.
Reduction. As observed above, when performing operations with PSD models, the resulting base point matrix might be of the form X ⌦ 1m (e.g. if we couple the multiplication in (6) with marginalization as in a Markov transition, see the corollary below). In these cases we can reduce the PSD model to have only n base points (rather than nm), as follows
f(x ; A,X ⌦ 1m, ⌘) = f(x ; B,X, ⌘) with B = (In ⌦ 1>m)A(In ⌦ 1m), (7)
where A 2 Snm+ and In 2 Rn⇥n is the n⇥ n identity matrix. The reduction operation is useful to avoid the dimensionality of the PSD model grow unnecessarily. This is for instance the case of a Markov transition. Corollary 3 (Markov Transition). Let X 2 Rn⇥d1 , Y 2 Rn⇥d2 , Y 0 2 Rm⇥d2 , A 2 Sn+, B 2 Sm+ and ⌘1 2 Rd1++, ⌘2, ⌘02 2 Rd2++. Then
Z f(x, y ; A, [X,Y ], (⌘1, ⌘2))f(y ; B, Y 0, ⌘02) dy = f(x ; C,X, ⌘1), (8)
with C 2 Sn+ obtained by applying in order, Prop. 2, Prop. 1 and reduction (7).
We remark that the result of a Markov transition retains the same base point matrix X and parameters ⌘ of the transition kernel. This is thanks to the reduction operation in (7), which avoids the resulting matrix C to be nm⇥nm. This fact is particularly useful in applications that require multiple Markov transitions, such as in hidden Markov models (see also Sec. 4).

heading: 2.2 Compression of a PSD model
text: From Prop. 2 we note that the multiplication operation can rapidly yield a large number of base points and thus incur in high computational complexity in some settings. It might therefore be useful to have a method to reduce the number of base points while retaining essentially the same model. To this purpose, here we propose a dimensionality reduction strategy. In particular, given a set of points x̃1, . . . , x̃m 2 Rd, we leverage the representation of the PSD model in terms of reproducing kernel Hilbert spaces [3] to use powerful sketching techniques as Nyström projection (see, e.g., [40]), to project the PSD model on a new PSD model now based only on the new points. Given A 2 Sn+, X 2 Rn⇥d, ⌘ 2 Rd++. Let X̃ 2 Rm⇥d be the base point matrix whose j-rows corresponds to the point x̃j . The compression of f(· ; A,X, ⌘) corresponds to
f(x ; Ã, X̃, ⌘) with Ã = BAB> 2 Sm+ , B = K 1eX, eX,⌘K eX,X,⌘ 2 R m⇥n, (9)
which it is still a PSD model since the matrix BAB> 2 Sm+ . We not that even with a rather simple strategy to choose the new base points x̃1, . . . , x̃m – such as uniform sampling – compression is an effective tool to reduce the computational complexity of costly operations. In particular, in Sec. 3.3 we show that a compressed model with only O(t polylog(1/")) centers (instead of t2) can "-approximate the product of two PSD models with t points each.

heading: 3 Representation power of PSD models
text: In this section we study the theoretical properties of Gaussian PSD models. We start by showing that they admit concise approximations of the target density (in the sense discussed in the introduction to this paper and of Table 1). We then proceed to studying the setting in which we aim to learn an unknown probability from i.i.d. samples. We conclude the section by characterizing the approximation properties of the compression operation introduced in Sec. 2.2.

heading: 3.1 Approximation properties of Gaussian PSD models
text: We start the section recalling that Gaussian PSD models are universal approximators for probability densities. In particular, the following result restates [Thm. 2 14] for the case of probabilities. Proposition 4 (Universal consistency – Thm. 2 in [14]). The Gaussian PSD family is a universal approximator for probabilities that admit a density.
The result above is not surprising since Gaussian PSD models generalize classical Gaussian mixtures (see Remark 1), which are known to be universal [5]. We now introduce a mild assumption that will enable us to complement Prop. 4 with approximation and learning results. In the rest of the section we assume that X = ( 1, 1)d (or more generally an open bounded subset of Rd with Lipschitz boundary). Here L1(X ) and L2(X ) denote respectively the space of essentially bounded and square-integrable functions over X , while W 2 (X ) denotes the Sobolev space of functions whose weak derivatives up to order are square-integrable on X (see [1] or Appendix A for more details). Assumption 1. Let > 0, q 2 N. There exists f1, . . . , fq 2 W 2 (X )\L1(X ), such that the density p : X ! R satisfies
p(x) = qX
j=1
fj(x) 2, 8x 2 X . (10)
The assumption above is quite general and satisfied by a wide family of probabilities, as discussed in the following proposition. The proof is reported in Appendix D.1 Proposition 5 (Generality of Assumption 1). The assumption above is satisfied by
(a) any probability density p 2 W 2 (X ) \ L 1(X ) and strictly positive on [ 1, 1]d,
(b) any exponential model p(x) = e v(x) with v 2 W 2 (X ) \ L 1(X ),
(c) any mixture model of Gaussians or, more generally, of exponential models from (b), (d) any p that is + 2-times differentiable on [ 1, 1]d, with a finite set of zeros, all in ( 1, 1)d,
and with positive definite Hessian in each zero. E.g. p(x) / x2e x 2 .
Moreover if p is -times differentiable over [ 1, 1]d, then it belongs to W 2 (X ) \ L 1(X ).
We note that in principle the class C ,d = W 2 (X )\L1(X ) is larger than the Nikolskii class usually considered in density estimation [10] when 2 N, but Assumption 1 imposes a restriction on it. However, Prop. 5 shows that the restriction imposed by Assumption 1 is quite mild, since it the resulting space includes all probabilities that are -times differentiable and strictly positive (it also contains probabilities that are not strictly positive but have some zeros, see [24]). We can now proceed to the main result of this work, which characterizes the approximation capabilities of PSD models. Theorem 6 (Conciseness of PSD Approximation). Let p satisfy Assumption 1. Let " > 0. There exists a Gaussian PSD model of dimension m 2 N, i.e., p̂m(x) = f(x ; Am, Xm, ⌘m), with Am 2 Sm+ and Xm 2 Rm⇥d and ⌘m 2 Rd++, such that
kp p̂mkL2(X )  ", with m = O(" d/ (log 1" ) d/2). (11)
The proof of Thm. 6 is reported in Appendix D.3. According to the result, the number of base points needed for a PSD model to approximate a density up to precision " depends on its smoothness (the smoother the better) and matches the bound m = " d/ that is optimal for function interpolation [18], corresponding to models allowing for negative weights, and is also optimal for convex combinations of oscillating kernels [36].

heading: 3.2 Learning a density with PSD models
text: In this section we study the capabilities of PSD models to estimate a density from n samples. Let X = ( 1, 1)d and let p be a probability on X . Denote by x1, . . . , xn the samples independently and identically distributed according to p, with n 2 N. We consider a Gaussian PSD estimator p̂n,m = f(x ; Â, X̃, ⌘) that is built on top of m additional points x̃1, . . . , x̃m, sampled independenly and uniformly at random in X . In particular, ⌘ 2 Rd++, X̃ 2 Rm⇥d is the base point matrix whose j-th row corresponds to the point x̃j and Â 2 Sm+ is trained as follows
Â = argmin A2Sm+
Z
X
f(x ; A, X̃, ⌘)2dx 2 n
nX
i=1
f(xi ; A, X̃, ⌘) + kK 1/2AK1/2k2 F , (12)
where K = K X̃,X̃,⌘ . Note that the functional is constituted by two parts. The first two elements are an empirical version of kf pk2
L2(X ) modulo a constant independent of f (and so not affecting the optimization problem), since xi are identically distributed according to p and so 1n P n
i=1 f(xi) ⇡R f(x)p(x)dx. The last term is a regularizer and corresponds to kK1/2AK1/2k2
F = Tr(AKAK),
i.e. the Frobenius norm of M = P m
ij=1 Aij ⌘(x̃i) ⌘(x̃j). The problem in (12) corresponds to a quadratic problem with a semidefinite constraint and can be solved using techniques such as Newton method [24] or first order dual methods [14]. We are now ready to state our result. Theorem 7. Let n,m 2 N, > 0, ⌘ 2 Rd++ and p be a density satisfying Assumption 1. With the definitions above, let p̂n,m be the model p̂n,m(x) = f(x ; Â, X̃, ⌘), with Â the minimizer of (12). Let ⌘ = n 2 2 +d 1d and = n 2 +2d 2 +d . When m C 0n d 2 +d (log n)d log(C 00n(log n)), the following holds with probability at least 1 ,
kp p̂n,mkL2(X )  Cn
2 +d (log n)d/2, (13)
where constant C depends only on , d and p and the constants C 0, C 00 depend only on , d.
The proof of Thm. 7 is reported in Appendix E.2. The theorem guarantees that under Assumption 1, Gaussian PSD models can achieve the rate O(n /(2 +d)) – that is optimal for the -times differentiable densities – while admitting a concise representation. Indeed, it needs a number m = O(nd/(2 +d)) of base points, matching the optimal rate in [36]. When d, a model with m = O(n1/3) centers achieves optimal learning rates.

heading: 3.3 The Effect of compression
text: We have seen in the previous section that Gaussian PSD models achieve the optimal learning rates, with concise models. However, we have seen in the operations section that multiplying two PSD
models of m centers leads to a PSD model with m2 centers. Here we study the effect of compression, to show that it is possible to obtain an "-approximation of the product via a compressed model with O(m polylog(1/")) centers. In the following theorem we analyze the effect in terms of the L1 distance on a domain [ 1, 1]d, induced by the compression, when using points taken independently and uniformly at random from the same domain.
Let A 2 Sn+, X 2 Rn⇥d, ⌘ 2 Rd++, we want to study the compressibility of the PSD model p(x) = f(x ; A,X, ⌘). Let X̃ 2 Rm⇥d be the base point matrix whose j-rows corresponds to the point x̃j with x̃1, . . . , x̃m be sampled independently and uniformly at random from [ 1, 1]d. Denote by p̃m(x) the PSD model p̃m(x) = f(x ; Ã, X̃, ⌘) where Ã is the compression of A via (9). We have the following theorem. Theorem 8 (Compression of Gaussian PSD models). Let 2 (0, 1], ⌘+ = max(1,maxi=1,...,d ⌘i). When m satisfies
m O ⇣ ⌘1/2+ log kAkn
"
d log 1 ⌘ , (14)
then the following holds with probability at least 1 ,
|p(x) p̃m(x)|  " 2 + " p p(x), 8x 2 [ 1, 1]d, (15)
The proof of the theorem above is in Appendix C.1. To understand its relevance, let p̂1 be a PSD model trained via (12) on n points sampled from p1 and p̂2 trained from n points sampled from p2, where both p1, p2 satisfy Assumption 1 for the same and m, , ⌘ are chosen as Thm. 7, in particular m = nd/(2 +d) and ⌘ = ⌘+1d, ⌘+ = n2/(2 +d). Consider the model p̂ = p̂1 · p̂2. By construction p̂ has m2 = n2d/(2 +d) centers, since it is the pointwise product of p̂1, p̂2 (see Prop. 2) and approximates p1 · p2 with error " = n /(2 +d) polylog(n), since both p̂1, p̂2 are "-approximators of p1, p2. Instead, by compressing p̂, we obtain an estimator p̄, that according to Thm. 8, achieves error " with a number of center
m0 = O(⌘d/2+ polylog(1/")) = O(n d/(2 +d) polylog(1/")) = O(m polylog(1/")). (16)
Then p̄ approximates p1 · p2 at the optimal rate n /(2 +d), but with a number of centers m0 that is only O(m polylog(n)), instead of m2. This means that p̄ is essentially as good as if we learned it from n samples taken directly from p1 · p2. This renders compression a suitable method to reduce the computational complexity of costly inference operations as the product rule.

heading: 4 Applications
text: PSD models are a strong candidate in a variety of probabilistic settings. On the one hand, they are computationally amenable to performing key operations such as sum and product rules, similarly to mixture models (Sec. 2.1). On the other hand, they are remarkably flexible and can approximate/learn (coincisely) a wide family of target probability densities (Sec. 3). Building on these properties, in this section we consider different possible applications of PSD models in practice.

heading: 4.1 PSD Models for Decision Theory
text: Decision theory problems (see e.g. [5] and references therein) can be formulated as a minimization
✓⇤ = argmin ✓2⇥ L(✓) = Ex⇠p `(✓, x), (17)
where ` is a loss function, ⇥ is the space of target parameters (decisions) and p is the underlying data distribution. When we can sample directly from p – e.g. in supervised or unsupervised learning settings – we can apply methods such as stochastic gradient descent to efficently solve (17). However, in many applications, sampling from p is challenging or computationally unfeasible. This is for instance the case when p has been obtained via inference (e.g. it is the t-th estimate in a hidden Markov model, see Sec. 4.3) or it is fully known but has a highly complex form (e.g. the dynamics of a physical system). In contexts where sampling cannot be performed efficiently, it is advisable to consider alternative approaches. Here we propose a strategy to tackle (17) when p can be modeled (or well-approximated) by a PSD model. Our method hinges on the following result.
Proposition 9. Let p(x) = f(x ; A,X, ⌘) with X 2 Rn⇥d, A 2 Sn+, ⌘ 2 Rd++. Let g : Rd ! R and define cg,⌘(z) = R g(x)e ⌘kx zk 2 dx for any z 2 Rd. Then
Ex⇠p g(x) = Tr( (A KX,X,⌘/2)G ) with Gij = cg,2⌘ xi+xj
2
. (18)
Thanks to Prop. 9 we can readily compute several quantities related to a PSD model such as its mean Ep[x] = X>b with b = (A KX,⌘/2)1n, its covariance or its characteristic function (see Appendix F for the explicit formulas and derivations). However, the result above is particularly useful to tackle the minimization in (17). In particular, since rL(✓) = Ex⇠pr✓`(✓, x), we can use Prop. 9 to directly compute the gradient of the objective function: it is sufficient to know how to evaluate (or approximate) the integral cr`✓,⌘(z) = R r✓`(✓, x)e ⌘kx zk 2
dx for any ✓ 2 ⇥ and z 2 Rd. Then, we can use first order optimization methods, such as gradient descent, to efficiently solve (17). Remarkably, this approach works well also when we approximate p with a PSD model p̂. If ` is convex, since p̂ is non-negative, the resulting L̂(✓) = Ex⇠p̂ `(✓, x) is still a convex functional (see also the discussion on structured prediction in Sec. 4.2). This is not the case if we use more general estimators of p that do not preserve non-negativity.

heading: 4.2 PSD Models for Estimating Conditional Probabilities
text: In supervised learning settings, one is typically interested in solving decision problems of the form min✓2⇥ E(x,y)⇠p `(h✓(x), y) where p is a probability over the joint input-output space X ⇥ Y and h✓ : X ! Y is a function parameterized by ✓. It is well-known [see e.g. 34] that the ideal solution of this problem is the ✓⇤ such that for any x 2 X the function h✓⇤(x) = argminz2Y Ey⇠p(·|x) `(z, y) is the minimizer with respect to z 2 Y of the conditional expectation of `(z, y) given x. This leads to target functions that capture specific properties of p, such as moments. For instance, when ` is the squared loss, h✓⇤(x) = Ey⇠p(·|x)y corresponds to the conditional expectation of y given x, while for ` the absolute value loss, h✓⇤ recovers the conditional median of p.
In several applications, associating an input x to a single quantity h✓(x) in output is not necessarily ideal. For instance, when p(y|x) is multi-modal, estimating the mean or median might not yield useful predictions for the given task. Moreover, estimators of the form h✓ require access to the full input x to return a prediction, and therefore cannot be used when some features are missing (e.g. due to data corruption). In these contexts, an alternative viable strategy is to directly model the conditional probability. When using PSD models, conditional estimation can be performed in two steps, by first modeling the joint distribution p(y, x) = f(y, x ; A, [Y,X], (⌘, ⌘0)) (e.g. by learning it as suggested Sec. 3) and then use the operations in Sec. 2.1 to condition it with respect to x0 2 X as
p(y|x0) = p(y, x0)
p(x0) = f(y ; B, Y, ⌘) with B =
A (KX,x0,⌘0K > X,x0,⌘ 0)
c2⌘Tr A (KX,x0,⌘0K > X,x0,⌘ 0)KY,Y,⌘
. (19)
In case of data corruption, it is sufficient to first marginalize p(y, x) on the missing variables and then apply (19). Below we discuss a few applications of the conditional estimator.
Conditional Expectation. Conditional mean embeddings [31] are a well-established tool to efficiently compute the conditional expectation Ey⇠p(·|x0) g(y) of a function g : Y ! R. However, although they enjoy good approximation properties [16], they to not guarntee the resulting estimator to take only non-negative values. In contrast, when p is a PSD model (or an approximation), we can apply Prop. 9 to p(·|x0) in (19) and evaluate the conditional expectation of any g for which we know how to compute (or approximate) the integral cg,⌘(z) = R g(y)e ⌘ky zk 2
dy. In particular we have Ey⇠p(·|x) g(y) = Tr((B KY,Y,⌘/2)G) with B as in (19) and G the matrix with entries Gij = cg,2⌘( yi+yj 2 ). Remarkably, differently from conditional mean embeddings estimators, this strategy allows us to compute the conditional expectations also of functions g not in H⌘ .
Structured Prediction. Structured prediction identifies supervised learning problems where the output space Y has complex structures so that it is challenging to find good parametrizations for h✓ : X ! Y [4, 19]. In [7], a strategy was proposed to tackle these settings by first learning an approximation ✓(z, x) ⇡ Ey⇠p(·|x) `(z, y) and then model h✓(x) = argminz2Y ✓(z, x). However, the resulting function ✓(·, x) is not guaranteed to be convex, even when ` is convex. In contrast, by combining the conditional PSD estimator in (19) with the reasoning in Sec. 4.1, we have a strategy that overcomes this issue: when p(·|x) is a PSD model approximating p, we can compute its gradient
Algorithm 1 PSD Hidden Markov Model Input: Transition ⌧̂(x+, x) = f(x+, x ; B, [X+, X], (⌘+, ⌘)), initial p̂(x0) = f(x0 ; A0, X0, ⌘0) and observation !̂(y, x) = f(y, x ; C, [Y,X 0], (⌘obs, ⌘0)) distributions. X̃ as in Prop. 10. ⌘̃ =
⌘(⌘0+⌘+) ⌘+⌘0+⌘+ .
X̃ 0 = ⌘⌘+⌘0+⌘+X ⌦ 1nm + ⌘0+⌘+ ⌘+⌘0+⌘+ 1n ⌦ X̃, and ⌘̃0 = ⌘ 0⌘+ ⌘0+⌘+
For any new observation yt: Ct = C (K>Y,yt,⌘obsKY,yt,⌘obs) // Partial evaluation !̂t(x+) = !̂(y = yt, x+) Bt = (B ⌦At 1) vec(KX,X̃,⌘̃)vec(KX,X̃,⌘̃)
> // Product ̂t(x+, x) = ⌧̂(x+, x)p̂(x|y1:t 1) Dt = (In ⌦ 1>nm)(Bt KX̃0,X̃0, ⌘̃2 )(In ⌦ 1nm) // Marginalization ̂t(x+) = R t(x+, x) dx Et = (Ct ⌦Dt) vec(K
X+,X0, ⌘̃0 2 )vec(K X+,X0, ⌘̃0 2 )>
// Product ⇡̂t(x+) = !̂t(x+)̂t(x+)
At = Et/ct, with ct = c2(⌘0+⌘+)Tr(EtKX̃,X̃, ⌘ 0+⌘+ 2 ) // Normalization p̂(x+|y1:t) = ⇡̂t(x+)R ⇡̂(x+) dx
Return f(xt ; At, X̃, ⌘0 + ⌘+)
Ey⇠p(·|x) rz`(z, y) as mentioned in Sec. 4.1 using Prop. 9. Moreover, if `(·, y) is convex, the term Ey⇠p(·|x) `(·, y) is also convex, and we can use methods such as gradient descent to find h(x) exactly.
Mode Estimation. When the output distribution p(y|x) is multimodal, having access to an explicit form for the conditional density can be useful to estimate its modes. This problem is typically non-convex, yet, when the output y belongs to a small dimensional space (e.g. in classification or scalar-valued regression settings), efficient approximations exist (e.g. bisection).

heading: 4.3 Inference on Hidden Markov Models
text: We consider the problem of performing inference on hidden Markov models (HMM) using PSD models. Let (xt)t2N and (yt)t2N denote two sequences of states and observations respectively. For each t 1, we denote by x0:t = x0, . . . , xt and y1:t = y1, . . . , yt and we assume that p(xt|x0:t 1, y1:t 1) = p(xt|xt 1) = ⌧(xt, xt 1) and p(yt|x0:t, y1:t 1) = p(yt|xt) = !(yt, xt) with ⌧ : X ⇥ X ! R+ and ! : Y ⇥ X ! R+ respectively the transition and observation functions. Our goal is to infer the distribution of possible states xt at time t, given all the observations y1:t and a probability p(x0) on the possible initial states. We focus on this goal for simplicity, but other forms of inferences are possible (e.g. estimating xm+t, namely m steps into the future or the past). We assume that the transition and observation functions can be well approximated by PSD models ⌧̂ and !̂ (e.g. by learning them or known a-priori for the problem’s dynamics). Then, given a PSD model estimate p̂(x0) of the initial state p(x0), we can recursively define the sequence of estimates
p̂(xt|yt:1) = ⌧̂(yt, xt) R !(xt, xt 1)p̂(xt|y1:t 1) dxt 1R
⌧̂(yt, xt)!(xt, xt 1)p̂(xt|y1:t 1) dxtdxt 1 . (20)
Note that when ⌧̂ , !̂, p̂(x0) correspond to the real transition, observation and initial state probability, the formula above yields the exact distribution p(xt|y1:t) over the states at time t (this follows directly by subsequent applications of Bayes’ rule. See also e.g. [5]). If ⌧̂ , !̂, p̂(x0) are PSD models, then each of the p̂(xt|yt:1) is a PSD model recursively defined only in terms of the previous estimate and the operations introduced in Sec. 2.1. In particular, we have the following result. Proposition 10 (PSD Hidden Markov Models (HMM)). Let X0 2 Rn0⇥d, X+, X 2 Rn⇥d, X 0 2 Rm⇥d, Y 2 Rm⇥d0 , A0 2 Sn0+ , A 2 Sn+, B 2 Sm+ and ⌘0, ⌘, ⌘0, ⌘+ 2 Rd++, ⌘obs 2 Rd 0 ++. Let
⌧̂(x+, x) = f(x+, x ; B, [X+, X], (⌘+, ⌘)), !̂(y, x) = f(y, x ; C, [Y,X 0], (⌘obs, ⌘ 0)), (21)
be approximate transition and observation functions. Then, given the initial state probability p̂(x0) = f(x0 ; A0, X0, ⌘0), for any t 1, the estimate p̂ in (20) is a PSD model of the form
p̂(xt|yt:1) = f(xt ; At, X̃, ⌘ 0 + ⌘+), (22)
where X̃ = ⌘ 0
⌘0+⌘+ X 0 ⌦ 1n +
⌘+
⌘0+⌘+ 1m ⌦X+ and At is recursively obtained from At 1 as in Alg. 1.
Remark 4 (Sum-product Algorithm). Eq. (20) is an instance of the so-called sum-product algorithm, a standard inference method for graphical models [5] (of which HMMs are a special case). The application of the sum-product algorithm relies mainly on sum and product rules for probabilities (as is the case for HMMs in (20)). Hence, according to Sec. 2.1, it is highly compatible with PSD models.

heading: 5 Discussion
text: In this work we have shown that PSD models are a strong candidate in practical application related to probabilistic inference. They satisfy both requirements for an ideal probabilistic model: i) they perform exact sum and product rule in terms of efficient matrix operations; ii) we proved that they can concisely approximate a wide range of probabilities.
Future Directions. We identify three main directions for future work: i) when performing inference on large graphical models (see Remark 4) the multiplication of PSD models might lead to an inflation in the number of base points. Building on our compression strategy, we plan to further investigate low-rank approximations to mitigate this issue. ii) An interesting problem is to understand how to efficiently sample from a PSD model. A first answer to this open question was recently given in [15]. iii) The current paper has a purely theoretical and algorithmic focus. In the future, we plan to investigate the empirical behavior of PSD models on the applications introduced in Sec. 4. Related to this, we plan to develop a library for operations with PSD models and make it available to the community.
Acknowledgments. A.R. acknowleges the support of the French government under management of Agence Nationale de la Recherche as part of the “Investissements d’avenir” program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute) and the support of the European Research Council (grant REAL 947908). C.C. acknowledges the support of the Royal Society (grant SPREM RGS\R1\201149) and Amazon.com Inc. (Amazon Research Award – ARA 2020)

Review summaries:
The authors discuss an application of recent work in nonnegative function estimation to the problem of density estimation.  In particular, they argue that under mild assumptions, PSD models can concisely model a large range of interesting densities and that these densities can be (more or less) efficiently learned.  The work concludes with a qualitative comparison against existing density estimation strategies.

This paper is based on the positive semi-definite model from reference [15], and propose a positive semi-definite (PSD) models for probability distribution. The paper then focuses on the family of Gaussian PSD, which can be written as a sum of weighted pairwise multiplicative kernels over sample points. Most importantly, the weights do not need to be positive, and it only requires that the matrix of weights is a symmetric positive semidefinite matrix to ensure that the function is positive everywhere. The authors then show that the PSD model facilitates the sum and product rules; and as the number of sample grows large, one can use projection to compress the model. Further, the authors show the approximation capacity of the PSD model and demonstrate the usefulness of the PSD model on some simple tasks.