# Numerical Reconstruction of the Covariance Matrix of a Spherically Truncated Multinormal Distribution

Authors: Filippo Palombi; Simona Toti; Romina Filippini
Journal: Journal of Probability and Statistics
 (2017)

---

## Abstract

We relate the matrixSB of the second moments of a spherically truncated normal multivariate to its full covariance matrix Σ and present an algorithm to invert the relation and reconstruct Σ from SB. While the eigenvectors of Σ are left invariant by the truncation, its eigenvalues are nonuniformly damped. We show that the eigenvalues of Σ can be reconstructed from their truncated counterparts via a fixed point iteration, whose convergence we prove analytically. The procedure requires the computation of multidimensional Gaussian integrals over an Euclidean ball, for which we extend a numerical technique, originally proposed by Ruben in 1962, based on a series expansion in chi-square distributions. In order to study the feasibility of our approach, we examine the convergence rate of some iterative schemes on suitably chosen ensembles of Wishart matrices. We finally discuss the practical difficulties arising in sample space and outline a regularization of the problem based on perturbation theory.


---

## Body

## 1. Introduction

It has been more than forty years since Tallis [1] worked out the moment-generating function of a normal multivariate X≡{Xk}k=1v~Nv(0,Σ), subject to the conditional event(1)X∈Evρ;Σ,Evρ;Σ≡x∈Rv:xTΣ-1x≤ρ.The perfect match between the symmetries of the ellipsoid Ev(ρ;Σ) and those of Nv(0,Σ) allows for an exact analytic result, from which the complete set of multivariate truncated moments can be extracted upon differentiation. Consider, for instance, the matrix SE(ρ;Σ) of the second truncated moments, expressing the covariances among the components of X within Ev(ρ;Σ). From Tallis’ paper it turns out that(2)SEρ;Σ=cTρΣ,cTρ≡Fv+2ρFvρ,with Fv denoting the cumulative distribution function of a χ2-variable with v degrees of freedom. Inverting (2)—so as to express Σ as a function of SE—is trivial, since cT(ρ) is a scalar damping factor independent of Σ. In this paper, we shall refer to such inverse relation as thereconstruction of Σ from SE. Unfortunately, life is not always so easy. In general, the effects produced on the expectation of functions of X by cutting off the probability density outside a generic domain D⊂Rv can be hardly calculated in closed form, especially if the boundary of D is shaped in a way that breaks the ellipsoidal symmetry of Nv(0,Σ). Thus, for instance, unlike (2), the matrix of the second truncated moments is expected in general to display a nonlinear/nontrivial dependence upon Σ.In the present paper, we consider the case whereD is a v-dimensional Euclidean ball with center in the origin and square radius ρ. Specifically, we discuss the reconstruction of Σ from the matrix SB of the spherically truncated second moments. To this aim, we need to mimic Tallis’ calculation, with (1) being replaced by the conditional event(3)X∈Bvρ,Bvρ≡x∈Rv:xTx≤ρ.This is precisely an example of the situation described in the previous paragraph: although Bv(ρ) has a higher degree of symmetry than Ev(ρ;Σ), still there is no possibility of obtaining a closed-form relation between Σ and SB, since Bv(ρ) breaks the ellipsoidal symmetry of Nv(0,Σ): technically speaking, in this case we cannot perform any change of variable under the defining integral of the moment-generating function, which may reduce the dimensionality of the problem, as in Tallis’ paper.In spite of that, the residual symmetries characterizing the truncated distribution help simplify the problem in the following respects: (i) the reflection invariance of the whole setup still yieldsE[Xk∣X∈Bv(ρ)]=0∀k and (ii) the rotational invariance of Bv(ρ) preserves the possibility of defining the principal components of the distribution just like in the unconstrained case. In particular, the latter property means that SB and Σ share the same orthonormal eigenvectors. In fact, the reconstruction of Σ from SB amounts to solving a system of nonlinear integral equations, having the eigenvalues λ≡λkk=1v of Σ as unknown variables and the eigenvalues μ≡μkk=1v of SB as input parameters. In a lack of analytic techniques to evaluate exactly the integrals involved, we resort to a numerical algorithm, of which we investigate feasibility, performance, and optimization.The paper is organized as follows. In Section2, we describe a few examples illustrating the occurrence of spherical truncations in practical situations. In Section 3, we show that the aforementioned integral equations have the analytic structure of a fixed point vector equation; that is to say, λ=T(λ). This suggests achieving the reconstruction of λ numerically via suitably chosen iterative schemes. In Section 4, we prove the convergence of the simplest of them by inductive arguments, the validity of which relies upon the monotonicity properties of ratios of Gaussian integrals over Bv(ρ). In Section 5, we review some numerical techniques for the computation of Gaussian integrals over Bv(ρ) with controlled systematic error. These are based on and extend a classic work by Ruben [2] on the distribution of quadratic forms of normal variables. For the sake of readability, we defer proofs of statements made in this section to the appendix. In Section 6, we report on our numerical experiences: since the simplest iterative scheme, namely, the Gauss–Jacobi iteration, is too slow for practical purposes, we investigate the performance of its improved version based on overrelaxation; as expected, we find that the latter has a higher convergence rate; yet it still slows down polynomially in 1/ρ as ρ→0 and exponentially in v as v→∞; in order to reduce the slowing down, we propose an acceleration technique, which boosts the higher components of the eigenvalue spectrum. A series of Monte Carlo simulations enables us to quantify the speedup. In Section 7 we discuss the problems arising when μ is affected by statistical uncertainty and propose a regularization technique based on perturbation theory. To conclude, we summarize our findings in Section 8.


## 2. Motivating Examples

Spherical truncations of multinormal distributions may characterize different kinds of experimental devices and may occur in various problems of statistical and convex analysis. In this section, we discuss two motivating examples.


### 2.1. A Two-DimensionalGedanken Experiment in Classical Particle Physics

Consider the following ideal situation. An accelerator physicist prepares an elliptical beam of classical particles with Gaussian transversal profile. The experimenter knowsa priori the spatial distribution of the beam, that is, the covariance matrix Σ of the two-dimensional coordinates of the particles on a plane orthogonal to their flight direction. We can assume with no loss of generality that the transversal coordinate system has origin at the maximum of the beam intensity and axes along the principal components of the beam; thus it holds Σ=diag⁡(λ1,λ2). The beam travels straightforward until it enters a linear coaxial pipeline with circular profile, schematically depicted in Figure 1, where the beam is longitudinally accelerated. While the outer part of the beam is stopped by an absorbing wall, the inner part propagates within the pipeline. At the end of the beam flight the physicist wants to know if the transversal distribution of the particles is changed, due to unknown disturbance factors arisen within the pipeline. Accordingly, he measures again the spatial covariance matrix of the beam. Unfortunately, the absorbing wall has cut off the Gaussian tail, thus damping the covariance matrix and making it no more comparable to the original one. To perform such a comparison in the general case λ1≠λ2, the covariance matrix of the truncated circular beam has to go through the reconstruction procedure described in next sections.Figure 1
A classical particle beam with elliptical transversal profile is cut off upon entering a circular coaxial pipeline.


### 2.2. A Multivariate Example: Connections to Compositional Data Analysis

Compositional Data Analysis (CoDA) has been the subject of a number of papers, pioneered by Aitchison [3] over the past forty years. As a methodology of statistical investigation, it finds application in all cases where the main object of interest is a multivariate with strictly positive continuous components to be regarded as portions of a total amount κ (the normalization κ=1 is conventionally adopted in the mathematical literature). In other words, compositional variates belong to the κ-simplex:(4)Sv=z∈R+v:z1=κ,v≥2,with z1=∑k=1vzk the taxi-cab norm of z, while compositions with different norms can be always projected onto Sv via the closure operator C·x≡{κx1/x1,…,κxv/x1}. There are countless types of compositional data, whose analysis raises problems of interest for statistics [4], for example, geochemical data, balance sheet data, and election data. The simplex constraint induces a kind of dependency among the parts of a composition that goes beyond the standard concept of covariance. This invalidates many ordinary techniques of statistical analysis.In order to measure distances onSv, Aitchison introduced a positive symmetric function dA:Sv×Sv→R+, explicitly defined by(5)dAx,y=12v∑i,k=1vlog⁡xixk-log⁡yiyk2.The Aitchison distance is a key tool in CoDA. It is scale invariant in both its first and second arguments; that is, it is left invariant by redefinitions z→{αz1,…,αzv} with α∈R+. Accordingly, its support can be extended to R+v×R+v by imposing(6)dAx,y≡dAC·x,C·y,x,y∈R+v.It was proved in [5] that dA is a norm-induced metric on Sv, provided the latter is given an appropriate normed vector space structure. Owing to the compositional constraint ·1=κ, it holds dim⁡Sv=v-1. Accordingly, the description of Sv in terms of v components is redundant: an essential representation requires compositions to be properly mapped onto (v-1)-tuples. Among various possibilities, the Isometric Log-Ratio (ILR) transform introduced in [5] is the only known map of this kind leaving dA invariant. More precisely, the ILR fulfills(7)dAx,y=dEilr⁡x,ilr⁡y,dEu,v≡∑k=1v-1uk-vk2.It is known from [6] that if X~log⁡Nv(μ,Σ) is a log-normal v-variate, then C·X~Lv(μ′,Σ′) is a logistic-normal v-variate (the reader should notice that in [6] the simplex is defined by Sv={z∈R+v:z1<1}; thus property 2.2 of [6] is here reformulated so as to take into account such difference), with a known relation between (μ,Σ) and (μ′,Σ′). Analogously, it is not difficult to show that ilr⁡C·X~Nv-1(μ′′,Σ′′) is a normal (v-1)-variate, with (μ′′,Σ′′) related to (μ′,Σ′) via the change of basis matrices derived in [5]. Just to sum up, it holds(8)X~log⁡Nvμ,Σ⟹C·X~Lvμ′,Σ′⟺ilr⁡C·X~Nv-1μ′′,Σ′′.Now, suppose that (i)X fulfills (8) and has a natural interpretation as a composition, (ii) a representative set DX of observations of X is given, and (iii) we wish to select from DX those units which are compositionally closer to the center of the distribution, according to the Aitchison distance. To see that the problem is well posed, we first turn DX into a set DC·X≡{y:y=C·x, x∈DX} of compositional observations of Y=C·X. Then, we consider the special point cen⁡[Y]=C·exp⁡{E⁡[ln⁡Y]}, representing the center of the distribution of Y in a compositional sense: cen⁡[Y] minimizes the expression E[dA2(Y,cen⁡[Y])] over Sv and fulfills cen⁡[Y]=ilr-1(E⁡[ilr⁡(Y)]), see [7]. By virtue of (8) this entails ilr⁡(cen⁡[Y])=E⁡[ilr⁡(Y)]=μ′′. In order to select the observations which are closer to cen[Y], we set a threshold δ>0 and consider only those elements y∈DC·X fulfilling dA2(y,cen⁡[y])<δ, with cen⁡[y] being a sample estimator of cen⁡[Y] on DC·X. Such selection rule operates a well-defined truncation of the distribution of Y. Moreover, in view of (7) and (8), we have(9)PdA2Y,cen⁡Y<δ∣Y~Lvμ′,Σ′=PdE2Z,μ′′<δ∣Z~Nv-1μ′′,Σ′′,with Z=ilr⁡(C·X). As a consequence, we see that a compositional selection rule based on the Aitchison distance and (8) is equivalent to a spherical truncation of a multinormal distribution. Obviously, once Z has been spherically truncated, the covariance matrix of the remaining data is damped; thus an estimate of the full covariance matrix requires the reconstruction procedure described in next sections.


### 2.3. Compositional Outlier Detection via Forward Search Techniques

Outlier detection in CoDA is a practical problem where the considerations made in Section2.2 find concrete application. Most of the established methods for outlier detection make use of the Mahalanobis metric. This is however unfit to describe the distance between compositions. Log-ratio transforms allow to get rid of the compositional constraint and make it possible to apply standard statistical methods [8]. Here we consider an innovative approach, namely, the Forward Search Algorithm (FSA), introduced in [9] and thoroughly discussed in [10]. The FSA admits an elegant extension to CoDA, of which the covariance reconstruction is a key step. In the next few lines we sketch the original algorithm and outline its extension to CoDA.


#### 2.3.1. Construction of the Signal

In its standard formulation the FSA applies to normal data. It assumes a datasetDX=x(k)k=1N with x(k)∈R+(v) for k=1,…,N. The null hypothesis is that all the elements of DX are independent realizations of a stochastic variable X~N(μ0,Σ0). The FSA consists of a sequence of steps where data are recursively sorted. Along the recursion a signal is built and tested.As a preliminary step,m0 observations are randomly chosen from the bulk of DX. Let S(m0) be the set of these observations. S(m0) is used to provide initial estimates μ(m0), Σ(m0) of the true mean μ0 and the true covariance matrix Σ0, respectively. For m=m0+1,m0+2,…, the (m-m0)th step of the algorithm goes as follows:(i)
sort the elements ofDX according to the increasing values of the square Mahalanobis distance:(10)dmx2=x-μm-1TΣm-1-1x-μm-1;(ii)
take themth element x(m) of the newly sorted dataset and regard sm=dmx(m)2 as the (m-m0)th component of the signal;(iii)
letS(m) be the set of the first m observations of the newly sorted dataset;(iv)
useS(m) to provide new estimates μ(m), Σ(m) of the true mean and the true covariance matrix, respectively.Notice that S(m) is a truncated dataset. Therefore, Σ(m) must include the Tallis’ correction factor, (2) with ρ=sm. While the recursion proceeds, the inliers of DX populate progressively S(m). The recursion stops at the m-th step, when the first outlier x(m-) of DX produces a statistically relevant discontinuity in the signal.


#### 2.3.2. Statistical Test of the Signal

Statistical tests are needed to assess the relevance of discontinuities observed in the signal. At each step of the algorithm a new test is actually performed. Specifically, at themth step, sm is computed together with the lowest and highest values, respectively, δsm,α and δsm,1-α, which are theoretically admissible for sm under the null hypothesis at (1-α) significance level for some α. These values are nothing but the α- and (1-α)-percentage points of sm. Their envelopes for m-m0>0 generate two curves that surround the signal when plotted versus m. More precisely, one curve lies below it and the other lies above, provided the null hypothesis is not broken. The violation of one of the curves by the signal is interpreted as the result of the entrance of an outlier into S(m). Although the distribution of sm cannot be calculated in closed form, its percentage points are obtained from a general formula, first derived in [11], yielding(11)δsm,α=χv2-1mm+N-m-1f2N-m-1,2m;1-α,with fa,b;α the α-percentage point of the Fisher distribution with parameters (a,b). Equation (11) holds with decent approximation, as confirmed by numerical simulations.


#### 2.3.3. Extension of the Forward Search to CoDA

WhenDX is a compositional dataset, it is unnatural to assume that its elements are realizations of a multivariate X~N(μ0,Σ0). In this case the use of the FSE as outlined above does not make sense at all. Sometimes it is reasonable to assume X~Lv(μ0,Σ0), as first argued in [6]. In this case we can use the FSA to find outliers, provided that we first modify the algorithm in two respects:(i)
we replace the Mahalanobis distance by the Aitchison one;(ii)
we perform statistical tests consistently with the change of null hypothesis.Specifically, at the mth step of the algorithm, we sort DX according to the increasing values of the square Aitchison distance:(12)dAx2=12v∑i,k=1vlog⁡xixk-log⁡cm-1icm-1k2,where cm-1=cen⁡[y∣y∈S(m-1)] is the center of S(m-1). Analogously, given the mth element x(m) of the newly sorted dataset, we regard sm=dAx(m)2 as the mth component of the signal. The percentage points of sm are obtained from ilrS(m) by using the probability correspondence established in (9). Since ilr⁡(S(m)) is a spherically truncated dataset, the estimate of the covariance matrix Σ(m) derived from it must undergo the reconstruction procedure described in next sections.


### 2.4. General Covariance Reconstruction Problem

The examples discussed in the previous subsections are special cases of a more general inverse problem, namely, the reconstruction of the covariance matrixΣ of a normal multivariate X on the basis of the covariance matrix SD of X truncated to some (convex) region D. This is the simplest yet nontrivial inverse problem, which can be naturally associated with the normal distribution. The case D=Bv(ρ) corresponds to a setup where theoretical and practical aspects of the problem can be investigated with relatively modest mathematical effort. It is certainly a well-defined framework where to study regularization techniques for nonlinear inverse problems in statistics, for which there is still much room for interesting work [12, 13].


## 2.1. A Two-DimensionalGedanken Experiment in Classical Particle Physics

Consider the following ideal situation. An accelerator physicist prepares an elliptical beam of classical particles with Gaussian transversal profile. The experimenter knowsa priori the spatial distribution of the beam, that is, the covariance matrix Σ of the two-dimensional coordinates of the particles on a plane orthogonal to their flight direction. We can assume with no loss of generality that the transversal coordinate system has origin at the maximum of the beam intensity and axes along the principal components of the beam; thus it holds Σ=diag⁡(λ1,λ2). The beam travels straightforward until it enters a linear coaxial pipeline with circular profile, schematically depicted in Figure 1, where the beam is longitudinally accelerated. While the outer part of the beam is stopped by an absorbing wall, the inner part propagates within the pipeline. At the end of the beam flight the physicist wants to know if the transversal distribution of the particles is changed, due to unknown disturbance factors arisen within the pipeline. Accordingly, he measures again the spatial covariance matrix of the beam. Unfortunately, the absorbing wall has cut off the Gaussian tail, thus damping the covariance matrix and making it no more comparable to the original one. To perform such a comparison in the general case λ1≠λ2, the covariance matrix of the truncated circular beam has to go through the reconstruction procedure described in next sections.Figure 1
A classical particle beam with elliptical transversal profile is cut off upon entering a circular coaxial pipeline.


## 2.2. A Multivariate Example: Connections to Compositional Data Analysis

Compositional Data Analysis (CoDA) has been the subject of a number of papers, pioneered by Aitchison [3] over the past forty years. As a methodology of statistical investigation, it finds application in all cases where the main object of interest is a multivariate with strictly positive continuous components to be regarded as portions of a total amount κ (the normalization κ=1 is conventionally adopted in the mathematical literature). In other words, compositional variates belong to the κ-simplex:(4)Sv=z∈R+v:z1=κ,v≥2,with z1=∑k=1vzk the taxi-cab norm of z, while compositions with different norms can be always projected onto Sv via the closure operator C·x≡{κx1/x1,…,κxv/x1}. There are countless types of compositional data, whose analysis raises problems of interest for statistics [4], for example, geochemical data, balance sheet data, and election data. The simplex constraint induces a kind of dependency among the parts of a composition that goes beyond the standard concept of covariance. This invalidates many ordinary techniques of statistical analysis.In order to measure distances onSv, Aitchison introduced a positive symmetric function dA:Sv×Sv→R+, explicitly defined by(5)dAx,y=12v∑i,k=1vlog⁡xixk-log⁡yiyk2.The Aitchison distance is a key tool in CoDA. It is scale invariant in both its first and second arguments; that is, it is left invariant by redefinitions z→{αz1,…,αzv} with α∈R+. Accordingly, its support can be extended to R+v×R+v by imposing(6)dAx,y≡dAC·x,C·y,x,y∈R+v.It was proved in [5] that dA is a norm-induced metric on Sv, provided the latter is given an appropriate normed vector space structure. Owing to the compositional constraint ·1=κ, it holds dim⁡Sv=v-1. Accordingly, the description of Sv in terms of v components is redundant: an essential representation requires compositions to be properly mapped onto (v-1)-tuples. Among various possibilities, the Isometric Log-Ratio (ILR) transform introduced in [5] is the only known map of this kind leaving dA invariant. More precisely, the ILR fulfills(7)dAx,y=dEilr⁡x,ilr⁡y,dEu,v≡∑k=1v-1uk-vk2.It is known from [6] that if X~log⁡Nv(μ,Σ) is a log-normal v-variate, then C·X~Lv(μ′,Σ′) is a logistic-normal v-variate (the reader should notice that in [6] the simplex is defined by Sv={z∈R+v:z1<1}; thus property 2.2 of [6] is here reformulated so as to take into account such difference), with a known relation between (μ,Σ) and (μ′,Σ′). Analogously, it is not difficult to show that ilr⁡C·X~Nv-1(μ′′,Σ′′) is a normal (v-1)-variate, with (μ′′,Σ′′) related to (μ′,Σ′) via the change of basis matrices derived in [5]. Just to sum up, it holds(8)X~log⁡Nvμ,Σ⟹C·X~Lvμ′,Σ′⟺ilr⁡C·X~Nv-1μ′′,Σ′′.Now, suppose that (i)X fulfills (8) and has a natural interpretation as a composition, (ii) a representative set DX of observations of X is given, and (iii) we wish to select from DX those units which are compositionally closer to the center of the distribution, according to the Aitchison distance. To see that the problem is well posed, we first turn DX into a set DC·X≡{y:y=C·x, x∈DX} of compositional observations of Y=C·X. Then, we consider the special point cen⁡[Y]=C·exp⁡{E⁡[ln⁡Y]}, representing the center of the distribution of Y in a compositional sense: cen⁡[Y] minimizes the expression E[dA2(Y,cen⁡[Y])] over Sv and fulfills cen⁡[Y]=ilr-1(E⁡[ilr⁡(Y)]), see [7]. By virtue of (8) this entails ilr⁡(cen⁡[Y])=E⁡[ilr⁡(Y)]=μ′′. In order to select the observations which are closer to cen[Y], we set a threshold δ>0 and consider only those elements y∈DC·X fulfilling dA2(y,cen⁡[y])<δ, with cen⁡[y] being a sample estimator of cen⁡[Y] on DC·X. Such selection rule operates a well-defined truncation of the distribution of Y. Moreover, in view of (7) and (8), we have(9)PdA2Y,cen⁡Y<δ∣Y~Lvμ′,Σ′=PdE2Z,μ′′<δ∣Z~Nv-1μ′′,Σ′′,with Z=ilr⁡(C·X). As a consequence, we see that a compositional selection rule based on the Aitchison distance and (8) is equivalent to a spherical truncation of a multinormal distribution. Obviously, once Z has been spherically truncated, the covariance matrix of the remaining data is damped; thus an estimate of the full covariance matrix requires the reconstruction procedure described in next sections.


## 2.3. Compositional Outlier Detection via Forward Search Techniques

Outlier detection in CoDA is a practical problem where the considerations made in Section2.2 find concrete application. Most of the established methods for outlier detection make use of the Mahalanobis metric. This is however unfit to describe the distance between compositions. Log-ratio transforms allow to get rid of the compositional constraint and make it possible to apply standard statistical methods [8]. Here we consider an innovative approach, namely, the Forward Search Algorithm (FSA), introduced in [9] and thoroughly discussed in [10]. The FSA admits an elegant extension to CoDA, of which the covariance reconstruction is a key step. In the next few lines we sketch the original algorithm and outline its extension to CoDA.


### 2.3.1. Construction of the Signal

In its standard formulation the FSA applies to normal data. It assumes a datasetDX=x(k)k=1N with x(k)∈R+(v) for k=1,…,N. The null hypothesis is that all the elements of DX are independent realizations of a stochastic variable X~N(μ0,Σ0). The FSA consists of a sequence of steps where data are recursively sorted. Along the recursion a signal is built and tested.As a preliminary step,m0 observations are randomly chosen from the bulk of DX. Let S(m0) be the set of these observations. S(m0) is used to provide initial estimates μ(m0), Σ(m0) of the true mean μ0 and the true covariance matrix Σ0, respectively. For m=m0+1,m0+2,…, the (m-m0)th step of the algorithm goes as follows:(i)
sort the elements ofDX according to the increasing values of the square Mahalanobis distance:(10)dmx2=x-μm-1TΣm-1-1x-μm-1;(ii)
take themth element x(m) of the newly sorted dataset and regard sm=dmx(m)2 as the (m-m0)th component of the signal;(iii)
letS(m) be the set of the first m observations of the newly sorted dataset;(iv)
useS(m) to provide new estimates μ(m), Σ(m) of the true mean and the true covariance matrix, respectively.Notice that S(m) is a truncated dataset. Therefore, Σ(m) must include the Tallis’ correction factor, (2) with ρ=sm. While the recursion proceeds, the inliers of DX populate progressively S(m). The recursion stops at the m-th step, when the first outlier x(m-) of DX produces a statistically relevant discontinuity in the signal.


### 2.3.2. Statistical Test of the Signal

Statistical tests are needed to assess the relevance of discontinuities observed in the signal. At each step of the algorithm a new test is actually performed. Specifically, at themth step, sm is computed together with the lowest and highest values, respectively, δsm,α and δsm,1-α, which are theoretically admissible for sm under the null hypothesis at (1-α) significance level for some α. These values are nothing but the α- and (1-α)-percentage points of sm. Their envelopes for m-m0>0 generate two curves that surround the signal when plotted versus m. More precisely, one curve lies below it and the other lies above, provided the null hypothesis is not broken. The violation of one of the curves by the signal is interpreted as the result of the entrance of an outlier into S(m). Although the distribution of sm cannot be calculated in closed form, its percentage points are obtained from a general formula, first derived in [11], yielding(11)δsm,α=χv2-1mm+N-m-1f2N-m-1,2m;1-α,with fa,b;α the α-percentage point of the Fisher distribution with parameters (a,b). Equation (11) holds with decent approximation, as confirmed by numerical simulations.


### 2.3.3. Extension of the Forward Search to CoDA

WhenDX is a compositional dataset, it is unnatural to assume that its elements are realizations of a multivariate X~N(μ0,Σ0). In this case the use of the FSE as outlined above does not make sense at all. Sometimes it is reasonable to assume X~Lv(μ0,Σ0), as first argued in [6]. In this case we can use the FSA to find outliers, provided that we first modify the algorithm in two respects:(i)
we replace the Mahalanobis distance by the Aitchison one;(ii)
we perform statistical tests consistently with the change of null hypothesis.Specifically, at the mth step of the algorithm, we sort DX according to the increasing values of the square Aitchison distance:(12)dAx2=12v∑i,k=1vlog⁡xixk-log⁡cm-1icm-1k2,where cm-1=cen⁡[y∣y∈S(m-1)] is the center of S(m-1). Analogously, given the mth element x(m) of the newly sorted dataset, we regard sm=dAx(m)2 as the mth component of the signal. The percentage points of sm are obtained from ilrS(m) by using the probability correspondence established in (9). Since ilr⁡(S(m)) is a spherically truncated dataset, the estimate of the covariance matrix Σ(m) derived from it must undergo the reconstruction procedure described in next sections.


## 2.3.1. Construction of the Signal

In its standard formulation the FSA applies to normal data. It assumes a datasetDX=x(k)k=1N with x(k)∈R+(v) for k=1,…,N. The null hypothesis is that all the elements of DX are independent realizations of a stochastic variable X~N(μ0,Σ0). The FSA consists of a sequence of steps where data are recursively sorted. Along the recursion a signal is built and tested.As a preliminary step,m0 observations are randomly chosen from the bulk of DX. Let S(m0) be the set of these observations. S(m0) is used to provide initial estimates μ(m0), Σ(m0) of the true mean μ0 and the true covariance matrix Σ0, respectively. For m=m0+1,m0+2,…, the (m-m0)th step of the algorithm goes as follows:(i)
sort the elements ofDX according to the increasing values of the square Mahalanobis distance:(10)dmx2=x-μm-1TΣm-1-1x-μm-1;(ii)
take themth element x(m) of the newly sorted dataset and regard sm=dmx(m)2 as the (m-m0)th component of the signal;(iii)
letS(m) be the set of the first m observations of the newly sorted dataset;(iv)
useS(m) to provide new estimates μ(m), Σ(m) of the true mean and the true covariance matrix, respectively.Notice that S(m) is a truncated dataset. Therefore, Σ(m) must include the Tallis’ correction factor, (2) with ρ=sm. While the recursion proceeds, the inliers of DX populate progressively S(m). The recursion stops at the m-th step, when the first outlier x(m-) of DX produces a statistically relevant discontinuity in the signal.


## 2.3.2. Statistical Test of the Signal

Statistical tests are needed to assess the relevance of discontinuities observed in the signal. At each step of the algorithm a new test is actually performed. Specifically, at themth step, sm is computed together with the lowest and highest values, respectively, δsm,α and δsm,1-α, which are theoretically admissible for sm under the null hypothesis at (1-α) significance level for some α. These values are nothing but the α- and (1-α)-percentage points of sm. Their envelopes for m-m0>0 generate two curves that surround the signal when plotted versus m. More precisely, one curve lies below it and the other lies above, provided the null hypothesis is not broken. The violation of one of the curves by the signal is interpreted as the result of the entrance of an outlier into S(m). Although the distribution of sm cannot be calculated in closed form, its percentage points are obtained from a general formula, first derived in [11], yielding(11)δsm,α=χv2-1mm+N-m-1f2N-m-1,2m;1-α,with fa,b;α the α-percentage point of the Fisher distribution with parameters (a,b). Equation (11) holds with decent approximation, as confirmed by numerical simulations.


## 2.3.3. Extension of the Forward Search to CoDA

WhenDX is a compositional dataset, it is unnatural to assume that its elements are realizations of a multivariate X~N(μ0,Σ0). In this case the use of the FSE as outlined above does not make sense at all. Sometimes it is reasonable to assume X~Lv(μ0,Σ0), as first argued in [6]. In this case we can use the FSA to find outliers, provided that we first modify the algorithm in two respects:(i)
we replace the Mahalanobis distance by the Aitchison one;(ii)
we perform statistical tests consistently with the change of null hypothesis.Specifically, at the mth step of the algorithm, we sort DX according to the increasing values of the square Aitchison distance:(12)dAx2=12v∑i,k=1vlog⁡xixk-log⁡cm-1icm-1k2,where cm-1=cen⁡[y∣y∈S(m-1)] is the center of S(m-1). Analogously, given the mth element x(m) of the newly sorted dataset, we regard sm=dAx(m)2 as the mth component of the signal. The percentage points of sm are obtained from ilrS(m) by using the probability correspondence established in (9). Since ilr⁡(S(m)) is a spherically truncated dataset, the estimate of the covariance matrix Σ(m) derived from it must undergo the reconstruction procedure described in next sections.


## 2.4. General Covariance Reconstruction Problem

The examples discussed in the previous subsections are special cases of a more general inverse problem, namely, the reconstruction of the covariance matrixΣ of a normal multivariate X on the basis of the covariance matrix SD of X truncated to some (convex) region D. This is the simplest yet nontrivial inverse problem, which can be naturally associated with the normal distribution. The case D=Bv(ρ) corresponds to a setup where theoretical and practical aspects of the problem can be investigated with relatively modest mathematical effort. It is certainly a well-defined framework where to study regularization techniques for nonlinear inverse problems in statistics, for which there is still much room for interesting work [12, 13].


## 3. Definitions and Setup

LetX∈Rv be a random vector with jointly normal distribution Nv(0,Σ) in v≥1 dimensions. The probability that X falls within Bv(ρ) is measured by the Gaussian integral(13)αρ;Σ≡PX∈Bvρ=12πv/2Σ1/2∫Bvρdvxexp⁡-12xTΣ-1x.Since Σ is symmetric positive definite, it has orthonormal eigenvectors Σv(i)=λiv(i). Let us denote by R≡vi(j)i,j=1v the orthogonal matrix having these vectors as columns and by Λ≡diag⁡(λ)=RTΣR the diagonal counterpart of Σ. From the invariance of Bv(ρ) under rotations, it follows that α depends upon Σ just by way of λ. Accordingly, we rename the Gaussian probability content of Bv(ρ) as(14)αρ;λ≡∫Bvρdvx∏m=1vδxm,λm,δy,η=12πηexp⁡-y22η.Note that (14) is not sufficient to fully characterize the random vector X under the spherical constraint, for which we need to calculate the distribution law P[X∈A∣X∈Bv(ρ)] as a function of A⊂Rv. Alternatively, we can describe X in terms of the complete set of its truncated moments(15)mk1⋯kvρ;Σ≡EX1k1⋯XvkvX∈Bvρ,kii=1,…,v∈Nv.As usual, these can be all obtained from the moment-generating function(16)αmt=12πv/2Σ1/2∫Bvρdvxexp⁡tTx-12xTΣ-1x,t∈Rv,by differentiating the latter an arbitrary number of times with respect to the components of t; namely,(17)mk1⋯kvρ;Σ=∂k1+⋯+kvmt∂t1k1⋯∂tvkvt=0.It will be observed that m(t) is in general not invariant under rotations of t. Therefore, unlike α, the moments mk1⋯kv depend effectively on both λ and R. For instance, for the matrix of the second moments SB≡∂2m/∂ti∂tjt=0i,j=1v such dependence amounts to(18)αSBij=∑k,l=1vRkiRlj∫Bvρdvxxkxl∏m=1vδxm,λm.By parity, the only nonvanishing terms in the above sum are those with k=l. Hence, it follows that Σ and SB share R as a common diagonalizing matrix. In other words, if M≡diag⁡(μ) is the diagonal matrix of the eigenvalues of SB, then M=RTSBR. Moreover, μk is related to λk by(19)μk=λkαkα,αkρ;λ≡∫Bvρdvxxk2λk∏m=1vδxm,λm,k=1,…,v.The ratios αk/α are naturally interpreted as adimensional correction factors to the eigenvalues of Σ, so they play the same role as cT(ρ) in (2). However, αk/α depends explicitly on the subscript k; thus each eigenvalue is damped differently from the others as a consequence of the condition X∈Bv(ρ).Remark 1.
In practical terms, (18)-(19) tell us that estimating the sample covariance matrix of X~Nv(0,Σ) from a spherically truncated population x(m)m=1M, made of M independent units, via the classical estimator Qij=M-1-1∑m=1Mxim-x~imxjm-x~jm, being x~=M-1∑m=1Mx(m) the sample mean, yields a damped result. Nonetheless, the damping affects only the eigenvalues of the estimator, whereas its eigenvectors are left invariant.


### 3.1. Monotonicity Properties of Ratios of Gaussian Integrals

Equations (14) and (19) suggest introducing a general notation for the Gaussian integrals over Bv(ρ), under the assumption Σ=Λ. So, we define(20)αklm⋯ρ;λ≡∫Bvρdvxxk2λkxl2λlxm2λm⋯∏n=1vδxn,λn,with each subscript q on theleft-hand side addressing an additional factor of xq2/λq under the integral sign on theright-hand side (no subscripts means α). Several analytic properties of such integrals are discussed in [14]. In the following proposition, we lay emphasis on some issues concerning the monotonicity trends of the ratios αk/α.Proposition 2 (monotonicities).
Letλ(k)≡λii=1,…,vi≠k denote the set of the full eigenvalues without λk. The ratios αk/α fulfill the following properties:(p1)
λkαk/α(ρ;λ) is a monotonic increasing function of λk at fixed ρ and λ(k);(p2)
αk/α(ρ;λ) is a monotonic decreasing function of λk at fixed ρ and λ(k);(p3)
αk/α(ρ;λ) is a monotonic decreasing function of λi(i≠k) at fixed ρ and λ(i),where an innocuous abuse of notation has been made on writing αk/α(ρ;λ) in place of αk(ρ;λ)/α(ρ;λ).Proof.
Let the symbol∂k≡∂/∂λk denote a derivative with respect to λk. In order to prove property (p1), we apply the chain rule of differentiation to λkαk/α and then we pass ∂k under the integral sign in ∂kα and ∂kαk. In this way, we obtain(21)∂kλkαkα=12αkkα-αk2α2=12λk2EXk4∣X∈Bvρ-EXk2∣X∈Bvρ2=12λk2var⁡Xk2∣X∈Bvρ≥0.Moreover, since the truncated marginal density of Xk2 is positive within a set of nonzero measure in R, the monotonic trend of λkαk/α in λk is strict.
Properties (p2) and (p3) are less trivial than (p1). Indeed, the same reasoning as above now yields on the one hand(22)λk∂kαkα=∂kλkαkα-αkα=12λk2var⁡Xk2∣X∈Bvρ-2λkEXk2∣X∈Bvρ≤0,and on the other(23)λi∂iαkα=12αikα-αiαkα2=12λiλkcov⁡Xi2,Xk2∣X∈Bvρ≤0i≠k.Despite being nota priori evident, theright-hand side of both (22) and (23) is negative (and vanishes in the limit ρ→∞). The inequalities var⁡(Xk2)≤2λkEXk2 within Euclidean balls have been first discussed in [14], while the inequalities cov(Xj2,Xk2) within generalized Orlicz balls have been discussed in [15, 16] for the case where the probability distribution of X is flat instead of being normal. More recently, a complete proof of both inequalities has been given in [17]. Despite the technical difficulties in proving them, their meaning should be intuitively clear. The variance inequality quantifies the squeezing affecting Xk2 as a consequence of the truncation (in the unconstrained case it would be var⁡(Xk2)=2λk2). The covariance inequality follows from the opposition arising among the square components in proximity of the boundary of Bv(ρ). Indeed, if Xj2↗ρ, then Xk2↘0∀k≠j in order for X to stay within Bv(ρ).


### 3.2. Definition Domain of the Reconstruction Problem

A consequence of Proposition2 is represented by the following.Corollary 3.
Givenv, ρ, and λ, μk is bounded by(24)ρrv,ρ/2λk≤μk≤ρ3,rv,z≡2v+1Mv,v+1/2,zMv,v+3/2,z,with M denoting the Kummer function; namely,(25)Ma,b,z=∑n=0∞1n!anbnzn,xn≡Γx+nΓx.Proof.
The upper bound of (24) corresponds to the value of μk in the v-tuple limit λk→∞, λ(k)→{0,…,0}. This is indeed the maximum possible value allowed for μk according to properties (p1) and (p3) of Proposition 2. In order to perform this limit, we observe that(26)limη→0+⁡δy,η=δy,with the δ symbol on theright-hand side representing the Dirac delta function (the reader who is not familiar with the theory of distributions may refer, for instance, to [18] for an introduction). Accordingly,(27)limλk→∞⁡limλk→0,…,0⁡μk=∫-ρ+ρdxkxk2∫-ρ+ρdxk=ρ3.The lower bound corresponds instead to the value taken by μk as λ(k)→{∞,…,∞} and λk is kept fixed. In this limit, all the Gaussian factors in the probability density function except the kth one flatten to one. Hence,(28)limλk→∞,…,∞⁡μk=limλk→∞,…,∞⁡∫-ρ+ρdxkxk2δxk,λkαv-1ρ-xk2;λk∫-ρ+ρdxkδxk,λkαv-1ρ-xk2;λk=∫-ρ+ρdxkxk2e-xk2/2λkρ-xk2v-1∫-ρ+ρdxke-xk2/2λkρ-xk2v-1=ρ∫01dxkxk2e-ρ/2λkxk21-xk2v-1∫01dxke-ρ/2λkxk21-xk2v-1.Numerator and denominator of the rightmost ratio are easily recognized to be integral representations of Kummer functions (see, e.g., [19, ch. 13]).The upper bound of (24) can be sharpened, as clarified by the following.Proposition 4 (bounds on the truncated moments).
Letv, ρ, and λ be given. If i1,…,iv is a permutation of 1,…,v such that μi1≤μi2≤⋯≤μiv, then the following upper bounds hold:(29)i∑k=1vμk≤ρ;iiμik≤ρv-k+1,k=1,…,v.Proof.
The overall upper bound on the sum of truncated moments follows from(30)∑k=1vμk=1α∑k=1vλkαk=1α∫Bvρdvx∑k=1vxk2∏m=1vδxm,λm≤ρ.At the same time, the sum can be split and bounded from below by(31)∑k=1vμik=∑k=1nμik+∑k=n+1vμik≥∑k=1nμik+v-nμin+1,n=0,1,…,v-1.The single upper bounds on the μk’s are then obtained from (30)-(31). It will be noted that (29) ii is sharper than the upper bound of (24) only for v>3 and k<v-2.From now on, we shall assume, with no loss of generality, that the eigenvalues ofΣ are increasingly ordered, namely, 0<λ1≤⋯≤λv (we can always permute the labels of the coordinate axes, so as to let this be the case). An important aspect related to the eigenvalue ordering is provided by the following.Proposition 5 (eigenvalue ordering).
Letv, ρ, and λ be given. If λ1≤λ2≤⋯≤λv, then μ1≤μ2≤⋯≤μv holds as well.Proof.
In order to show that the spherical truncation does not violate the eigenvalue ordering, we make repeated use of the monotonicity properties of Proposition2. Specifically, if i<j, then(32)μi=λiαiαρ;λ1,…,λi,…,λj,…,λv≤λjαiαρ;λ1,…,λj,…,λj,…,λv↭increasing  monotonicity  ofλiαiα=λjαjαρ;λ1,…,λj,…,λj,…,λv↭exchange  symmetryi⟷j≤λjαjαρ;λ1,…,λi,…,λj,…,λv↭decreasing  monotonicity  ofαjα=μj,where the symbol “↭” is used to explain where the inequality sign preceding it comes from and the “exchange symmetry” refers to the formal property of the one-index Gaussian integrals over Bv(ρ) to fulfill αi(ρ;{λ1,…,λi,…,λj,…,λv})=αj(ρ;{λ1,…,λj,…,λi,…,λv}).Let us now focus on (19). They have to be solved in order to reconstruct λ from μ. Formally, if we introduce a family of truncation operators τρ:R+v→R+v (parametrically depending on ρ), such that(33)τρ·λk≡λkαkαρ;λ,k=1,…,v,then the reconstruction of λ from μ amounts to calculating λ=τρ-1·μ. One should be aware that τρ is not a surjective operator in view of Corollary 3 and Proposition 4. Therefore, τρ-1 is only defined within a bounded domain Dτρ-1. If we define(34)D0=μ∈R+v:μ1≤⋯≤μv,μk=λkαkαfork=1,…,v,forsomeλ∈R+v,then we have D(τρ-1)={μ:μ=σ·μ0forsomeμ0∈D0,σ∈Sv}, where Sv is the set of permutations of v elements. From Proposition 4 we conclude that D0⊆Hv(ρ), being(35)Hvρ≡x∈R+v:xk≤min⁡ρ3,ρv-k+1,∑k=1vxk≤ρ,∀k.In fact, there are vectors μ∈R+v fulfilling μ∈Hv(ρ) and μ∉D0; thus we conclude that D0 is a proper subset of Hv(ρ). Numerical experiences based on the techniques discussed in the next sections show indeed that(36)Dτρ-1=⋂k=1vμ∈R+v:∑j≠k1⋯vμj+3μk≤ρ.A graphical representation of (36) in v=2 and v=3 dimensions is depicted in Figure 2. The reader should note that until Section 7 we shall always assume that μ comes from the application of τρ to some λ; thus μ∈D(τρ-1) by construction.Figure 2
(a) Numerical reconstruction ofD(τρ-1) in v=2 dimensions. (b) Numerical reconstruction of D(τρ-1) in v=3 dimensions.
(a) (b)Now, we observe that (19) can be written in the equivalent form(37)λ=Tλ;μ;ρ,(38)T:R+v×R+v×R+⟶R+v;Tkλ;μ;ρ=μkααkρ;λ,k=1,…,v.Since ρ and μ are (nonindependent) input parameters for the covariance reconstruction problem (and in order to keep the notation light), in the sequel we shall leave the dependence of T upon ρ and μ implicitly understood; that is, we shall write (37) as λ=T(λ). Hence, we see that the full eigenvalue spectrum λ is a fixed point of the operator T. This suggests obtaining it as the limit of a sequence(39)λ0=μ,λn+1=Tλn,n=0,1,…,(40)λ=limn→∞⁡λn,provided that this can be shown to converge. Note that since αk<α, it follows that Tk(λ(n))>μk∀n, so the sequence is bounded from below by μ. In particular, this holds for n=0. Therefore, the sequence moves to the right direction at least at the beginning. A formal proof of convergence, based on the monotonicity properties stated by Proposition 2, is given in the next section.


## 3.1. Monotonicity Properties of Ratios of Gaussian Integrals

Equations (14) and (19) suggest introducing a general notation for the Gaussian integrals over Bv(ρ), under the assumption Σ=Λ. So, we define(20)αklm⋯ρ;λ≡∫Bvρdvxxk2λkxl2λlxm2λm⋯∏n=1vδxn,λn,with each subscript q on theleft-hand side addressing an additional factor of xq2/λq under the integral sign on theright-hand side (no subscripts means α). Several analytic properties of such integrals are discussed in [14]. In the following proposition, we lay emphasis on some issues concerning the monotonicity trends of the ratios αk/α.Proposition 2 (monotonicities).
Letλ(k)≡λii=1,…,vi≠k denote the set of the full eigenvalues without λk. The ratios αk/α fulfill the following properties:(p1)
λkαk/α(ρ;λ) is a monotonic increasing function of λk at fixed ρ and λ(k);(p2)
αk/α(ρ;λ) is a monotonic decreasing function of λk at fixed ρ and λ(k);(p3)
αk/α(ρ;λ) is a monotonic decreasing function of λi(i≠k) at fixed ρ and λ(i),where an innocuous abuse of notation has been made on writing αk/α(ρ;λ) in place of αk(ρ;λ)/α(ρ;λ).Proof.
Let the symbol∂k≡∂/∂λk denote a derivative with respect to λk. In order to prove property (p1), we apply the chain rule of differentiation to λkαk/α and then we pass ∂k under the integral sign in ∂kα and ∂kαk. In this way, we obtain(21)∂kλkαkα=12αkkα-αk2α2=12λk2EXk4∣X∈Bvρ-EXk2∣X∈Bvρ2=12λk2var⁡Xk2∣X∈Bvρ≥0.Moreover, since the truncated marginal density of Xk2 is positive within a set of nonzero measure in R, the monotonic trend of λkαk/α in λk is strict.
Properties (p2) and (p3) are less trivial than (p1). Indeed, the same reasoning as above now yields on the one hand(22)λk∂kαkα=∂kλkαkα-αkα=12λk2var⁡Xk2∣X∈Bvρ-2λkEXk2∣X∈Bvρ≤0,and on the other(23)λi∂iαkα=12αikα-αiαkα2=12λiλkcov⁡Xi2,Xk2∣X∈Bvρ≤0i≠k.Despite being nota priori evident, theright-hand side of both (22) and (23) is negative (and vanishes in the limit ρ→∞). The inequalities var⁡(Xk2)≤2λkEXk2 within Euclidean balls have been first discussed in [14], while the inequalities cov(Xj2,Xk2) within generalized Orlicz balls have been discussed in [15, 16] for the case where the probability distribution of X is flat instead of being normal. More recently, a complete proof of both inequalities has been given in [17]. Despite the technical difficulties in proving them, their meaning should be intuitively clear. The variance inequality quantifies the squeezing affecting Xk2 as a consequence of the truncation (in the unconstrained case it would be var⁡(Xk2)=2λk2). The covariance inequality follows from the opposition arising among the square components in proximity of the boundary of Bv(ρ). Indeed, if Xj2↗ρ, then Xk2↘0∀k≠j in order for X to stay within Bv(ρ).


## 3.2. Definition Domain of the Reconstruction Problem

A consequence of Proposition2 is represented by the following.Corollary 3.
Givenv, ρ, and λ, μk is bounded by(24)ρrv,ρ/2λk≤μk≤ρ3,rv,z≡2v+1Mv,v+1/2,zMv,v+3/2,z,with M denoting the Kummer function; namely,(25)Ma,b,z=∑n=0∞1n!anbnzn,xn≡Γx+nΓx.Proof.
The upper bound of (24) corresponds to the value of μk in the v-tuple limit λk→∞, λ(k)→{0,…,0}. This is indeed the maximum possible value allowed for μk according to properties (p1) and (p3) of Proposition 2. In order to perform this limit, we observe that(26)limη→0+⁡δy,η=δy,with the δ symbol on theright-hand side representing the Dirac delta function (the reader who is not familiar with the theory of distributions may refer, for instance, to [18] for an introduction). Accordingly,(27)limλk→∞⁡limλk→0,…,0⁡μk=∫-ρ+ρdxkxk2∫-ρ+ρdxk=ρ3.The lower bound corresponds instead to the value taken by μk as λ(k)→{∞,…,∞} and λk is kept fixed. In this limit, all the Gaussian factors in the probability density function except the kth one flatten to one. Hence,(28)limλk→∞,…,∞⁡μk=limλk→∞,…,∞⁡∫-ρ+ρdxkxk2δxk,λkαv-1ρ-xk2;λk∫-ρ+ρdxkδxk,λkαv-1ρ-xk2;λk=∫-ρ+ρdxkxk2e-xk2/2λkρ-xk2v-1∫-ρ+ρdxke-xk2/2λkρ-xk2v-1=ρ∫01dxkxk2e-ρ/2λkxk21-xk2v-1∫01dxke-ρ/2λkxk21-xk2v-1.Numerator and denominator of the rightmost ratio are easily recognized to be integral representations of Kummer functions (see, e.g., [19, ch. 13]).The upper bound of (24) can be sharpened, as clarified by the following.Proposition 4 (bounds on the truncated moments).
Letv, ρ, and λ be given. If i1,…,iv is a permutation of 1,…,v such that μi1≤μi2≤⋯≤μiv, then the following upper bounds hold:(29)i∑k=1vμk≤ρ;iiμik≤ρv-k+1,k=1,…,v.Proof.
The overall upper bound on the sum of truncated moments follows from(30)∑k=1vμk=1α∑k=1vλkαk=1α∫Bvρdvx∑k=1vxk2∏m=1vδxm,λm≤ρ.At the same time, the sum can be split and bounded from below by(31)∑k=1vμik=∑k=1nμik+∑k=n+1vμik≥∑k=1nμik+v-nμin+1,n=0,1,…,v-1.The single upper bounds on the μk’s are then obtained from (30)-(31). It will be noted that (29) ii is sharper than the upper bound of (24) only for v>3 and k<v-2.From now on, we shall assume, with no loss of generality, that the eigenvalues ofΣ are increasingly ordered, namely, 0<λ1≤⋯≤λv (we can always permute the labels of the coordinate axes, so as to let this be the case). An important aspect related to the eigenvalue ordering is provided by the following.Proposition 5 (eigenvalue ordering).
Letv, ρ, and λ be given. If λ1≤λ2≤⋯≤λv, then μ1≤μ2≤⋯≤μv holds as well.Proof.
In order to show that the spherical truncation does not violate the eigenvalue ordering, we make repeated use of the monotonicity properties of Proposition2. Specifically, if i<j, then(32)μi=λiαiαρ;λ1,…,λi,…,λj,…,λv≤λjαiαρ;λ1,…,λj,…,λj,…,λv↭increasing  monotonicity  ofλiαiα=λjαjαρ;λ1,…,λj,…,λj,…,λv↭exchange  symmetryi⟷j≤λjαjαρ;λ1,…,λi,…,λj,…,λv↭decreasing  monotonicity  ofαjα=μj,where the symbol “↭” is used to explain where the inequality sign preceding it comes from and the “exchange symmetry” refers to the formal property of the one-index Gaussian integrals over Bv(ρ) to fulfill αi(ρ;{λ1,…,λi,…,λj,…,λv})=αj(ρ;{λ1,…,λj,…,λi,…,λv}).Let us now focus on (19). They have to be solved in order to reconstruct λ from μ. Formally, if we introduce a family of truncation operators τρ:R+v→R+v (parametrically depending on ρ), such that(33)τρ·λk≡λkαkαρ;λ,k=1,…,v,then the reconstruction of λ from μ amounts to calculating λ=τρ-1·μ. One should be aware that τρ is not a surjective operator in view of Corollary 3 and Proposition 4. Therefore, τρ-1 is only defined within a bounded domain Dτρ-1. If we define(34)D0=μ∈R+v:μ1≤⋯≤μv,μk=λkαkαfork=1,…,v,forsomeλ∈R+v,then we have D(τρ-1)={μ:μ=σ·μ0forsomeμ0∈D0,σ∈Sv}, where Sv is the set of permutations of v elements. From Proposition 4 we conclude that D0⊆Hv(ρ), being(35)Hvρ≡x∈R+v:xk≤min⁡ρ3,ρv-k+1,∑k=1vxk≤ρ,∀k.In fact, there are vectors μ∈R+v fulfilling μ∈Hv(ρ) and μ∉D0; thus we conclude that D0 is a proper subset of Hv(ρ). Numerical experiences based on the techniques discussed in the next sections show indeed that(36)Dτρ-1=⋂k=1vμ∈R+v:∑j≠k1⋯vμj+3μk≤ρ.A graphical representation of (36) in v=2 and v=3 dimensions is depicted in Figure 2. The reader should note that until Section 7 we shall always assume that μ comes from the application of τρ to some λ; thus μ∈D(τρ-1) by construction.Figure 2
(a) Numerical reconstruction ofD(τρ-1) in v=2 dimensions. (b) Numerical reconstruction of D(τρ-1) in v=3 dimensions.
(a) (b)Now, we observe that (19) can be written in the equivalent form(37)λ=Tλ;μ;ρ,(38)T:R+v×R+v×R+⟶R+v;Tkλ;μ;ρ=μkααkρ;λ,k=1,…,v.Since ρ and μ are (nonindependent) input parameters for the covariance reconstruction problem (and in order to keep the notation light), in the sequel we shall leave the dependence of T upon ρ and μ implicitly understood; that is, we shall write (37) as λ=T(λ). Hence, we see that the full eigenvalue spectrum λ is a fixed point of the operator T. This suggests obtaining it as the limit of a sequence(39)λ0=μ,λn+1=Tλn,n=0,1,…,(40)λ=limn→∞⁡λn,provided that this can be shown to converge. Note that since αk<α, it follows that Tk(λ(n))>μk∀n, so the sequence is bounded from below by μ. In particular, this holds for n=0. Therefore, the sequence moves to the right direction at least at the beginning. A formal proof of convergence, based on the monotonicity properties stated by Proposition 2, is given in the next section.


## 4. Convergence of the Fixed Point Equation

We split our argument into three propositions, describing different properties of the sequenceλn. They assert, respectively, that (i) the sequence is componentwise monotically increasing; (ii) the sequence is componentwise bounded from above by any fixed point of T; and (iii) if T has a fixed point, this must be unique. Statements (i) and (ii) are sufficient to guarantee the convergence of the sequence to a finite limit (the unconstrained spectrum is a fixed point of T). In addition, the limit is easily recognized to be a fixed point of T. Hence, statement (iii) guarantees that the sequence converges to the unconstrained eigenvalue spectrum. We remark that all the monotonicities discussed in Proposition 2 are strict; that is, the ratios αk/α have no stationary points at finite ρ and λ, which is crucial for the proof.Proposition 6 (increasing monotonicity).
Givenv, ρ, and μ∈D(τρ-1), the sequence λ(0)=μ, λ(n+1)=T(λ(n)), n=0,1,…, is monotically increasing; namely, λk(n+1)>λk(n)∀k=1,…,v.Proof.
The proof is by induction. We first notice that(41)λk1=Tkλ0=Tkμ=μkααkρ;μ>μk=λk0,k=1,…,v;the inequality follows from αk(ρ;μ)<α(ρ;μ). Suppose now that the property of increasing monotonicity has been checked off up to thenth element of the sequence. Then,(42)λkn+1=μkααkρ;λn>μkααkρ;λn-1=λkn;the inequality follows this time from the inductive hypothesis and from properties (p2) and (p3) of Proposition 2.Proposition 7 (boundedness).
Givenv, ρ, and μ∈D(τρ-1), the sequence λ(0)=μ, λ(n+1)=T(λ(n)), n=0,1,…, is bounded from above; namely, λk(n)<λk∗∀k=1,…,v, λ∗ being a fixed point of T.Proof.
We proceed again by induction. We first notice that(43)λk0=μk<μkααkρ;λ∗=λk∗,k=1,…,v;the inequality follows as previously from αk(ρ;λ∗)<α(ρ;λ∗). Suppose now that the property of boundedness has been checked off up to thenth element of the sequence. Then,(44)λkn+1=μkααkρ;λn=λk∗αkαρ;λ∗ααkρ;λn<λk∗αkαρ;λnααkρ;λn=λk∗;the inequality follows for the last time from the inductive hypothesis and from properties (p2) and (p3) of Proposition 2.According to Propositions6 and 7, the sequence converges. Now, let λ~=limn→∞λ(n) be the limit of the sequence. Effortlessly, we prove that λ~ is a fixed point of T. Indeed,(45)λ~k=limn→∞⁡λkn=limn→∞⁡μkααkρ;λn-1=μkααkρ;limn→∞⁡λn-1=Tkλ~.Note that passing the limit over n under the integral sign is certainly allowed for Gaussian integrals.Proposition 8 (uniqueness of the fixed point).
Letλ′=T(λ′) and λ′′=T(λ′′) be two fixed points of T, corresponding to the same choice of v, ρ, and μ∈D(τρ-1). Then, it must be that λ′=λ′′.Proof.
According to the hypothesis,λ′ and λ′′ fulfill the equations(46)λk′=μkααkρ;λ′⟹μk=λk′αkαρ;λ′,λk′′=μkααkρ;λ′′⟹μk=λk′′αkαρ;λ′′.Hence,(47)0=μk-μk=λk′αkαρ;λ′-λk′′αkαρ;λ′′=∑l=1v∫01dtJklρ;λ′′+tλ′-λ′′λl′-λl′′,where J denotes the Jacobian matrix of τρ and is given by (48)Jklρ;λ=∂kλlαlαρ;λ=12λlλkαklα-αkαlα2=Λ-1Ωρ;λΛkl,having set Ωkl≡(1/2)(αkl/α-αkαl/α2). It will be noted that Ω={Ωkl}k,l=1v is essentially the covariance matrix of the square components of X under spherical truncation (we have come across its matrix elements in (21)–(23)). As such, Ω is symmetric positive definite. Indeed,(49)Ωkl=12λkλlcov⁡Xk2,Xl2∣X∈Bvρ=12λkλlEXk2-EXk2∣X∈BvρXl2-EXl2∣X∈Bvρ∣X∈Bvρ.On setting Zk=(Xk2-E[Xk2∣X∈Bv(ρ)])/2λk, we can represent Ω as Ω=E[ZZT∣X∈Bv(ρ)]. If x∈Rv is not the null vector, then xTΩx=E[xTZZTx∣X∈Bv(ρ)]=E[(xTZ)2∣X∈Bv(ρ)]>0. Moreover, the eigenvalues of Ω fulfill the secular equation(50)0=det⁡Ω-ϕIv=det⁡Λ-1Ω-ϕIvΛ=det⁡Λ-1ΩΛ-ϕIv=det⁡J-ϕIv,whence it follows that J is positive definite as well (though it is not symmetric). Since the sum of positive definite matrices is positive definite, we conclude that ∫01‍dtJρ;λ′′+tλ′-λ′′ is positive definite too. As such, it is nonsingular. Therefore, from (47), we conclude that λ′=λ′′.


## 5. Numerical Computation of Gaussian Integrals overBvρ

Let us now see how to computeαklm⋯ with controlled precision. Most of the relevant work has been originally done by Ruben in [2], where the case of α is discussed. We extend Ruben’s technique to Gaussian integrals containing powers of the integration variable. Specifically, it is shown in [2] that α(ρ;λ) can be represented as a series of chi-square cumulative distribution functions:(51)αρ;λ=∑m=0∞cms;λFv+2mρs.The scale factor s has the same physical dimension as ρ and λ. It is introduced in order to factorize the dependence of α upon ρ and λ at each order of the expansion. The series on theright-hand side of (51) converges uniformly on every finite interval of ρ. The coefficients cm are given by (52)cms;λ=1m!sv/2+mΛ1/2Γv/2+mΓv/2M-Qm,m=0,1,…,having defined Q(x)≡xT[Λ-1-s-1Iv]x for x∈Rv and M as the uniform average operator on the (v-1)-sphere ∂Bv(1)≡{u∈Rv:uTu=1}; namely,(53)Mϕ≡Γv/22πv/2∫∂Bv1duϕu,∀ϕ∈C0∂Bv1a.e.Unfortunately, (52) is not particularly convenient for numerical computations, since M[(-Q)m] is only given in integral form. However, it is also shown in [2] that the coefficients cm can be extracted from the Taylor expansion (at z0=0) of the generating function(54)ψz=∏k=1vsλk1/21-1-sλkz-1/2,that  isψz=∑m=0∞cms;λzm.This series converges uniformly for z<mini⁡1-s/λi-1. On evaluating the derivatives of ψ(z), it is then shown that cm’s fulfill the recursion:(55)c0=∏m=1vsλm;cn=12n∑r=0n-1gn-rcr,n=1,2,…;gn≡∑m=1v1-sλmn.Finally, the systematic error produced on considering only the lowest k terms of the chi-square series of (51) is estimated by(56)Rnρ;λ≡∑m=n∞cms;λFv+2mρs≤c0s;λΓv2+nΓv2ηnn!1-η-v2+nFv+2n1-ηρs≡Rn,with η=maxi⁡1-s/λi.Now, as mentioned, it is possible to extend the above expansion to all Gaussian integralsαklm⋯. Here, we are interested only in αk and αjk, since these are needed in order to implement the fixed point iteration and to compute the Jacobian matrix of τρ. The extension is provided by the following.Theorem 9 (Ruben’s expansions).
The integralsαk and αjk admit the series representations:(57)αkρ;λ=∑m=0∞ck;ms;λFv+2m+1ρs,(58)αjkρ;λ=∑m=0∞cjk;ms;λFv+2m+2ρs,with s being an arbitrary positive constant. The series coefficients are given, respectively, by (59)ck;ms;λ=sλkv+2mm!sv/2+mΛ1/2Γv/2+mΓv/2M-Qmuk2,(60)cjk;ms;λ=1+2δjksλjsλkv+2m+2v+2mm!sv/2+mΛ1/2Γv/2+mΓv/2M-Qmuj2uk2,with δjk denoting the Kronecker symbol. The series on theright-hand side of (57)-(58) converge uniformly on every finite interval of ρ. The functions(61)ψkz=sλk3/21-1-sλkz-3/2∏i≠ksλi1/21-1-sλiz-1/2,(62)ψkkz=3sλk5/21-1-sλkz-5/2∏i≠ksλi1/21-1-sλiz-1/2,(63)ψjkz=sλjsλk3/21-1-sλjz1-1-sλkz-3/2×∏i≠j,ksλi1/21-1-sλiz-1/2j≠kare generating functions, respectively, for the coefficients ck;m, ckk;m and cjk;m(j≠k); that is, they fulfill(64)ψkz=∑m=0∞ck;ms;λzm,ψjkz=∑m=0∞cjk;ms;λzm,for z<mini⁡1-s/λi-1. Finally, the coefficients ck;m,  ckk;m, and cjk;m(j≠k) can be obtained iteratively from the recursions(65)ck;0=sλkc0;ck;m=12m∑r=0m-1gk;m-rck;r;gk;m≡∑i=1vek;i1-sλim,m≥1;(66)cjk;0=1+2δjksλjsλkc0;cjk;m=12m∑r=0m-1gjk;m-rcjk;r;gjk;m≡∑i=1vejk;i1-sλim,m≥1,where the auxiliary coefficients ek;i and ejk;i are defined by(67)ek;i=3ifi=k1otherwise,(68)ekk;i=5ifi=k1otherwise,ejk;i=3ifi=jork1otherwisej≠k.It is not difficult to further generalize this theorem, so as to provide a chi-square expansion for any Gaussian integralαklm⋯. The proof follows closely the original one given by Ruben. We reproduce it in the appendix for αk, just to highlight the differences arising when the Gaussian integral contains powers of the integration variable.Analogously to (56), it is possible to estimate the systematic error produced when considering only the lowest k terms of the chi-square series of αk and αjk. Specifically, we find that(69)Rk;n≡∑m=n∞ck;ms;λFv+2m+1ρs≤ck;0ηnn!1-η-v/2+n+1Γv/2+n+1Γv/2Fv+2n+11-ηρs≡Rk;n,Rjk;n≡∑m=n∞cjk;ms;λFv+2m+2ρs≤cjk;0ηnn!1-η-v/2+n+2Γv/2+n+2Γv/2Fv+2n+21-ηρs≡Rjk;n.In order to evaluate all Ruben series with controlled uncertainty, we first set (see once more [2] for an exhaustive discussion on how to choose s)   s=2λ1λv/(λ1+λv); then we choose a unique threshold ε representing the maximum tolerable systematic error; for example, εdp=1.0×10-14 (roughly corresponding to double floating-point precision), for all α, αk, and αjk, and finally for each αX we compute the integer(70)kth≡minn≥1⁡n:RX;n<ε,providing the minimum number of chi-square terms, for which the upper bound RX;n to the residual sum RX;n lies below ε. Of course, this procedure overshoots the minimum number of terms really required for the R’s to lie below ε, since we actually operate on the R’s instead of the R’s. Nevertheless, the computational overhead is acceptable, as it will be shown in the next section. For the sake of completeness, it must be said that typically the values of kth for α, αk, and αjk with the same ϵ (and ρ, λ) are not much different from each other.To conclude, we notice thatkth depends nontrivially upon λ. By contrast, since Fv(x) is monotically increasing in x, we clearly see that kth is monotically increasing in ρ. Now, should one evaluate α and the like for a given λ at several values of ρ, say ρ1≤ρ2≤⋯≤ρmax, it is advisable to save computing resources and work out Ruben coefficients just once, up to the order kth corresponding to ρmax, since kth(ρ1)≤⋯≤kth(ρmax). We made use of this trick throughout our numerical experiences, as reported in the sequel.


## 6. Numerical Analysis of the Reconstruction Process

The fixed point (39) represents the simplest iterative scheme that can be used in order to reconstruct the solution λ=τρ-1·μ. In the literature of numerical methods, this scheme is known as a nonlinear Gauss–Jacobi (GJ) iteration (see, e.g., [20]). Accordingly, we shall rewrite it as λGJ,k(n+1)=TkλGJn. As we have seen, the sequence λGJ(n) converges with no exception as n→∞, provided μ∈D(τρ-1). Given ϵT>0, the number of steps nit needed for an approximate convergence with relative precision ϵT, that is,(71)nit≡minn≥1⁡n:λGJn-λGJn-1∞λGJn-1∞<ϵT,depends not only upon ϵT, but also on ρ and μ (note that the stopping rule is well conditioned, since λ(n)∞>0∀n and also limn→∞⁡λ(n)∞>0). In order to characterize statistically the convergence rate of the reconstruction process, we must integrate out the fluctuations of nit due to changes of μ; that is, we must average nit by letting μ fluctuate across its own probability space. In this way, we obtain the quantity n-it≡Eμnit∣ϵT,ρ, which better synthesizes the cost of the reconstruction for given ϵT and ρ. It should be evident that carrying out this idea analytically is hard, for on the one hand nit depends upon μ nonlinearly and on the other hand μ has a complicated distribution, as we briefly explain below.


### 6.1. Choice of the Eigenvalue Ensemble

Sinceλ is the eigenvalue spectrum of a full covariance matrix, it is reasonable to assume its distribution to be a Wishart Wv(p,Σ0) for some scale matrix Σ0 and for some number of degrees of freedom p≥v. In the sequel, we shall make the ideal assumption Σ0=p-1·Iv, so that the probability measure of λ is (see, e.g., [21])(72)dwvp;λ=pp+v2-1/2πv2/2∏k=1vλkp-v-1/2exp⁡-p/2∑k=1vλk∏k<jλj-λk2vp/2Γvp/2Γvv/2dvλ.Under this assumption, the probability measure of μ is obtained by performing the change of variable λ=τρ-1·μ in (72). Unfortunately, we have no analytic representation of τρ-1. Thus, we have neither an expression for the distribution of μ. However, μ can be extracted numerically as follows:(i)
generate randomlyΣ~Wv(p,p-1·Iv) by means of the Bartlett decomposition [22];(ii)
take the ordered eigenvalue spectrumλ of Σ;(iii)
obtainμ by applying the truncation operator τρ to λ.Note that since Wv(p,p-1·Iv) is only defined for p≥v, we need to rescale p as v increases. The simplest choice is to keep the ratio p/v fixed. The larger this ratio, the closer Σ fluctuates around Iv (recall that if Σ~Wv(p,p-1·Iv), then E[Σij]=δij and var⁡(Σij)=p-1[1+δij]). In view of this, large values of p/v are to be avoided, since they reduce the probability of testing the fixed point iteration on eigenvalue spectra characterized by large condition numbers ncond≡λv/λ1. For this reason, we have set p=2v in our numerical study.Having specified an ensemble of matrices from which the eigenvalue spectra are extracted, we are now ready to perform numerical simulations. To begin with, we report in Figure3 the marginal probability density function of the ordered eigenvalues λkk=1v and their truncated counterparts μkk=1v for the Wishart ensemble W10(20,20-1·I10) at ρ=1, as obtained numerically from a rather large sample of matrices (≃106 units). It will be noted that (i) the effect of the truncation is severe on the largest eigenvalues, as a consequence of the analytic bounds of Corollary 2.1 and Proposition 2.2; (ii) while the skewness of the lowest truncated eigenvalues is negative, it becomes positive for the largest ones. This is due to a change of relative effectiveness of (29) i with respect to (29) ii.Figure 3
Monte Carlo simulation of the probability density function of the ordered eigenvaluesλk (even rows) and their truncated counterparts μk at ρ=1 (odd rows) for the Wishart ensemble W10(20,20-1⋅I10). The last two plots (bottom right) display the distribution of the sum of eigenvalues.


### 6.2. Choice of the Simulation Parameters

In order to explore the dependence ofn-it upon ρ, we need to choose one or more simulation points for the latter. Ideally, it is possible to identify three different regimes in our problem: ρ≲λ1 (strong truncation regime), λ1≲ρ≲λv (crossover), and ρ≲λv (weak truncation regime). We cover all of them with the following set of points:(73)ρ∈Moλ1,…,Moλv∪12Moλ1,2Moλv,where Mo(·) stands for the mode. In principle, it is possible to determine Mo(λk) with high accuracy by using analytic representations of the marginal probability densities of the ordered eigenvalues [24]. In practice, the latter become computationally demanding at increasingly large values of v: for instance, the determination of the probability density of λ2 requires v!2 sums, which is unfeasible even at v~10. Moreover, to our aims, it is sufficient to choose approximate values, provided that these lie not far from the exact ones. Accordingly, we have determined the eigenvalue modes numerically from samples of N≃106 Wishart matrices. Our estimates are reported in Table 1 for v=3,…,10. They have been obtained from Grenander’s estimator [23]:(74)Moλkrs=12∑i=1N-rλki+λki+rλki-λki+r-s∑i=1N-rλki-λki+r-s,with properly chosen parameters r, s.Table 1
Numerical estimates of the mode of the ordered eigenvaluesλ1,…,λv of Σ∼Wv2v,2v-1⋅Iv with v=3,…,10. The estimates have been obtained from Grenander’s mode estimator [23].
v = 3 v = 4 v = 5 v = 6 v = 7 v = 8 v = 9 v = 10 M o ^ ( λ 1 ) 0.1568 0.1487 0.1435 0.1383 0.1344 0.1310 0.1269 0.1258 M o ^ ( λ 2 ) 0.6724 0.4921 0.4017 0.3424 0.3039 0.2745 0.2554 0.2399 M o ^ ( λ 3 ) 1.6671 1.0112 0.7528 0.6071 0.5138 0.4543 0.4048 0.3693 M o ^ ( λ 4 ) — 1.8507 1.2401 0.9621 0.7854 0.6684 0.5858 0.5288 M o ^ ( λ 5 ) — — 2.0150 1.4434 1.1269 0.9263 0.7956 0.7032 M o ^ ( λ 6 ) — — — 2.1356 1.5789 1.2559 1.0527 0.9111 M o ^ ( λ 7 ) — — — — 2.2190 1.6764 1.3673 1.1603 M o ^ ( λ 8 ) — — — — — 2.2763 1.7687 1.4624 M o ^ ( λ 9 ) — — — — — — 2.3210 1.8473 M o ^ ( λ 10 ) — — — — — — — 2.3775We are now in the position to investigate numerically how many terms in Ruben’s expansions must be considered asε is set to εdp=1.0×10-14, for our choice of the eigenvalue ensemble λ~Wv2v,(2v)-1·Iv and with ρ set as in Table 1. As an example, we report in Figure 4 the discrete distributions of kth for the basic Gaussian integral α at v=10, the largest dimension that we have simulated. As expected, we observe an increase of kth with ρ. Nevertheless, we see that the number of Ruben’s components to be taken into account for a double precision result keeps altogether modest even in the weak truncation regime, which proves the practical usefulness of the chi-square expansions.Figure 4
Monte Carlo simulation of the probability mass function of the parameterkth for the Gaussian probability content α. The histograms refer to the eigenvalue ensemble λ of Σ∼W10(20,20-1⋅I10), with ρ chosen as in Table 1 and ε=1.0×10-14.


### 6.3. Fixed Point Iteration at Work

The GJ iteration is too slow to be of practical interest. For instance, atv=10, ρ≃Mo(λ1) and ϵT=1.0×10-7 (corresponding to a reconstruction of λ with single floating-point precision), it is rather easy to extract realizations of μ which require nit≃15,000 to converge. An improvement of the GJ scheme is achieved via overrelaxation (GJOR); that is,(75)λGJOR,k0=μk,λGJOR,kn+1=λGJOR,kn+ωTkλGJORn-λGJOR,kn,k=1,…,v.Evidently, at ω=1, the GJOR scheme coincides with the standard GJ one. The optimal value ωopt of the relaxation factor ω is not obvious even in the linear Jacobi scheme, where ωopt depends upon the properties of the coefficient matrix of the system. For instance, if the latter is symmetric positive definite, it is demonstrated that the best choice is provided by ωopt≡2(1+1-σ2)-1, σ being the spectral radius of the Jacobi iteration matrix [25]. In our numerical tests with the GJOR scheme, we found empirically that the optimal value of ω at ρ≲λv is close to the linear prediction, provided that σ is replaced by J∞, J being defined as in Section 3 (note that J∞<1).By contrast, the iteration diverges after few steps with increasing probability asρ/λv→∞ if ω is kept fixed at ω=ωopt; in order to restore the convergence, ω must be lowered towards ω=1 as such limit is taken.To give an idea of the convergence rate of the GJOR scheme, we show in Figure5(a) a joint box-plot of the distributions of nit at v=10 and ϵT=1.0×10-7. From the plot we observe that the distribution of nit shifts rightwards as ρ decreases: clearly, the reconstruction is faster if ρ is in the weak truncation regime (where μ is closer to λ), whereas it takes more iterations in the strong truncation regime. The dependence of n-it upon ρ, systematically displayed in Figure 6, is compatible with a scaling law(76)log⁡n-itρ,v,ϵT=av,ϵT-bv,ϵTlog⁡ρ,apart from small corrections occurring at large ρ. Equation (76) tells us that n-it increases polynomially in 1/ρ at fixed v. In order to estimate the parameters a and b in the strong truncation regime (where the algorithm becomes challenging), we performed jackknife fits to (76) of data points with ρ≲1. Results are collected in Figure 5(b), showing that b is roughly constant, while a increases almost linearly in v. Thus, while the cost of the eigenvalue reconstruction is only polynomial in 1/ρ at fixed v, it is exponential in v at fixed ρ. The scaling law of the GJOR scheme is therefore better represented by n-it=Ceκv/ρb, with C being a normalization constant independent of ρ and v and κ representing approximately the slope of a as a function of v. Although the GJOR scheme improves the GJ one, the iteration reveals to be still inefficient in a parameter subspace, which is critical for the applications.Figure 5
(a) Box-plot ofnit in the GJOR scheme at v=10, with ϵT=1.0×10-7 and ρ chosen as in Table 1. The distributions have been reconstructed from a sample of N≃103 eigenvalue spectra extracted from W10(20,20-1⋅I10). The whiskers extend to the most extreme data point within (3/2)(75%−25%) data range. (b) Numerical estimates of the scaling parameters a and b of the GJOR scheme, as obtained from jackknife fits to (76) of data points with ρ≲1 and ϵT=1.0×10-7. We quote in parentheses the jackknife error.
(a) (b)Figure 6
Log-log plots ofn-it versus ρ in the GJOR scheme at ϵT=1.0×10-7. The parameter ρ has been chosen as in Table 1. The (red) dashed line in each plot represents our best jackknife linear fit to (76) of data points with ρ≲1.


### 6.4. Boosting the GJOR Scheme

A further improvement can be obtained by lettingω depend on the eigenvalue index in the GJOR scheme. Let us discuss how to work out such an adjustment. On commenting on Figure 3, we have already noticed that the largest eigenvalues are affected by the truncation to a larger extent than the smallest ones. Therefore, they must perform a longer run through the fixed point iteration, in order to converge to the untruncated values. This is a possible qualitative explanation for the slowing down of the algorithm as ρ→0. In view of it, we expect to observe some acceleration of the convergence rate, if ω is replaced, for instance, by(77)ω⟶ωk≡1+β·kωopt,β≥0,k=1,…,v.The choice β=0 corresponds obviously to the standard GJOR scheme. Any other choice yields ωk>ωopt. Therefore, the new scheme is also expected to display a higher rate of failures than the GJOR one at ρ≫λv, for the reason explained in Section 5.3. The componentwise overrelaxation proposed in (77) is only meant to enhance the convergence speed in the strong truncation regime and in the crossover, where the improvement is actually needed.In order to confirm this picture, we have explored systematically the effect ofβ on n-it by simulating the reconstruction process at v=3,…,10, with β varying from 0 to 2 in steps of 1/5. First of all, we have observed that the rate of failures at large ρ is fairly reduced if the first 30÷50 iterations are run with ωk=ωopt, and only afterwards β is switched on. Having minimized the failures, we have checked that for each value of β, the scaling law assumed in (76) is effectively fulfilled. Then, we have computed jackknife estimates of the scaling parameters a and b. These are plotted in Figure 7 as functions of v. Each trajectory (represented by a dashed curve) corresponds to a given value of β. Those with darker markers refer to smaller values of β and the other way round. From the plots we notice that(i)
all the trajectories withβ>0 lie below the one with β=0;(ii)
the trajectories ofa display a clear increasing trend with v; yet their slope lessens as β increases. By contrast, the trajectories of b develop a mild increasing trend with v as β increases, though this is not strictly monotonic;(iii)
the trajectories of botha and b seem to converge to a limit trajectory as β increases; we observe a saturation phenomenon, which thwarts the benefit of increasing β beyond a certain threshold close to βmax⋍2.We add that pushing β beyond βmax is counterproductive, as the rate of failures becomes increasingly relevant in the crossover and eventually also in the strong truncation regime. By contrast, if β≲βmax the rate of failures keeps very low for essentially all simulated values of ρ.Figure 7
Scaling parametersa and b of the modified GJOR scheme as functions of v at ϵT=1.0×10-7, with β varying in the range 0÷2 in steps of 1/5. Each trajectory (represented by a dashed curve) refers to a different value of β. Those with darker markers correspond to smaller values of β and the other way round.Our numerical results signal a strong reduction of the slowing down of the convergence rate. Indeed, (i) means qualitatively thatC and b are reduced as β increases. (ii) means that κ is reduced as β increases (this is the most important effect, as κ is mainly responsible for the exponential slowing down with v). The appearance of a slope in the trajectories of b as β increases indicates that a mild exponential slowing down is also developed at denominator of the scaling law n-it=Ceκv/ρb, but the value of b is anyway smaller than at β=0. Finally, (iii) means that choosing β>βmax has a minor impact on the performance of the algorithm. In Figure 8, we report a plot of the parameter κ (obtained from least-squares fits of data to a linear model a=a0+κ·v) as a function of β. We see that κ(β=0)/κ(β=2)≃4. This quantifies the maximum exponential speedup of the convergence rate, which can be achieved by our proposal. When β is close to βmax, n-it amounts to few hundreds at v=10 and ρ≃λ1/2.Figure 8
The parameterκ as a function of β. Estimates of κ are obtained from least-squares fits of data to a linear model a=a0+κ⋅v.


## 6.1. Choice of the Eigenvalue Ensemble

Sinceλ is the eigenvalue spectrum of a full covariance matrix, it is reasonable to assume its distribution to be a Wishart Wv(p,Σ0) for some scale matrix Σ0 and for some number of degrees of freedom p≥v. In the sequel, we shall make the ideal assumption Σ0=p-1·Iv, so that the probability measure of λ is (see, e.g., [21])(72)dwvp;λ=pp+v2-1/2πv2/2∏k=1vλkp-v-1/2exp⁡-p/2∑k=1vλk∏k<jλj-λk2vp/2Γvp/2Γvv/2dvλ.Under this assumption, the probability measure of μ is obtained by performing the change of variable λ=τρ-1·μ in (72). Unfortunately, we have no analytic representation of τρ-1. Thus, we have neither an expression for the distribution of μ. However, μ can be extracted numerically as follows:(i)
generate randomlyΣ~Wv(p,p-1·Iv) by means of the Bartlett decomposition [22];(ii)
take the ordered eigenvalue spectrumλ of Σ;(iii)
obtainμ by applying the truncation operator τρ to λ.Note that since Wv(p,p-1·Iv) is only defined for p≥v, we need to rescale p as v increases. The simplest choice is to keep the ratio p/v fixed. The larger this ratio, the closer Σ fluctuates around Iv (recall that if Σ~Wv(p,p-1·Iv), then E[Σij]=δij and var⁡(Σij)=p-1[1+δij]). In view of this, large values of p/v are to be avoided, since they reduce the probability of testing the fixed point iteration on eigenvalue spectra characterized by large condition numbers ncond≡λv/λ1. For this reason, we have set p=2v in our numerical study.Having specified an ensemble of matrices from which the eigenvalue spectra are extracted, we are now ready to perform numerical simulations. To begin with, we report in Figure3 the marginal probability density function of the ordered eigenvalues λkk=1v and their truncated counterparts μkk=1v for the Wishart ensemble W10(20,20-1·I10) at ρ=1, as obtained numerically from a rather large sample of matrices (≃106 units). It will be noted that (i) the effect of the truncation is severe on the largest eigenvalues, as a consequence of the analytic bounds of Corollary 2.1 and Proposition 2.2; (ii) while the skewness of the lowest truncated eigenvalues is negative, it becomes positive for the largest ones. This is due to a change of relative effectiveness of (29) i with respect to (29) ii.Figure 3
Monte Carlo simulation of the probability density function of the ordered eigenvaluesλk (even rows) and their truncated counterparts μk at ρ=1 (odd rows) for the Wishart ensemble W10(20,20-1⋅I10). The last two plots (bottom right) display the distribution of the sum of eigenvalues.


## 6.2. Choice of the Simulation Parameters

In order to explore the dependence ofn-it upon ρ, we need to choose one or more simulation points for the latter. Ideally, it is possible to identify three different regimes in our problem: ρ≲λ1 (strong truncation regime), λ1≲ρ≲λv (crossover), and ρ≲λv (weak truncation regime). We cover all of them with the following set of points:(73)ρ∈Moλ1,…,Moλv∪12Moλ1,2Moλv,where Mo(·) stands for the mode. In principle, it is possible to determine Mo(λk) with high accuracy by using analytic representations of the marginal probability densities of the ordered eigenvalues [24]. In practice, the latter become computationally demanding at increasingly large values of v: for instance, the determination of the probability density of λ2 requires v!2 sums, which is unfeasible even at v~10. Moreover, to our aims, it is sufficient to choose approximate values, provided that these lie not far from the exact ones. Accordingly, we have determined the eigenvalue modes numerically from samples of N≃106 Wishart matrices. Our estimates are reported in Table 1 for v=3,…,10. They have been obtained from Grenander’s estimator [23]:(74)Moλkrs=12∑i=1N-rλki+λki+rλki-λki+r-s∑i=1N-rλki-λki+r-s,with properly chosen parameters r, s.Table 1
Numerical estimates of the mode of the ordered eigenvaluesλ1,…,λv of Σ∼Wv2v,2v-1⋅Iv with v=3,…,10. The estimates have been obtained from Grenander’s mode estimator [23].
v = 3 v = 4 v = 5 v = 6 v = 7 v = 8 v = 9 v = 10 M o ^ ( λ 1 ) 0.1568 0.1487 0.1435 0.1383 0.1344 0.1310 0.1269 0.1258 M o ^ ( λ 2 ) 0.6724 0.4921 0.4017 0.3424 0.3039 0.2745 0.2554 0.2399 M o ^ ( λ 3 ) 1.6671 1.0112 0.7528 0.6071 0.5138 0.4543 0.4048 0.3693 M o ^ ( λ 4 ) — 1.8507 1.2401 0.9621 0.7854 0.6684 0.5858 0.5288 M o ^ ( λ 5 ) — — 2.0150 1.4434 1.1269 0.9263 0.7956 0.7032 M o ^ ( λ 6 ) — — — 2.1356 1.5789 1.2559 1.0527 0.9111 M o ^ ( λ 7 ) — — — — 2.2190 1.6764 1.3673 1.1603 M o ^ ( λ 8 ) — — — — — 2.2763 1.7687 1.4624 M o ^ ( λ 9 ) — — — — — — 2.3210 1.8473 M o ^ ( λ 10 ) — — — — — — — 2.3775We are now in the position to investigate numerically how many terms in Ruben’s expansions must be considered asε is set to εdp=1.0×10-14, for our choice of the eigenvalue ensemble λ~Wv2v,(2v)-1·Iv and with ρ set as in Table 1. As an example, we report in Figure 4 the discrete distributions of kth for the basic Gaussian integral α at v=10, the largest dimension that we have simulated. As expected, we observe an increase of kth with ρ. Nevertheless, we see that the number of Ruben’s components to be taken into account for a double precision result keeps altogether modest even in the weak truncation regime, which proves the practical usefulness of the chi-square expansions.Figure 4
Monte Carlo simulation of the probability mass function of the parameterkth for the Gaussian probability content α. The histograms refer to the eigenvalue ensemble λ of Σ∼W10(20,20-1⋅I10), with ρ chosen as in Table 1 and ε=1.0×10-14.


## 6.3. Fixed Point Iteration at Work

The GJ iteration is too slow to be of practical interest. For instance, atv=10, ρ≃Mo(λ1) and ϵT=1.0×10-7 (corresponding to a reconstruction of λ with single floating-point precision), it is rather easy to extract realizations of μ which require nit≃15,000 to converge. An improvement of the GJ scheme is achieved via overrelaxation (GJOR); that is,(75)λGJOR,k0=μk,λGJOR,kn+1=λGJOR,kn+ωTkλGJORn-λGJOR,kn,k=1,…,v.Evidently, at ω=1, the GJOR scheme coincides with the standard GJ one. The optimal value ωopt of the relaxation factor ω is not obvious even in the linear Jacobi scheme, where ωopt depends upon the properties of the coefficient matrix of the system. For instance, if the latter is symmetric positive definite, it is demonstrated that the best choice is provided by ωopt≡2(1+1-σ2)-1, σ being the spectral radius of the Jacobi iteration matrix [25]. In our numerical tests with the GJOR scheme, we found empirically that the optimal value of ω at ρ≲λv is close to the linear prediction, provided that σ is replaced by J∞, J being defined as in Section 3 (note that J∞<1).By contrast, the iteration diverges after few steps with increasing probability asρ/λv→∞ if ω is kept fixed at ω=ωopt; in order to restore the convergence, ω must be lowered towards ω=1 as such limit is taken.To give an idea of the convergence rate of the GJOR scheme, we show in Figure5(a) a joint box-plot of the distributions of nit at v=10 and ϵT=1.0×10-7. From the plot we observe that the distribution of nit shifts rightwards as ρ decreases: clearly, the reconstruction is faster if ρ is in the weak truncation regime (where μ is closer to λ), whereas it takes more iterations in the strong truncation regime. The dependence of n-it upon ρ, systematically displayed in Figure 6, is compatible with a scaling law(76)log⁡n-itρ,v,ϵT=av,ϵT-bv,ϵTlog⁡ρ,apart from small corrections occurring at large ρ. Equation (76) tells us that n-it increases polynomially in 1/ρ at fixed v. In order to estimate the parameters a and b in the strong truncation regime (where the algorithm becomes challenging), we performed jackknife fits to (76) of data points with ρ≲1. Results are collected in Figure 5(b), showing that b is roughly constant, while a increases almost linearly in v. Thus, while the cost of the eigenvalue reconstruction is only polynomial in 1/ρ at fixed v, it is exponential in v at fixed ρ. The scaling law of the GJOR scheme is therefore better represented by n-it=Ceκv/ρb, with C being a normalization constant independent of ρ and v and κ representing approximately the slope of a as a function of v. Although the GJOR scheme improves the GJ one, the iteration reveals to be still inefficient in a parameter subspace, which is critical for the applications.Figure 5
(a) Box-plot ofnit in the GJOR scheme at v=10, with ϵT=1.0×10-7 and ρ chosen as in Table 1. The distributions have been reconstructed from a sample of N≃103 eigenvalue spectra extracted from W10(20,20-1⋅I10). The whiskers extend to the most extreme data point within (3/2)(75%−25%) data range. (b) Numerical estimates of the scaling parameters a and b of the GJOR scheme, as obtained from jackknife fits to (76) of data points with ρ≲1 and ϵT=1.0×10-7. We quote in parentheses the jackknife error.
(a) (b)Figure 6
Log-log plots ofn-it versus ρ in the GJOR scheme at ϵT=1.0×10-7. The parameter ρ has been chosen as in Table 1. The (red) dashed line in each plot represents our best jackknife linear fit to (76) of data points with ρ≲1.


## 6.4. Boosting the GJOR Scheme

A further improvement can be obtained by lettingω depend on the eigenvalue index in the GJOR scheme. Let us discuss how to work out such an adjustment. On commenting on Figure 3, we have already noticed that the largest eigenvalues are affected by the truncation to a larger extent than the smallest ones. Therefore, they must perform a longer run through the fixed point iteration, in order to converge to the untruncated values. This is a possible qualitative explanation for the slowing down of the algorithm as ρ→0. In view of it, we expect to observe some acceleration of the convergence rate, if ω is replaced, for instance, by(77)ω⟶ωk≡1+β·kωopt,β≥0,k=1,…,v.The choice β=0 corresponds obviously to the standard GJOR scheme. Any other choice yields ωk>ωopt. Therefore, the new scheme is also expected to display a higher rate of failures than the GJOR one at ρ≫λv, for the reason explained in Section 5.3. The componentwise overrelaxation proposed in (77) is only meant to enhance the convergence speed in the strong truncation regime and in the crossover, where the improvement is actually needed.In order to confirm this picture, we have explored systematically the effect ofβ on n-it by simulating the reconstruction process at v=3,…,10, with β varying from 0 to 2 in steps of 1/5. First of all, we have observed that the rate of failures at large ρ is fairly reduced if the first 30÷50 iterations are run with ωk=ωopt, and only afterwards β is switched on. Having minimized the failures, we have checked that for each value of β, the scaling law assumed in (76) is effectively fulfilled. Then, we have computed jackknife estimates of the scaling parameters a and b. These are plotted in Figure 7 as functions of v. Each trajectory (represented by a dashed curve) corresponds to a given value of β. Those with darker markers refer to smaller values of β and the other way round. From the plots we notice that(i)
all the trajectories withβ>0 lie below the one with β=0;(ii)
the trajectories ofa display a clear increasing trend with v; yet their slope lessens as β increases. By contrast, the trajectories of b develop a mild increasing trend with v as β increases, though this is not strictly monotonic;(iii)
the trajectories of botha and b seem to converge to a limit trajectory as β increases; we observe a saturation phenomenon, which thwarts the benefit of increasing β beyond a certain threshold close to βmax⋍2.We add that pushing β beyond βmax is counterproductive, as the rate of failures becomes increasingly relevant in the crossover and eventually also in the strong truncation regime. By contrast, if β≲βmax the rate of failures keeps very low for essentially all simulated values of ρ.Figure 7
Scaling parametersa and b of the modified GJOR scheme as functions of v at ϵT=1.0×10-7, with β varying in the range 0÷2 in steps of 1/5. Each trajectory (represented by a dashed curve) refers to a different value of β. Those with darker markers correspond to smaller values of β and the other way round.Our numerical results signal a strong reduction of the slowing down of the convergence rate. Indeed, (i) means qualitatively thatC and b are reduced as β increases. (ii) means that κ is reduced as β increases (this is the most important effect, as κ is mainly responsible for the exponential slowing down with v). The appearance of a slope in the trajectories of b as β increases indicates that a mild exponential slowing down is also developed at denominator of the scaling law n-it=Ceκv/ρb, but the value of b is anyway smaller than at β=0. Finally, (iii) means that choosing β>βmax has a minor impact on the performance of the algorithm. In Figure 8, we report a plot of the parameter κ (obtained from least-squares fits of data to a linear model a=a0+κ·v) as a function of β. We see that κ(β=0)/κ(β=2)≃4. This quantifies the maximum exponential speedup of the convergence rate, which can be achieved by our proposal. When β is close to βmax, n-it amounts to few hundreds at v=10 and ρ≃λ1/2.Figure 8
The parameterκ as a function of β. Estimates of κ are obtained from least-squares fits of data to a linear model a=a0+κ⋅v.


## 7. On the Ill-Posedness of the Reconstruction in Sample Space

So far we have discussed the covariance reconstruction under the assumption thatμ=τρ·λ represents the exact truncated counterpart of some λ∈Rv and we have looked at the algorithmic properties of the iteration schemes which operatively define τρ-1. Such analysis is essential in order to characterize τρ-1 mathematically; yet it is not sufficient in real situations, specifically when μ is perturbed by statistical noise.In this section, we examine the difficulties arising when performing the covariance reconstruction in sample space. We first recall that, according to Hadamard [26], a mathematical problem is well posed provided that the following conditions are fulfilled:(H1)
there exists always a solution to the problem;(H2)
the solution is unique;(H3)
the solution depends smoothly on the input data.Inverse problems are often characterized by violation of one or more of them; see, for instance, [13]. In such cases, the standard practice consists in regularizing the inverse operator, that is, in replacing it by a stable approximation. With regard to our problem, the reader will recognize that (H1) is violated (and the problem becomes ill-posed) as soon as the space of the input data is allowed to be a superset of D(τρ-1): once clarifying how μ is concretely estimated in the applications (Sections 7.1 and 7.2), we propose a perturbative regularization of τρ-1, which improves effectively the fulfillment of (H1) (Section 7.3). By contrast, Proposition 8 guarantees that whenever a solution exists, it is also unique; thus (H2) is never of concern. Finally, the fulfillment of (H3) depends on how the statistical noise on μ is nonlinearly inflated by the action of τρ-1. For the sake of conciseness, in the present paper, we just sketch the main ideas underlying the perturbative regularization of τρ-1, whereas a technical implementation of it and a discussion of (H3) are deferred to a separate paper [27].


### 7.1. Definition of the Sample Truncated Covariance Matrix

The examples of Section2 assume that (i) spherical truncations are operated on a representative sample PN=x(k)k=1N of X~Nv(0,Σ) with finite size N, (ii) ρ is known exactly, and (iii) the input budget for the covariance reconstruction is given by the subset(78)QM=x∈PN:x2<ρ,withQM=M≤N.As usual in the analysis of stochastic variables in sample space, we assume that the observations x(k) are realizations of i.i.d. stochastic variables X(k)~Nv(0,Σ), k=1,…,N. Thus, M is itself a stochastic variable in sample space, where it reads(79)M=∑k=1NIBvρXk≡∑k=1NIk,with IBvρ· denoting the characteristic function of Bv(ρ) and Ik≡IBvρXk being just a shortcut for its extended counterpart. It is easily recognized that M~B(N,α) is a binomial variate. If we indeed denote by E the sample expectation operator (i.e.,the integral with respect to the product measure of the joint variables X(k)k=1N), then a standard calculation yields(80)EM=E∑k=1NIk=∑k=1NEIk=∑k=1Nα=αN,varM=EM2-EM2=∑k=1NEIk2+∑k,s:k≠s1⋯NEIkIs-α2N2=∑k=1NEIk+∑k,s:k≠s1⋯NEIkEIs-α2N2=αN+α2NN-1-α2N2=α1-αN.Hence, we see that the relative dispersion of M is O(N-1/2). Now, the simplest way to measure Σ and SB, respectively, from the sets PN and QM is via the classical estimators(81)Σ^ij=1N-1∑x∈PNx-x-ix-x-j,x-i=1N∑x∈PNxi,(82)S^Bij=1M-1∑x∈QMx-x~ix-x~j,x~i=1M∑x∈QMxi.We define the sample estimates λ^ and μ^, respectively, of λ and μ as the eigenvalue spectra of Σ^ and S^B. By symmetry arguments we see that x~i is unbiased. Indeed, it holds(83)Ex~i=∑k=1NEXikIk∑s=1NIs.Theright-hand side of (83) makes only sense if we conventionally define the integrand to be zero in the integration subdomain {X(k)∉Bv(ρ),∀k} or equivalently if we interpret E[x~i] as the conditional one E[x~i∣M>0] (the event M>0 occurs a.s. only as N→∞). Since the sample measure is even under X(k)→-X(k) while the integrand is odd, we immediately conclude that bias[x~i]=0.


### 7.2. Bias of the Sample Truncated Covariance Matrix

The situation gets somewhat less trivial withS^B: the normalization factor M-1-1, which has been chosen in analogy with (81), is not sufficient to remove completely the bias of S^B at finite N, though we aim at showing here that the residual bias is exponentially small and asymptotically vanishing. In order to see this, we observe that(84)ES^Bij=∑k=1NEXik-X~iXjk-X~jIk∑s=1NIs-1=∑l,r=1vRilRjr∑k=1NEdiagXlk-X~lXrk-X~rIk∑s=1NIs-1,with Ediag denoting the sample expectation corresponding to a multinormal measure with diagonal covariance matrix Λ=diag⁡(λ)=RTΣR, conditioned to M>1. Having diagonalized the product measure, we observe that the integrand on theright-hand side is odd for l≠r and even for l=r under the joint change of variables Xl(k)→-Xl(k) for k=1,…,N, similarly to what we did in Section 2. As a consequence, it holds(85)ES^Bij=∑l=1vRilRjl∑k=1NEdiagXlk-X~l2Ik∑s=1NIs-1,whence we infer that the matrix E[S^B] is diagonalized by the same matrix R as Σ. From (85) we also conclude that(86)bias⁡S^B=Rdiag⁡wRT,wi≡∑k=1NEdiagXik-X~i2Ik∑s=1NIs-1-μi,i=1,…,v. It should be observed that in general wi≠bias⁡[μ^i] since the computation of μ^i requires the diagonalization of S^B, which is in general performed by a diagonalizing matrix R^≠R. Nevertheless, if w vanishes then bias[S^B] vanishes too. Now, we observe that wi splits into three contributions:(87)wi1=∑k=1NEdiagXik2Ik∑s=1NIs-1-μi,wi2=-2∑k=1NEdiagXikX~iIk∑s=1NIs-1,wi3=∑k=1NEdiagX~i2Ik∑s=1NIs-1,which can be exactly calculated and expressed in terms of μi, α, and N. For instance,(88)wi1=NEdiag1M-1Xi12I1∣M>1-μi=N∑m=2N1m-1EdiagXi12I1∣M=m-μi=N∑m=2N1m-1αμiN-1m-1αm-11-αN-m-μi=μi∑m=2Nmm-1Nmαm1-αN-m-μi.Analogously, we have(89)wi2=-2μi∑m=2N1m-1Nmαm1-αN-m,wi3=μi∑m=2N1m-1Nmαm1-αN-m.Hence, it follows that(90)wi=-μi1+αN-11-αN-1.Since α>0, we see that limN→∞wi=0. Thus, we conclude that S^B is asymptotically unbiased.A discussion of the variance of the sample truncated covariance matrix is beyond the scope of the present paper. We just observe that, apart from the above calculation, studying the sample properties of the truncated spectrum is made hard by the fact that eigenvalues and eigenvectors of a diagonalizable matrix are intimately related from their very definition; thus such study would require a careful consideration of the distribution of the sample diagonalizing matrixR^≠R of S^B.


### 7.3. Perturbative Regularization ofτρ-1

Whenμ is critically close to the internal boundary of D(τρ-1), a sample estimate μ^ may fall outside of it due to statistical fluctuations. In that case the iterative procedure described in the previous sections diverges. On the quantitative side, the ill-posedness of the reconstruction problem is measured by the failure probability(91)pfailρ,Σ,N=Pμ^∉Dτρ-1∣Xk~Nv0;Σ,k=1,…,N,which is a highly nontrivial function of ρ, Σ, and N. An illustrative example of it is reported in Figure 9(a), which refers to a specific case with v=4 and Σ=diag(0.1,0.3,0.8,2.2). The plot suggests that the iterative procedure becomes severely ill-posed in the regime of strong truncation.Figure 9
(a) Numerical reconstruction of the failure probability of the iterative procedure asv=4 and Σ=diag⁡(0.1,0.3,0.8,2.2), for several values of ρ and for N=200,250,…,1000. (b) Failure probability of the perturbative regularization with the same parameters.
(a) (b)In order to regularize the problem, we propose to go back to (19) and consider it from a different perspective. Specifically, we move from the observation that a simplified framework occurs in the special circumstance when the eigenvalue spectra are fully degenerate, which is essentially equivalent to the setup of [1]. If μ1=⋯=μv≡μ~, by symmetry arguments it follows that λ1=⋯=λv≡λ~ and the other way round. Equation (19) reduces in this limit to(92)μ~=λ~Fv+2Fvρλ~≡Tρλ~.It can be easily checked that the function Tρ(λ~) is monotonically increasing in λ~. In addition, we have(93)ilimλ~→0⁡Tρλ~=0,iilimλ~→∞⁡Tρλ~=ρv+2;thus (92) can be surely (numerically) inverted provided that 0<μ~<ρ/(v+2). We can regard (92) as an approximation to the original problem (19). When μ is not degenerate, we must define μ~ in terms of the components of μ. One possibility is to average them, that is, to choose(94)μ~=1v∑i=1vμi.Subject to this, we expect λ~ to lie somewhere between λ1 and λv. Equation (92) can be thought of as the lowest order approximation of a perturbative expansion of (19) around the point λT={λ~,…,λ~}. If the condition number of Σ is not extremely large, such an expansion is expected to quickly converge, so that a few perturbative corrections to λT should be sufficient to guarantee a good level of approximation.As mentioned above, a technical implementation of the perturbative approach and a thorough discussion of its properties are deferred to a separate paper [27]. Here, we limit ourselves to observing that the definition domain of perturbation theory is ultimately set by its lowest order approximation, since corrections to (92) are all algebraically built in terms of it, with no additional constraints. Following (94), the domain of Tρ-1 comes to be defined as(95)DTρ-1=μ∈R+v:∑i=1vμi≤ρvv+2,and it is clear that D(τρ-1)⊂D(Tρ-1) (it is sufficient to sum term by term all the inequalities contributing to (36)). In Figure 10, we show the set difference DTρ-1∖D(τρ-1) in v=2 and v=3 dimensions. When μ∈D(τρ-1) but its estimate μ^∉D(τρ-1), it may well occur μ^∈D(Tρ-1); that is, the set difference acts as an absorbing shield of the statistical noise. Therefore, if we define the failure probability of the perturbative reconstruction as(96)qfailρ,Σ,N=Pμ^∉DTρ-1∣Xk~Nv0;Σ,k=1,…,N,we expect the inequality qfail(ρ,Σ,N)≪pfail(ρ,Σ,N) to generously hold. An example is given in Figure 9(b): we see that qfail becomes lower than pfail by orders of magnitude as soon as ρ and N are not exceedingly small. In this sense, the operator Tρ-1 can be regarded as the lowest order approximation of a regularizing operator for τρ-1.Figure 10
(a) Set differenceDTρ-1∖Dτρ-1 in v=2 dimensions. (b) Set difference DTρ-1∖Dτρ-1 in v=3 dimensions.
(a) (b)


## 7.1. Definition of the Sample Truncated Covariance Matrix

The examples of Section2 assume that (i) spherical truncations are operated on a representative sample PN=x(k)k=1N of X~Nv(0,Σ) with finite size N, (ii) ρ is known exactly, and (iii) the input budget for the covariance reconstruction is given by the subset(78)QM=x∈PN:x2<ρ,withQM=M≤N.As usual in the analysis of stochastic variables in sample space, we assume that the observations x(k) are realizations of i.i.d. stochastic variables X(k)~Nv(0,Σ), k=1,…,N. Thus, M is itself a stochastic variable in sample space, where it reads(79)M=∑k=1NIBvρXk≡∑k=1NIk,with IBvρ· denoting the characteristic function of Bv(ρ) and Ik≡IBvρXk being just a shortcut for its extended counterpart. It is easily recognized that M~B(N,α) is a binomial variate. If we indeed denote by E the sample expectation operator (i.e.,the integral with respect to the product measure of the joint variables X(k)k=1N), then a standard calculation yields(80)EM=E∑k=1NIk=∑k=1NEIk=∑k=1Nα=αN,varM=EM2-EM2=∑k=1NEIk2+∑k,s:k≠s1⋯NEIkIs-α2N2=∑k=1NEIk+∑k,s:k≠s1⋯NEIkEIs-α2N2=αN+α2NN-1-α2N2=α1-αN.Hence, we see that the relative dispersion of M is O(N-1/2). Now, the simplest way to measure Σ and SB, respectively, from the sets PN and QM is via the classical estimators(81)Σ^ij=1N-1∑x∈PNx-x-ix-x-j,x-i=1N∑x∈PNxi,(82)S^Bij=1M-1∑x∈QMx-x~ix-x~j,x~i=1M∑x∈QMxi.We define the sample estimates λ^ and μ^, respectively, of λ and μ as the eigenvalue spectra of Σ^ and S^B. By symmetry arguments we see that x~i is unbiased. Indeed, it holds(83)Ex~i=∑k=1NEXikIk∑s=1NIs.Theright-hand side of (83) makes only sense if we conventionally define the integrand to be zero in the integration subdomain {X(k)∉Bv(ρ),∀k} or equivalently if we interpret E[x~i] as the conditional one E[x~i∣M>0] (the event M>0 occurs a.s. only as N→∞). Since the sample measure is even under X(k)→-X(k) while the integrand is odd, we immediately conclude that bias[x~i]=0.


## 7.2. Bias of the Sample Truncated Covariance Matrix

The situation gets somewhat less trivial withS^B: the normalization factor M-1-1, which has been chosen in analogy with (81), is not sufficient to remove completely the bias of S^B at finite N, though we aim at showing here that the residual bias is exponentially small and asymptotically vanishing. In order to see this, we observe that(84)ES^Bij=∑k=1NEXik-X~iXjk-X~jIk∑s=1NIs-1=∑l,r=1vRilRjr∑k=1NEdiagXlk-X~lXrk-X~rIk∑s=1NIs-1,with Ediag denoting the sample expectation corresponding to a multinormal measure with diagonal covariance matrix Λ=diag⁡(λ)=RTΣR, conditioned to M>1. Having diagonalized the product measure, we observe that the integrand on theright-hand side is odd for l≠r and even for l=r under the joint change of variables Xl(k)→-Xl(k) for k=1,…,N, similarly to what we did in Section 2. As a consequence, it holds(85)ES^Bij=∑l=1vRilRjl∑k=1NEdiagXlk-X~l2Ik∑s=1NIs-1,whence we infer that the matrix E[S^B] is diagonalized by the same matrix R as Σ. From (85) we also conclude that(86)bias⁡S^B=Rdiag⁡wRT,wi≡∑k=1NEdiagXik-X~i2Ik∑s=1NIs-1-μi,i=1,…,v. It should be observed that in general wi≠bias⁡[μ^i] since the computation of μ^i requires the diagonalization of S^B, which is in general performed by a diagonalizing matrix R^≠R. Nevertheless, if w vanishes then bias[S^B] vanishes too. Now, we observe that wi splits into three contributions:(87)wi1=∑k=1NEdiagXik2Ik∑s=1NIs-1-μi,wi2=-2∑k=1NEdiagXikX~iIk∑s=1NIs-1,wi3=∑k=1NEdiagX~i2Ik∑s=1NIs-1,which can be exactly calculated and expressed in terms of μi, α, and N. For instance,(88)wi1=NEdiag1M-1Xi12I1∣M>1-μi=N∑m=2N1m-1EdiagXi12I1∣M=m-μi=N∑m=2N1m-1αμiN-1m-1αm-11-αN-m-μi=μi∑m=2Nmm-1Nmαm1-αN-m-μi.Analogously, we have(89)wi2=-2μi∑m=2N1m-1Nmαm1-αN-m,wi3=μi∑m=2N1m-1Nmαm1-αN-m.Hence, it follows that(90)wi=-μi1+αN-11-αN-1.Since α>0, we see that limN→∞wi=0. Thus, we conclude that S^B is asymptotically unbiased.A discussion of the variance of the sample truncated covariance matrix is beyond the scope of the present paper. We just observe that, apart from the above calculation, studying the sample properties of the truncated spectrum is made hard by the fact that eigenvalues and eigenvectors of a diagonalizable matrix are intimately related from their very definition; thus such study would require a careful consideration of the distribution of the sample diagonalizing matrixR^≠R of S^B.


## 7.3. Perturbative Regularization ofτρ-1

Whenμ is critically close to the internal boundary of D(τρ-1), a sample estimate μ^ may fall outside of it due to statistical fluctuations. In that case the iterative procedure described in the previous sections diverges. On the quantitative side, the ill-posedness of the reconstruction problem is measured by the failure probability(91)pfailρ,Σ,N=Pμ^∉Dτρ-1∣Xk~Nv0;Σ,k=1,…,N,which is a highly nontrivial function of ρ, Σ, and N. An illustrative example of it is reported in Figure 9(a), which refers to a specific case with v=4 and Σ=diag(0.1,0.3,0.8,2.2). The plot suggests that the iterative procedure becomes severely ill-posed in the regime of strong truncation.Figure 9
(a) Numerical reconstruction of the failure probability of the iterative procedure asv=4 and Σ=diag⁡(0.1,0.3,0.8,2.2), for several values of ρ and for N=200,250,…,1000. (b) Failure probability of the perturbative regularization with the same parameters.
(a) (b)In order to regularize the problem, we propose to go back to (19) and consider it from a different perspective. Specifically, we move from the observation that a simplified framework occurs in the special circumstance when the eigenvalue spectra are fully degenerate, which is essentially equivalent to the setup of [1]. If μ1=⋯=μv≡μ~, by symmetry arguments it follows that λ1=⋯=λv≡λ~ and the other way round. Equation (19) reduces in this limit to(92)μ~=λ~Fv+2Fvρλ~≡Tρλ~.It can be easily checked that the function Tρ(λ~) is monotonically increasing in λ~. In addition, we have(93)ilimλ~→0⁡Tρλ~=0,iilimλ~→∞⁡Tρλ~=ρv+2;thus (92) can be surely (numerically) inverted provided that 0<μ~<ρ/(v+2). We can regard (92) as an approximation to the original problem (19). When μ is not degenerate, we must define μ~ in terms of the components of μ. One possibility is to average them, that is, to choose(94)μ~=1v∑i=1vμi.Subject to this, we expect λ~ to lie somewhere between λ1 and λv. Equation (92) can be thought of as the lowest order approximation of a perturbative expansion of (19) around the point λT={λ~,…,λ~}. If the condition number of Σ is not extremely large, such an expansion is expected to quickly converge, so that a few perturbative corrections to λT should be sufficient to guarantee a good level of approximation.As mentioned above, a technical implementation of the perturbative approach and a thorough discussion of its properties are deferred to a separate paper [27]. Here, we limit ourselves to observing that the definition domain of perturbation theory is ultimately set by its lowest order approximation, since corrections to (92) are all algebraically built in terms of it, with no additional constraints. Following (94), the domain of Tρ-1 comes to be defined as(95)DTρ-1=μ∈R+v:∑i=1vμi≤ρvv+2,and it is clear that D(τρ-1)⊂D(Tρ-1) (it is sufficient to sum term by term all the inequalities contributing to (36)). In Figure 10, we show the set difference DTρ-1∖D(τρ-1) in v=2 and v=3 dimensions. When μ∈D(τρ-1) but its estimate μ^∉D(τρ-1), it may well occur μ^∈D(Tρ-1); that is, the set difference acts as an absorbing shield of the statistical noise. Therefore, if we define the failure probability of the perturbative reconstruction as(96)qfailρ,Σ,N=Pμ^∉DTρ-1∣Xk~Nv0;Σ,k=1,…,N,we expect the inequality qfail(ρ,Σ,N)≪pfail(ρ,Σ,N) to generously hold. An example is given in Figure 9(b): we see that qfail becomes lower than pfail by orders of magnitude as soon as ρ and N are not exceedingly small. In this sense, the operator Tρ-1 can be regarded as the lowest order approximation of a regularizing operator for τρ-1.Figure 10
(a) Set differenceDTρ-1∖Dτρ-1 in v=2 dimensions. (b) Set difference DTρ-1∖Dτρ-1 in v=3 dimensions.
(a) (b)


## 8. Conclusions

In this paper we have studied how to reconstruct the covariance matrixΣ of a normal multivariate X~Nv(0,Σ) from the matrix SB of the spherically truncated second moments, describing the covariances among the components of X when the probability density is cut off outside a centered Euclidean ball. We have shown that Σ and SB share the same eigenvectors. Therefore, the problem amounts to relating the eigenvalues of Σ to those of SB. Such relation entails the inversion of a system of nonlinear integral equations, which admits unfortunately no closed-form solution. Having found a necessary condition for the invertibility of the system, we have shown that the eigenvalue reconstruction can be achieved numerically via a converging fixed point iteration. In order to prove the convergence, we rely ultimately upon some probability inequalities, known in the literature assquare correlation inequalities, which have been recently proved in [17].In order to explore the convergence rate of the fixed point iteration, we have implemented some variations of the nonlinear Gauss–Jacobi scheme. Specifically, we have found that overrelaxing the basic iteration enhances the convergence rate by a moderate factor. However, the overrelaxed algorithm still slows down exponentially in the number of eigenvalues and polynomially in the truncation radius of the Euclidean ball. We have shown that a significant reduction of the slowing down can be achieved in the regime of strong truncation by adapting the relaxation parameter to the eigenvalue that is naturally associated with, so as to boost the higher components of the spectrum.We have also discussed how the iterative procedure works when the eigenvalue reconstruction is performed on sample estimates of the truncated covariance spectrum. Specifically, we have shown that the statistical fluctuations make the problem ill-posed. We have sketched a possible way out based on perturbation theory, which is thoroughly discussed in a separate paper [27].A concrete implementation of the proposed approach requires the computation of a set of multivariate Gaussian integrals over the Euclidean ball. For this, we have extended to the case of interest a technique, originally proposed by Ruben for representing the probability content of quadratic forms of normal variables as a series of chi-square distributions. In the paper, we have shown the practical feasibility of the series expansion for the integrals involved in our computations.



---

*Source: 6579537-2017-01-10.xml*