Unnormalized Density Estimation with Root Sobolev Norm Regularization

We thank the reviewer for the feedback.

Sobolev space is closely related to the RKHS with Matern kernel; see [1] below. Some supplementary comments on that relationship would be helpful for researchers interested in theoretical analysis of the learning algorithms.

Thank you for including the reference [1]. It contains a clear relevant results on the Matern kernel in context of Sobolev space and we will include it in the paper.
Please see the discussion in the global "General Reply" comment on the relation between the SDO kernel and Matern kernels.

The estimator is similar to the one proposed in Ferraccioli's JRSS paper. Please clarify the main difference between them.

We discuss the relation with the paper

Nonparametric density estimation over complicated domains, Ferraccioli et al, JRSS 2020

In the above paper, the authors also use a reparametrisation to ensure non-negativity of the density, and regularise using a version of a Sobolev space. However, as we discuss in our introduction and literature review sections, that general principle is known since the 1980s for one dimension. The challenge is how one can implement the principle in $d>1$.

The above mentioned paper is solely restricted to two dimensions, d=2. (both for applications and theory). In particular, to obtain their solutions, they discretize (triangulate) the domain of interest in $R^2$. This is clearly infeasible in any higher dimension.

In contrast, the approach we develop, including the natural gradient optimisation and the use of Fisher divergence for parameter tuning, allow us to obtain good anomaly detection performance on datasets with $d>100$. Moreover, we consider the fact that such methods are useful in high dimensional situations is a significant contribution by itself, since such methods were not considered in recent literature in high dimensions. We will emphasize this more in the discussion in the paper. Would the reviewer agree that this is indeed a contribution of interest?

Section 6.3 reports the ratio of negative values for the estimated function f. I'm not sure why the ratio of negative values is important to assess the computational properties.

This point is important and we now clarify it. As we discuss in the paragraph following Eq. (5), the main optimisation objective (5) is not convex as a function of $f$. Moreover, if one nevertheless tries to optimise (5) in $f$ (i.e., in $\alpha$, using $f_{\alpha}$) directly, one would obtain very poor results. The gradient descent will get stuck in a local optimum with poor values. We can provide an experimental demonstration of this if that would be of interest.

However, we observe that (5) is convex on the non-negative cone, i.e., $f$ s.t. $f>0$ everywhere. (It is also sufficient to require only $f(x_i)>0$ on data points $x_i$).
Thus, if one restricts to such cone, one can obtain the optimal solution on that cone. We observe that such solutions are empirically much better. Next, we note that solving a constrained optimisation problem, requiring the solutions to be restricted to the cone is computationally difficult. Moreover, since the cone is not closed (topologically), this is also not quite clear theoretically, since projection methods can not be applied.

To resolve this, we observe that natural gradient preserves the cone. This means that one can in fact do unconstrained optimisation. This is true for non-negative kernels, such as the Matern kernels (discussion following Equation (5) and Proposition 10), and is approximately true for the SDO kernel.

In the experiment in Section 6.3 we show that empirically, for the SDO kernel, natural gradient indeed tends to preserve positivity, while the regular gradient does not. This demonstrates that natural gradient is indeed better for optimisation purposes. Note that, theoretically, this is one of the very few settings were we know why the natural gradient is better (i.e. due to positivity).

Moreover, we can also directly show that the results in Fig 2, our main empirical results, would be worse if regular gradient descent is used instead of natural.

Does this answer make sense? We would be glad to provide additional clarifications.

Continued in the comment below.

Dear Reviewer CL4E, Thank you again for your feedback. In our earlier comments we have answered the questions posed in the review and addressed the concerns raised in the Weaknesses section.

Have all the concerns been addressed, or are there any issues still remaining?

Thank you for the feedback!

We agree with the point raised in this review that the specific choice of the SDO kernel requires more motivation. Please see a detailed discussion of our motivation, and of relations of SDO with other kernels, in the General Reply comment common to all reviews. In addition, in a few days we will include experiments that evaluate our algorithms with the Laplacian and Gaussian kernels, which correspond to a representative range of Sobolev related kernels and have analytic expressions. As we mention in the above discussion, all contributions of the paper besides the SDO kernel are relevant for these kernels too.

We comment on an additional point on equivalence of Sobolev spaces, which was mentioned in the review. It is true for instance that the spaces corresponding to the Laplacian and SDO kernels are equivalent for all values of the Laplacian parameter $\sigma$ and SDO parameter $a$ (see the discussion above). However, this only means that the spaces contain the same set of functions. This in turn implies that the norms are equivalent, but the equivalence constants diverge with the dimension. Since for regularization we are interested in an actual value of the norm, the resulting algorithms will be different between such spaces. Consider for instance that $\|\cdot \|_1$ and $\|\cdot\|_2$ norms on $R^d$ are equivalent norms, but with drastically different regularisation properties. Thus just norm equivalence is too coarse a notion to distinguish (or unify) between different regularisation algorithms.

Do these comments address the issue with the SDO choice?

(ii) the paper mentions on the top section of p. 2 that regularizing over $\|f\|_{L_2}^2$ is a desirable property: what does it mean?

Our intention here was to say by using the Sobolev type kernel, we automatically deal with functions that satisfy $\inf f^2(x)dx < \infty$. Such functions are densities up to a constant. Non integrable functions on the other hand can not represent densities (for the purposes of sampling from them, for instance). This is not a trivial condition, as there is no reason why an "arbitrary" function should be integrable.

The paper's main method is marked as "INER" in Figures 1 and 2, but the text makes no mention of what it stands for: is this a mistake?

Yes, it should be RSR. Will be fixed.

The paper mentions the consistency of the estimator in the main manuscript, but delays all discussions about consistency to the supplementary materials - I believe that it would be helpful to bring up the main theorems for consistency in the main text.

We will definitely consider this. The challenge here is that there seems to be no other part of the paper which is non crucial enough to be moved to the supplementary instead.

Thank you for the detailed feedback. We were glad to learn that the empirical results were found to be quite impressive.

In the following we discuss all points raised in the review. In particular, we discuss the technical issue of density positivity, the relation to the work [MFBR20] (thanks for bringing it to our attention!), and the purpose of Section 5. Please let us know whether some issues remain.

In addition, we believe the reviewer may find the material in the global "General Reply" comment above to be of interest. There, we attempt to clarify the relation between various Sobolev space related kernels.

First, we address the "positivity on the data points" question in the consistency result. This was raised in point 4 of the review and in the related Question 2.

positivity on the data points:

The reason we only look for solutions in the cone (104) with strictly positive $f(x_i)>0$ is that the gradient of the loss is not defined otherwise (Eq. (6) or (106)). Note that the functions in the cone satisfy $f(x_i)>0$ simply because we defined them to satisfy this, thus, there is no contradiction here.
In addition, we do want to evaluate the gradient of the true density $v$ at the sample points. We believe this was the source of the issue in the review, when $n\rightarrow \infty$. However, this also does not constitute a problem. Indeed, note that we do not require a lower bound on $f(x_i)$, we just need to know the gradient is well defined, i.e., $f(x_i)\neq 0$. To this end, for any continuous density $v^2$, we have the following:

Let $\{x_i\}$ be an infinite sequence of iid samples from $v^2$. Then $v(x_i)>0$ for all $i$ with probability 1. (since choosing $x$ with $f(x) =0$ is a 0 probability event).

As an example, let $v^2$ be either a Gaussain, or a uniform density on a sphere. In both cases samples will have $v^2(x_i)>0$. The infimum of such values will certainly be 0 for the Gaussian but there is no contradiction in that. Such points, approaching 0, will not be optimal, and again, we only need the gradient to be well defined.

Does this resolve the issue? We would be glad to discuss it further.

Next we address the other major question, i.e., the relation to work [MFBR20].

relation to [MFBR20]:

At a very high level, the approach of [MFBR20] is indeed similar to ours, since they consider quadratic functions and regularisation with a version of an RKHS norm. The differences are both in scope and in details of the construction that enable computation. From a computational viewpoint, they optimise in the space on non-negative matrices. This makes the problem convex without any need for positivity, but it requires $N^2$ parameters to encode such a matrix (their Theorem 1), and $O(N^3)$ computational cost per step (their top of page 6). Further, the optimisation is a constrained one, which is also significantly more difficult. It is unlikely their methods can work on something like ADBench benchmark, and indeed, they only test it on a 1 dim Gaussian with N=50.

Second, while [MFBR20] is concerned with other applications as well, we investigate the density estimation problem much more closely. We prove consistency, which is completely orthogonal to all their results, we prove the method is different from KDE, and we do provide evaluation on a state of the art benchmark of datasets.

Finally, they do treat the normalisation (their Prop 4, and paragraph below it). However, this is simply by a straight forward opening of the brackets. That is (in our case), if $f(x) = \sum_{i} \alpha_i k(x_i,x)$ then $\int f^2 (x)dx = \sum_{i,j} \alpha_i \alpha_j k(x_i,x)k(x_j,x)$.

Thus, if one knows $\int k(x_i,x)k(x_j,x) dx$ for every $i,j$ one can compute theoretically the normalisation constant. The issues with this are as follows: (a) This would involve summing $N^2$ floating point numbers, which would be inherently unstable for large N. (b) Computing $\int k(x_i,x)k(x_j,x) dx$ analytically is easy only for the Gaussian kernel. Even for exponential (Laplacian) kernel in $d>1$, we do not see any way to evaluate this analytically. Of course, for the SDO kernel, this also does not seem to be computable. Due to this reason, we have not pursued normalisation.

It is worth noting however, that while the work of [MFBR20] advocates for using the computation above as a constraint, if one uses the Sobolev norms as we do, the integral is guaranteed to be finite, and it is thus sufficient to normalise only once, when the optimisation is finished.

Continued in the comment below.

Thanks to the authors for their responses to my questions. The discussion of [MFBR20] was helpful to me; in particular, it seems that the optimization iterations of the proposed approach are $O(n^2)$, in contrast to the $O(n^3)$ iterations of [MFBR20], allowing the proposed method to scale to larger datasets.

However, I didn't understand the authors' respose to my question about the "positivity on the data points" condition. Please let me know if I'm missing something here.

In addition, we do want to evaluate the gradient of the true density

The parameters of the true density are fixed; what does it mean to evaluate "the gradient of the true density"? Doesn't (natural or standard) gradient descent actually require taking the gradient of $\alpha$ at the current iterate, not at the true value?

To be explicit, let's consider the initial value $\alpha_0$. I believe computing the natural gradient (Eq. (6)) requires that the estimated density is positive on the data points, i.e., $f_0 := \min_{i = 1,...,N} \sum_{j = 1}^N \alpha_{0,j} k_{x_j}(x_i) > 0$. I agree that this is satisfied if $f_0$ lies in the positive cone. However, the paper claims that enforcing this cone constraint is computationally intractable and therefore suggests an alternative approach in practice, namely initializing each $\alpha_{0,j} \geq 0$. As noted in the paper, for a potentially negative kernel, such as the SDO kernel, $\sum_{i = 1}^N \alpha_{0,i} k_{x_i}(x)$ may still take negative values at some $x$. So it seems necessary to me to justify that, for $\alpha_{0,j} \geq 0$, $\min_{i = 1,...,N} \sum_{j = 1}^N \alpha_{0,j} k_{x_j}(x_i) > 0$.

This seems unlikely to be true in general. In fact, it's easy to construct a counterexample when $N = 2$ (let $\alpha_{0,1} \gg \alpha_{0,2}$, and let $x_2$ be such that $k_{x_1}(x_2) < 0$). I think this creates a hole in the consistency argument. I was wondering, however, if there is a way to choose $\alpha$ such that this is true with probability approaching 1 for large $N$? If so, then this might not be a problem in practice, and showing this would solidify the consistency argument.

EDIT: Another possibility that just occured to me is that the consistency guarantee is intended to be for the "ideal" algorithm, which exactly enforces the cone constraint. If this is the case, then I guess the guarantee stands, but there is a bit of a gap between the theory and the practice. Is this the intent?

Regarding my point about the paper's motivation

I'm not convinced by this motivation, for a two reasons...

I fully agree with the authors' claim that

all methods have a certain bias, and there [is] no "best" method.

In fact, I did not mean to question the motivation underlying the search for a new method (sorry it came across in this way!). Rather, I am commenting on the specific sentences

While there is recent work for low dimensional (one or two dimensional) data... there still are very few non-parametric methods applicable in higher dimensions.

These sentences don't seem accurate or convincing to me, and I am suggesting that the authors re-write this paragraph to identify more specific strengths or motivations of their proposed method.

Thus, the Laplacian requires $d/2$ derivatives, a minimal requirement for the Sobolev RKHS, and the Gaussain requires infinite smoothness. The Matern kernels have (quite recently) been shown to interpolate between these extremes, requiring $\nu > d/2$ derivatives. The Laplacian and Gaussian are special instances of this class; see, for instance, Example 5.7 in [6].

The relation between these kernels is now clear. The Laplacian, Gaussian and Matern kernels enjoy closed analytic forms, but this comes at the expense of a quite complex weighting of different orders of derivatives, and the bandwidth parameters $\sigma$ affecting all of the orders, but with an exponential scaling. Imagine trying to explain this for instance to a medical doctor, in a context of interpretability discussion of your model. On the other hand, SDO as we use here is similar to the Laplace in requiring only a $d/2$-th smoothness, but we measure only the size of that derivative.

Finally, we would like to mention that for Matern kernels, the positivity preserving properties of the natural gradient apply exactly (see the discussion following Eq. (5)) rather than just approximately, since these are non-negative kernels.

We hope this clarifies the contribution. The empirical evaluation of RSR with Laplace and Gaussian kernels will be included in the few coming days. We would be glad to answer any further questions.

[1] Cox, Multivariate smoothing spline functions, SIAM Journal on Numerical Analysis, 1984

[2] Delecroix, Simioni Thomas-Agnan, Functional estimation under shape constraints, Journal of Nonparametric Statistics, 1996

[3] Novak, Ullrich, Wo{'z}niakowski,Zhang, Reproducing kernels of Sobolev spaces on R^d and applications to embedding constants and tractability, Analysis and Applications, 2018

[4] Berlinet, Thomas-Agnan, Reproducing kernel Hilbert spaces in probability and statistics,Springer Science, 2011

[5] Thomas-Agnan, Computing a family of reproducing kernels for statistical applications, Numerical Algorithms, 1996

[6] Fasshauer, Ye, Reproducing kernels of generalized Sobolev spaces via a Green function approach with distributional operators, Numerische Mathematik, 2011

[7] Rakhlin, Zhai, Consistency of interpolation with Laplace kernels is a high-dimensional phenomenon, Conference on Learning Theory, 2019

[8] Rasmussen, Williams, Gaussian processes for machine learning, 2006, Springer