Log-concave Density Estimation with Independent Components^††thanks: Submitted July 26, 2025. \fundingElina Robeva was supported by an NSERC Discovery Grant DGECR-2020-00338.

Shape-constrained density estimation has a long history. Classes of shape-constrained densities that have previously been studied include non-increasing [18], convex [19], $k$-monotone [5], $s$-concave [14] and quasi-concave densities [25]. The log-concave class of densities was first studied in dimension 1 by Walther [37], and a concise description of the log-concave MLE in higher dimensions was given by Cule et al. [11]. Various statistical properties of the procedure were established by Schuhmacher and Dümbgen [34], Carpenter et al. [7], Kur et al. [27], and Barber and Samworth [6] amongst others. Log-concavity has also gained popularity in applications [4, 33].

Subclasses of log-concave densities have also been studied in the recent years for the purposes of reducing the computational and statistical complexities of estimation, and also because they are important on their own. In the 1-dimensional setting, Kim et al. [23] study the set of log-concave densities that are exponentials of concave functions defined by at most $k$ linear hyperplanes. In this case, they show that the statistical rate becomes parametric. Xu and Samworth [39] consider log-concave densities whose super-level sets are scalar multiples of a fixed convex body in $\mathbb{R}^{d}$, and establish univariate rates for estimation. Robeva et al. [30] consider maximum likelihood estimation for totally positive and log-concave densities, which correspond to random vectors whose components are highly positively correlated. They prove existence and uniqueness of the MLE and propose an algorithm for computing it, although the statistical rates for this estimator are still unknown. Next, Kubjas et al. [26] study the family of log-concave densities that factorize according to an undirected graphical model, also providing existence and uniqueness results for the MLE (when the graph is given and chordal) together with an algorithm. Sample complexity rates for this model have not been studied either.

Here, we study the subfamily of log-concave densities on $\mathbb{R}^{d}$ obtained by transforming product densities through orthogonal maps. These correspond to random vectors whose coordinates become independent post some orthogonal transformation. The most similar family to ours was considered by Samworth and Yuan [32], who study the class of log-concave densities of random vectors whose coordinates become independent after a linear transformation (rather than just an orthogonal transformation). They show that the projection (in a Kullback-Leibler (KL) sense) of a (not necessarily log-concave) distribution satisfying the above independent components property, to the set of log-concave distributions, still satisfies the independent components property with the same linear transformation. They also provide a method to perform density estimation in this setting, although their algorithm does not guarantee convergence to a global optimum. On the statistical front, they prove consistency; to the best of our knowledge however, finite sample rates have not yet been established in this setting. By simplifying the model here, we are able to provide a very fast algorithm together with sample complexity upper bounds displaying univariate rates.

Towards applicability to real data for tasks such as clustering, we show that our estimator can be used in conjunction with the Expectation-Maximization (EM) algorithm [13] to estimate mixtures of densities from our subfamily – this is a semiparametric generalization of Gaussian mixtures. A similar application was proposed in Cule et al. [11], where they show that mixtures of general log-concave densities can work better than Gaussian mixtures. However, estimating a mixture of general log-concave densities can be very costly, particularly in higher dimensions. Plugging our proposed estimator instead, into the EM algorithm, yields a much faster computation.

Suppose $X$ is a zero-mean random vector in $\mathbb{R}^{d}$, with a (probability) density $p:\mathbb{R}^{d}\to[0,\infty)$. We say that $X$ (or $p$) satisfies the orthogonal independent components property if there exists an orthogonal matrix $\mathbf{W}\in\mathbb{R}^{d\times d}$ such that the components of $Z=\mathbf{W}X$ are independent. In this case, the joint density of $Z$ can be expressed as the product

where each $f_{i}:\mathbb{R}\to[0,\infty)$ is a univariate density. We focus mainly on the case where $p$ is log-concave, which is equivalent in this setting to requiring log-concavity of $f_{1},\dots,f_{d}$. Our goal is to efficiently estimate $p$ from independent and identically distributed (iid) samples, given that it satisfies the orthogonal independent components property with some unknown $\mathbf{W}$.

By writing $X=\mathbf{W}^{-1}Z$, we can interpret this setting as an independent component analysis (ICA) model, where $Z$ is the vector of independent sources and $\mathbf{W}$ is the so-called unmixing matrix. Although mixing matrices in ICA do not need to be orthogonal to begin with, a pre-whitening procedure [10] can transform the problem to such an orthogonal form.

Instead of estimating the “$d$-dimensional” density $p$ directly, we would like to leverage the underlying structure provided by (1), which allows us to simplify $p$ (using a change-of-variables formula) as

Here, $w_{1},\dots,w_{d}\in\mathbb{S}^{d-1}$ are the rows of $\mathbf{W}$, and we refer to them as independent directions of $p$. Note that the model $Z=\mathbf{W}X$ and the factorization in (2) are invariant under permuting the rows of $\mathbf{W}$ and the corresponding components of $Z$, as well as under flipping their respective signs.

If $\mathbf{W}$ is known beforehand – if it is provided by an oracle – the estimation problem is simplified, and one can use the following procedure:

1. 1.

   For each $i=1,\dots,d$, estimate the univariate density $f_{i}$ (say via log-concave maximum likelihood estimation [33]) using the samples $\{Z_{i}^{(\ell)}=w_{i}\cdot X^{(\ell)}:\ell\in[N]\}$, where $X^{(1)},\dots,X^{(N)}\overset{\mathrm{iid}}{\sim}p$. Denote the estimated density by $\hat{f_{i}}$.

2. 2.

   Plug-in the estimated densities into (2) to obtain the oracle-informed estimate of $p$:

   (3)

   $\displaystyle\hat{p}_{\text{oracle}}(x)=\prod_{i=1}^{d}\hat{f_{i}}(w_{i}\cdot x).$

In reality, the unmixing matrix $\mathbf{W}$ is unknown, and one needs to estimate it. If the covariance matrix of $X$ has distinct, well-separated eigenvalues, then this can be done by Principal Component Analysis (PCA), which involves diagonalizing the sample covariance matrix of $X$. Otherwise, certain ICA algorithms may be used to estimate $\mathbf{W}$ provided some non-Gaussianity requirements are satisfied [3, 17]. We refrain from using the same samples for the unmixing matrix estimation and the density estimation steps. Instead, we split the samples into two subsets, $\{X^{(1)},\dots,X^{(N)}\}$ and $\{Y^{(1)},\dots,Y^{(M)}\}$, and use different subsets for the two purposes in order to prevent potentially problematic cross-interactions.

Denote the (PCA or ICA) estimate of $\mathbf{W}$ by $\hat{\mathbf{W}}$. We always require $\hat{\mathbf{W}}$ to be an orthogonal matrix. Given enough samples, the rows $\hat{w}_{1},\dots,\hat{w}_{d}\in\mathbb{S}^{d-1}$ of $\hat{\mathbf{W}}$ can be shown to be close to those of $\mathbf{W}$ up to permutations and sign-flips. More precisely, there exists a permutation $\sigma:[d]\to[d]$ and signs $\theta_{i}\in\{-1,+1\}$ such that for some threshold $\epsilon>0$

As the factorization in (2) is invariant under such permutations and sign-flips, one can simply relabel $w_{i}\leftarrow\theta_{i}w_{\sigma(i)}$ and $Z_{i}\leftarrow\theta_{i}Z_{\sigma(i)}$ to proceed with the analysis.

From $\hat{\mathbf{W}}$, we can follow a procedure analogous to that used to obtain the oracle-informed estimate of $p$. Namely, we use the estimated independent directions $\hat{w}_{i}\in\mathbb{R}^{d}$ and the projected samples onto each such direction to produce univariate estimates $\hat{p}_{\hat{w}_{i}}$; then we take the product of these estimates over $i$ (see Definition 2.4).

The main contribution of this work is an algorithm for density estimation in the setting of Section 1.2, for which we prove finite sample guarantees and demonstrate numerical performance in terms of estimation errors and runtimes. The proposed algorithm is recorded as Algorithm 1.

The theoretical analysis in Section 3 yields sample complexity bounds for our algorithm under various assumptions, which are listed in Table 1 and Table 2. Briefly, if the true density $p$ is log-concave and satisfies some moment assumptions in addition to the orthogonal independent components property, then it can be estimated efficiently. In particular, $\tilde{\mathcal{O}}_{d}(\epsilon^{-4})$ samples (up to constants and log factors) are sufficient for our proposed estimator $\hat{p}$ to be within $\epsilon$ of $p$ in squared Hellinger distance with probability 0.99. If the density $p$ has some additional smoothness, the rate can be improved up to $\epsilon^{-2}$.

Since we referred to smoothness, a quick comparison with kernel density estimators (KDE) is in order. Recall that in dimension $d$, a KDE typically requires the true density $p$ to have “smoothness of order $d$” in order to achieve a dimension-independent statistical rate. More precisely, for $p$ in the Sobolev space $\mathcal{W}^{s,1}(\mathbb{R}^{d})$ with $s\geq 1$, achieving an expected total variation error of $\epsilon$ requires $n\gtrsim_{d}\epsilon^{-\frac{d+2s}{s}}$ samples [21] so that one needs $s\gtrsim d$ to counter the curse of dimensionality. On the contrary, our smoothness assumptions (see Smoothness assumption S1 and Smoothness assumption S2) are “order one” in that we only require control on the first-order (and if available, the second-order) derivatives.

Computational experiments are presented in Section 4, where we demonstrate improved estimation errors (in squared Hellinger distance) and improved runtimes as compared to usual log-concave MLE. A key consequence of our fast algorithm is easy scalability to higher dimensions – we are able to estimate densities in $d=30$ within no more than a few seconds.

We consider the oracle-informed estimator first. The analysis involved is also applicable to the proposed estimator.

Lemma 3.1 says that the problem of bounding $h_{d}^{2}(\hat{p}_{\text{oracle}},p)$ boils down to controlling the univariate estimation errors $h_{1}^{2}(\hat{p}_{w_{i}},p_{w_{i}})$. For log-concave densities, Kim and Samworth [24] and Carpenter et al. [7] have already established results of this kind. We stay in a general setting for now, by deriving the sample complexity of the oracle estimator from the sample complexity of any univariate density estimation procedure that it invokes.

Now, we can immediately express the sample complexity of the oracle estimator in terms of $N^{*}$ as a consequence of Lemma 3.1 and a union bound.

Having concluded the analysis of the oracle-informed estimator, we can now turn to the general setting. Here, we only have an estimate of $\mathbf{W}$, and hence, the estimation error needs to be taken into account. This contributes additional terms to the risk of the proposed estimator, as the theorem below shows. The proof can be found in Section A.1.

Here, $\hat{\mathbf{R}}$ is an orthogonal matrix, because $\mathbf{W}$ and $\hat{\mathbf{W}}$ are. The term $\max_{i\in[d]}h_{1}(\hat{p}_{\hat{w}_{i}},p_{\hat{w}_{i}})$ in the bound of Theorem 3.5 warrants an analysis similar to the oracle-informed estimator. Since $\hat{p}_{\hat{w}_{i}}$ is a univariate density estimate of $p_{\hat{w}_{i}}$ conditional on $\hat{w}_{i}$, bounds on univariate density estimation control this term, when applied in a conditional sense.

The second and third terms in the bound of Theorem 3.5 call for a notion of stability; $h_{d}(\hat{\mathbf{R}}_{\sharp}p,p)$ measures how far the ‘rotated’ density $\hat{\mathbf{R}}_{\sharp}p$ gets from the original density $p$. If Hellinger distance is stable with respect to small rotations of the density, then $h_{d}(\hat{\mathbf{R}}_{\sharp}p,p)$ can be bounded effectively in terms of the operator norm $\|\mathbf{I}-\hat{\mathbf{R}}\|_{\mathrm{op}}$, where $\mathbf{I}$ is the $d\times d$ identity matrix. This operator norm is then well-controlled by the unmixing matrix estimation error as

The next lemma demonstrates such stability for zero-mean log-concave densities.

The proof involves the log-concave projection of Dümbgen et al. [15], particularly making use of its local Hölder continuity [6]. It is delayed to Section A.2. Note that Lemma 3.6 is quite general, as it only uses log-concavity of $p$ (with no additional smoothness requirements). However, Lemma 3.6 inherits a constant $K_{d}$ that typically depends at least exponentially in $d$. Also, the exponent $1/2$ of $\|\mathbf{I}-\mathbf{A}\|_{\mathrm{op}}$ contributes a $1/\epsilon^{4}$ dependence to our sample complexity bounds (see the first two rows of Table 2). Later in Section 3.4, we show that smoothness assumptions on $p$ can give stronger bounds.

The last bit to consider is the unmixing matrix estimation error. In the same spirit as the discussion on univariate density estimation above, we consider a generic estimator for now, with sample complexity defined as below.

Now, we state the main result of this section, which expresses the sample complexity of the proposed estimator (Definition 2.4) in terms of the sample complexities in Definition 3.3 and Definition 3.7. The proof is presented in Section A.3.

To obtain concrete sample complexities, one needs to choose specific estimation procedures. In Section 3.2, we consider univariate log-concave MLE and its boosted variant, and provide corresponding bounds on $N^{*}$. In Section 3.3, we present standard bounds on $M^{*}$ for PCA- and ICA-based estimators of $\mathbf{W}$.

Recall that for each $i=1,\dots,d$, we need to estimate the densities $p_{\hat{w}_{i}}$ from samples $\{\hat{Z}_{i}^{(j)}=\hat{w}_{i}\cdot X^{(j)}:j\in[M]\}$. When $p$ is log-concave, an efficient as well as convenient option is log-concave MLE [33], which is in fact minimax optimal.

In terms of Definition 3.3, we can restate this as $N^{*}(\epsilon,\gamma;f)\leq C_{1}\epsilon^{-5/4}\log^{5/4}\left(C_{2}/\gamma\right)\vee C_{3}/\gamma$. Notice that the right-hand side does not depend on $f$, and as a result, the same upper bound also holds for the term $\sup_{s\in\mathbb{S}^{d-1}}N^{*}(\epsilon,\gamma;p_{s})$ from Theorem 3.8

The statement above is a “high-probability” variant of Theorem 5 of Kim and Samworth [24], and we provide a proof in Section A.4. Note that while the $\epsilon^{-5/4}$ dependence is optimal, the $1/\gamma$ dependence in the sample complexity could be improved. A standard way to do this is via boosting. The following definition uses the boosting technique from Algorithm A of Lovász and Vempala [28].

As the next proposition show, the sample complexity of the boosted log-concave MLE depends only logarithmically on $1/\gamma$, which is much milder. See Section A.5 for the proof.

In particular, by multiplying the bounds for $n$ and $q$, we get that the sample complexity
$N^{*}(\epsilon,\gamma;f)\leq C_{1}C_{2}\epsilon^{-5/4}\log(1/\gamma)$. As before, the right-hand side does not depend on $f$.

Here, we discuss standard techniques for estimating $\mathbf{W}$, and show closeness of $\hat{w}_{i}$ and $w_{i}$. Recall that the methods used here for estimating $\mathbf{W}$ – PCA or ICA – require different assumptions on the moments of $X$ (or $Z=\mathbf{W}X$) in order to yield bounds on the estimation errors.

Consider regular PCA first. The PCA estimate of $\mathbf{W}$ was obtained by computing the eigenvectors of the sample covariance matrix $\hat{\mathbf{\Sigma}}$. Bounding the sample complexity of this estimation procedure is standard, and can be broken down into two simple steps: (i) bounding the distance between the empirical covariance matrix $\hat{\mathbf{\Sigma}}$ and the true covariance matrix $\mathbf{\Sigma}$, and (ii) bounding the distances between their corresponding eigenvectors. The result is summarized in the following proposition, the proof of which is delayed to A.6.

The right-hand side of (9) is an upper bound on sample complexity $M^{*}(\epsilon,\gamma;p)$ of the PCA-based estimator, for $0<\gamma\leq 1/e$. The term $\tilde{M}(d,\gamma)$ is needed in order to allow the result above to hold for all positive $\epsilon$, instead of just small values of $\epsilon$.

As discussed earlier, assumptions on strict eigenvalue separation can be stringent, and possibly unrealistic in practical settings. In such cases, $\mathbf{W}$ should be estimated by ICA instead, for which we have a result similar to Proposition 3.12.

Recall the definitions of $S$ and $\mathbf{A}$ from Section 2, and express $X=\mathbf{A}S$. The algorithm of Auddy and Yuan [3] then estimates the columns of $\mathbf{A}$, which coincide with $w_{1},\dots,w_{d}\in\mathbb{S}^{d-1}$ up to scalings and signs.

The proposition below is a direct consequence of Theorem 3.4 of Auddy and Yuan [3], and provides an upper bound on $M^{*}(\epsilon,\gamma;p)$ for $3d^{-3}\leq\gamma<1$. The proof is delayed to Section A.7.

The result above relies on non-zero excess kurtosis of the components of $Z$, which may be interpreted as a measure of non-Gaussianity. However, there exist distributions which differ from Gaussians only in higher cumulants. Techniques such as Fourier PCA [17] allow for estimating $\mathbf{W}$ in such cases. As discussed earlier, our theory is agnostic to specific choices of ICA algorithms.

While the statement of Theorem 3.8 is quite general, the inequality in (8b) contains a constant $K^{\prime}_{d}>0$ that typically depends at least exponentially in $d$. This is owed to the stability result Lemma 3.6, which uses only the log-concavity of $p$. Stronger stability inequalities are possible when $p$ is smooth.

In what follows, we consider two possible assumptions on the smoothness of $p$. The choice of assumptions then allows us to quantify stability in either Kullback-Leibler (KL) divergence or total variation distance first, before translating this to stability in Hellinger distance. Let $p$ be a probability density on $\mathbb{R}^{d}$ with mean zero.