{'fraction_non_alphanumeric': 0.02676622267219287, 'fraction_numerical': 0.010377448033663429, 'mean_word_length': 4.767883211678832, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal Polyspectra Estimation

We show that full inverse covariance weighting can be naturally incorporated into modal estimation methods making them optimal for the CMB power spectrum, bispectrum and trispectrum, as well as in other 3D applications. Modal estimation methods are highly efficient requiring inversion of only an nmax x nmax covariance matrix, where nmax is the limited number of modes needed to describe both the theoretical model and the noise which affects it (in contrast to the full lmax^2 x lmax^2 inverse covariance weighting). A WMAP resolution implementation for modal bispectrum estimation is discussed and shown to yield optimal constraints for local and other models at lmax =500 and nmax<= 50. We also discuss further improvements by characterising and preconditioning the covariance matrix using a single bispectrum variance shape induced by the noise and mask. We describe optimal modal estimation in general terms applicable to polyspectra of any order in a formalism that can be extended beyond the isotropic case. We use the example of anisotropic noise in bispectrum estimation to show the completeness of the isotropic modal basis and the necessity of subtracting a linear estimator term for optimality, removing cross-terms in the variance. We briefly discuss applications to the CMB power spectrum and trispectrum, as well as large-scale structure.


I. INTRODUCTION
A separable eigenmode approach to CMB bispectrum and trispectrum estimation has been developed and implemented [1,2], providing a reconstruction of the WMAP CMB bispectrum and a wide variety of constraints on theoretical bispectra [3] and trispectra [4]. In this paper, we demonstrate how this modal methodology can be made fully optimal by incorporating inverse covariance weighting in the presence of anisotropic noise and other systematic effects. It can be efficiently implemented provided the eigenmode expansion utilised has sufficient resolution to describe, first, the theoretical polyspectrum shape being investigated and, secondly, the noise, mask and other systematic contributions that are correlated with it. The generality and simplicity of this modal covariance inversion, combined with the efficiency of the separable modal estimator, opens up the possibility of fully optimal analysis of high resolution CMB experiments such as the Planck satellite, as well as higher order correlators in the three-dimensional distribution of large-scale structure.
Full inverse covariance weighting for CMB power spectrum analysis has a long history going back to the study of the COBE data using direct Cholesky decomposition with O(l 6 max ) operations (see, for example, [5,6]). This was not feasible for higher resolution experiments such as WMAP, but could be improved to O(l 4 max ) with conjugate gradient methods and iterative approximations [7]. The importance of inverse covariance weighting for the CMB bispectrum has also been discussed for some time (see, for example, [8]) but was only tackled at WMAP resolution more recently using multigrid preconditioning and conjugate gradient methods [9,10]. Achieving full inversion at Planck resolution remains formidable goal for a realistic noise model with non-trivial correlations (see, for example, ref. [11]). The significant advance of the present work is that we show that it is not necessary to invert the full covariance matrix to obtain an optimal variance. Instead, we invert a much smaller modal covariance matrix from a subspace with sufficient degrees of freedom to describe the theoretical model, as well as the masking and noise modes that actually contaminate it. For example, we will investigate an isotropic bispectrum b l1l2l3 affected by anisotropic instrument noise or a cut-sky mask, discussing the projection of the pixel-variance noise map, variance and bispectrum which can affect the estimation. We have previously demonstrated, in the isotropic case, that a highly efficient and general modal bispectrum estimator can be implemented projecting the CMB data into a much smaller subspace spanned by the relevant bispectrum modes [1,2]. What is remarkable is the small number of modes required to characterise almost all popular theoretical bispectrum models n max = O(10 − 100) [12]. For this isotropic subspace, we then noted that the required mode number is not noticeably increased when characterising the relevant noise and mask effects [2]. Our expectation, therefore, when seeking an isotropic theoretical polyspectra estimator on the full anisotropic space that is that we will require few additional modes to incorporate relevant noise and masking effects (this will be the subject of a future publication). Returning to the full isotropic implementation in a WMAP-realistic context, we are able to significantly improve the optimality of the variance using a much simpler inverse modal covariance weighting, getting about 30% closer to that predicted by a Fisher matrix analysis than previous isotropic bispectrum estimators.
Finally we comment on other modal estimators such as wavelets [13], bins [14] and needlets [15] which differ in there choice of initial basis but are otherwise identical. Important distinctions remain, including the ability to analytically calculate modes, the number of modes required for convergence and the ability to project between primordial and CMB(or LSS)-spectra. Early versions of the wavelet method [13], claimed optimality while neglecting to incorporate the linear term in the estimator which is required to minimise the variance [16], though this has now been incorporated [17]. We demonstrate the crucial role of this linear term for optimal estimation in the realistic situation where rotational invariance is broken.
The plan of the paper is as follows. We will first outline a general polyspectrum estimation methodology describing projections from the full polyspectrum space into the modal subspace of interest. In this context we show how the full covariance matrix is projected down into a much smaller modal covariance matrix which can be easily inverted. Next we discuss the modal bispectrum estimator in detail, deriving the relationship between the full and modal covariance matrices. WMAP results are presented comparing analytic optimal Fisher matrix forecasts with the variance obtained from Gaussian simulations for the general isotropic modal estimator. We calculate the bispectrum variance from anisotropic noise and we note that identifying the primary bispectrum noise shape contributing to this variance can be used to efficiently precondition the covariance matrix for inversion. Finally, we give the modal estimation components required for optimal CMB power spectrum and trispectrum analysis and briefly discuss large-scale structure correlators.

II. MODAL CMB POLYSPECTRA ESTIMATION
The CMB polyspectrum a ℘ of degree p for a given non-Gaussian model is the expectation value of the product of p temperature multipole a lm 's obtained from an ensemble of universes, that is, where the index ℘ represents the 2p degrees of freedom ℘ = {l 1 , m 1 , l 2 , m 2 , ..., l p , m p } which cover the full domain V of possible polyspectra (up to a given resolution l ≤ l max ). Here we assume a ℘ is the connected part of the correlator we seek. Most primordial theories predict an isotropic polyspectrum a ℓ which resides on a smaller d-dimensional subspace V I , after integrating out (or summing over) the irrelevant anisotropic modes. The isotropic a ℓ has an equivalent counterpart in the full space V found by multiplying by where ℘ = {ℓ, ℓ ⊥ } with the ℓ index covering the d isotropic degrees of freedom and ℓ ⊥ the remaining 2p−d anisotropic degrees. For concreteness, we note that the power spectrum C l has d = 1 with , as we work with a lm a * lm , while for the bispectrum we have d = 3 with the Gaunt integral , m 2 , m 3 } (and similar for higher 'nondiagonal' polyspectra).
When investigating CMB polyspectra we are typically seeking to match a statistically isotropic prediction a ℘i to the CMB data, with an appropriate weighting to reflect our knowledge of the signal-to-noise for the specific experiment. The estimatorĒ determines the goodness of fit between these through the inner product where the C −1 lpmp,l ′ p m ′ 2 p represents the inverse covariance across all p multipoles. Here, a lin ℘ represents the 'linear term' used to subtract out known systematic effects; for the power spectrum, a lin ℘ = C l1m1,l2,m2 and for the bispectrum a lin ℘ = 3C l1m1,l2,m2 a l3m3 . In practice, the linear term is calculated from Monte Carlo simulations of Gaussian maps with the correct experimental effects (mask, beam, noise) included, see [2]. We also note that a ℘ is understood to be the theoretical prediction modified to include the relevant experimental effects.
If we define the weighted vectors A, B and the matrix C then the estimator can be written in matrix notation as Modal decompositions have proven to be very useful for estimating general CMB bispectra and trispectra. Suppose we have a basis R n℘ where n is the mode number. We will take it to be orthonormal so or in matrix notation RR T = I. We can then represent our theoretical polyspectra in this basis using mode coefficients α defined by The α can then be calculated easily via α = RA.
Most theoretical cosmological models have rotational invariance at their foundation, so the polyspectra A th they predict are generically isotropic. For this reason, our starting point must be to include sufficient isotropic mode functions in our basis to accurately characterise them. Most theoretical bispectra A th have smooth or simple forms and so for a suitable choice of R they can in fact be represented by surprisingly few of these isotropic modes (for example, most scale-invariant bispectrum models require only 15 modes at WMAP resolution [12]). This motivates us to work with a truncated basis containing n max modes which span the space of interest which is determined by the class of theoretical models under consideration modified by experimental effects, i.e. masking and noise. For an isotropic bispectrum estimator, we have shown that the experimental effects of the mask and noise can be incorporated without increasing n max substantially [2]. For anisotropic contributions from experimental effects which yield a partially off-diagonal A, we need to supplement our basis functions with additional anisotropic modes. In line with the remarkable data compression possible in the isotropic case, we anticipate only increasing n max by multiplying by a factor of a few provided the extra anisotropic R are carefully chosen.
Given that n max is so small relative to the vast number 2p of multipole modes, the R become highly rectangular n max ×2p matrices. We can now define a projection operator from the full space V to the subspace V P P ≡ R T R Here, we assume we have truncated our basis such that our theory and distorting experimental effects remain accurately described, ie. PA = A. Using this projection we define the modal counterparts of the data on the same subspace V P , as well as the covariance The projected covariance is related to the covariance of β by p! as we now show.
As we know the a lm are very close to Gaussian we can expand a ℘ a ℘ ′ into the product a lm a l ′ m ′ , a lm a lm and a l ′ m ′ a l ′ m ′ . In general there will be p! ways to have all pairs consist of one primed and one unprimed a lm plus (2p)!/2p! − p! combinations with at least one pairing between two unprimed and two primed a lm . The linear term is always chosen such that it removes all contributions of the second type. As a result we have BB T = p!C and so Using these quantities we can define a new estimator where we have used (RCR T ) −1 = RC −1 R T + Z ⊥ where Z ⊥ is an arbitrary term orthogonal to the projection which we can ignore. Strictly C is non-invertible, formally because it contains duplicate data due to permutations, but also because of loss of information from a cut-sky mask. However, this problem can be circumvented. For example, we can restrict ℘ to count only combinations where l 1 ≤ l 2 ≤ l 3 , with a suitable restriction on m i when the l i are equal, and weight them appropriately to remove this complication. Once it is inverted we can then reinsert duplicates to return to the previous ordering for ℘. Given the estimator expression (14), it can be easily seen that its variance is given by We now wish to determine how closely the standard estimator (5) and the modal estimator (14) coincide. In the idealised case of homogeneous noise and no mask the covariance matrix is diagonal and so C = I and they are clearly identical. Also in the limit of completeness n max → 2p, we have P → I and again the estimators are equivalent. In a realistic context with inhomogeneous noise, a mask and a truncated set of modes the difference can be quantified by considering where we define the vector X ≡ PX to be the component of X in the subspace V P and the vector X ⊥ ≡ (I − P)X to be orthogonal to it and, similarly, for the matrix we can decompose it into M ≡ PMP, M ⊥ ≡ (I − P)M(I − P) and M × ≡ PM(I − P). In this case, we can represent the standard estimator (5) in the form,Ē whereas our new modal estimator (14) is The difference between the estimators is the inverse covariance projection of the perpendicular part of the data onto the parallel space. We should however pause to note that the normalisation in both cases is identical. As the variance we expect for the estimator is the reciprocal of the normalisation (with a factor for permutations), we can see the effect of the cross-projection by comparing the variance we obtain from simulations with this prediction. If they differ we must deduce that our covariance weighting is sub-optimal and the cross-term is important in which case we must increase n max and expand our basis to include the relevant parts.
It should be noted that as the covariance matrix ζ is symmetric positive definite it is possible to find its Cholesky decomposition, whereλ is lower triangular. We can then absorb these factors into our α and β, This is equivalent to re-orthonormalising our basis with respect to the full covariance matrix, and hence the β are now uncorrelated. This is obviously a desirable property and has several applications which we explore elsewhere in future publications. The most obvious is that the estimator takes a very simple form as by construction ζ ′ = I. The computational advantages of the modal estimator (14) are considerable, especially when we note that direct evaluation of the full optimal estimator (5) is extremely challenging at high resolution. This is because, first, it represents a huge sum over l 2p max quantities and, secondly, the inverse covariance weighting requires a huge l 2 max × l 2 max matrix to be calculated and then inverted. As we noted in the introduction, there do exist methods which make this tractable in specific cases, however for large l max and general noise models at Planck resolution this has not yet been achieved. Calculation of the mode coefficients α, β may appear to be troublesome, but it is made highly efficient using a separable form which has been discussed at length previously [2]. Instead, here, we focus on the improvements offered by the modal approach for inverse covariance weighting, reducing this to the inversion of a small n max × n max matrix.

III. CMB BISPECTRUM ESTIMATION
A. Inverse covariance weighting with the modal bispectrum estimator Using the CMB bispectrum, we will now demonstrate a concrete example and implementation of modal polyspectra estimation with inverse covariance weighting. The estimator for the bispectrum takes the general form where B l1l2l3 m1m2m3 = a l1m1 a l2m2 a l3m3 is the full bispectrum, modified from being purely isotropic by experimental effects, and N is the appropriate normalisation. This estimator has been shown to be optimal in the ideal case without mask or noise in the absence of a linear term [18]. The linear term was proposed in [16] to minimise the variance in a realistic setting when rotational invariance is broken. To construct the modal bispectrum estimator, the quantities A, B, C in (4) take the form For simplicity and to relate the discussion more easily to previous work, let us restrict attention to an isotropic subspace, though the expressions derived below can be converted for the general anisotropic case. We shall take the usual approximation that the primary contribution to the isotropic bispectrum is the predicted theoretical model itself, that is, where b l1l2l3 is the predicted reduced bispectrum. Here we note that we are always working with quantities that include the relevant experimental effects. For example, as we are only considering the isotropic case we can approximate b l1l2l3 ≈ f sky b l1 b l2 b l3 b clean l1l2l3 and C l ≈ f sky b 2 l C clean l + N l where b l is the beam window function, N l is the noise power spectrum and f sky is the fraction of the sky remaining after masking.
An isotropic set of basis functions can be represented in the full anisotropic bispectrum space by where v l is a weight function chosen to mimic the scaling of G, the Gaunt integral. The simple orthogonality condition RR T = I in the full space V, yields a weighted orthogonality condition on R in the projected isotropic subspace V P , where h l1l2l3 is a geometric factor defined by h l1l2l3 = (2l 1 + 1)(2l 2 + 1)(2l 3 + 1) 4π and decomposing the A, B, C into mode space yields the α n , β n coefficients and the covariance matrix ζ nn ′ respectively, that is, Calculation of α n , β n is straightforward using methods described in detail in ref. [2] and, as we have chosen the linear term correctly to remove cross terms in the variance, ζ = 1/6 ββ R as we demonstrate below Our inverse covariance weighted estimator (18) becomes This result applies equally to the general anisotropic case where we have supplemented with basis functions R l1 l2 l3 m1m2m3 which cannot be simply represented in terms of the Gaunt integral; the derivation above is unaffected. Using (36), we note that the variance σ 2 of the optimal modal estimator (36) becomes simply

B. Bivector representation of covariance matrix
Further efficiencies are made possible by identifying the underlying bispectrum 'noise' and 'mask' shapes that contribute to the covariance matrix. We can identify such a shape because the noise is correlated in harmonic space. Empirically in our previous isotropic WMAP analysis we have found that the noise and mask contribute in tandem to produce a single shape vector β = {β n } which we could identify [3]. This noise/mask vector is significantly correlated with local non-Gaussianity which explains the difficulty in obtaining optimal constraints on this model. The dominant effect of β leads to the remarkable result that we can accurately represent the isotropic covariance matrix simply as the identity matrix perturbed by this bivector, that is, For this, there is a simple inverse covariance For the isotropic estimator (36) in the presence of this noise/mask bivector the variance from (37) becomes where we have described the degree of correlation between the two vector directions through the cosine From (40), we see that the variance is made up of two components. First, there is the ideal or 'optimal' variance 1/|α| 2 that can be achieved in the presence of the underlying Gaussian CMB signal and homogenous instrument noise, given the degree of sky coverage and beam resolution (incorporated into the signal-to-noise definition of α). Secondly, there is the contribution from the inhomogeneous noise and the anisotropic mask which arises if the noise/mask vectorβ is closely aligned with the model α under investigation. Significantly, this second contribution vanishes if the model being studied has little correlation withβ, in which case nearoptimal variance can be approached even in this isotropic case. In principle, in the full anisotropic case, the general inverse covariance weighting will approach optimality given sufficient modes for any model.
Another key point in a practical implementation is that the single noise/mask vectorβ can be estimated accurately from relatively few simulations. Usingβ with the covariance matrix in the form (38) allows accurate and efficient inversion because we only need to determineβ. This can be achieved with n sim 1000 simulations, while direct averaging of the covariance matrix takes n sim ≈ 10000 to converge to percent levels (making it reliably invertible). We remark that early studies of these contributions in a WMAP-realistic context showed that the local shape is strongly correlated with the spurious contributions from the mask and inhomogenous noise, with the noise/mask vectorβ showing an 70% correlation with the local α loc direction. Separately, there was about an 70% local-mask correlation and a 47% local-noise correlation. In a forthcoming publication, we shall discuss the more general case where several 'noise' vectors can be obtained using a principal component analysis (PCA), separating the noise, mask and other contaminant vectors from model directions of theoretical interest.
Finally, we note that for the simple bivector form of the inverse covariance matrix (39), the Cholesky decomposition of ζ can also be expressed analytically. Here, we have ζ =λλ T withλ the lower triangular matrix:λ Rotation with this matrixλ using (20), offers an orthonormal frame in which the bispectrum estimator (36) becomes simply E = α ′ · β ′ / |α ′ | 2 , that is, the effects of the noise and mask are incorporated into the definitions of α ′ and β ′ from the outset.

C. Isotropic inverse covariance implementation
The practical implementation of the inverse covariance weighting in (18) has been incorporated in the isotropic case for the modal bispectrum estimator. Bispectrum constraints for a wide variety of models using WMAP data have been presented elsewhere [3]. Here, our purpose is to focus on the improvement in optimality for the variance ∆f NL and we only present results for local, equilateral and orthogonal non-Gaussian models. We compare the variance achieved for the modal isotropic estimator (which agrees well with previous isotropic KSW estimators) to show that it can be improved with modal inverse covariance weighting in this isotropic case. Of course, we do not achieve full optimality, which can be seen by comparing to the Fisher matrix forecast, because relevant anisotropic effects deriving chiefly from the mask cannot be described using only an isotropic basis.
For the purposes of direct comparison, we present the variance ∆F NL using the integrated bispectrum F NL defined in ref. [3], that is, the total bispectrum signal-to-noise normalised relative to the local model (for which f loc NL = F NL ). For an arbitrary predicted theoretical bispectrum b th l1l2l3 , we use the estimator (22) to define where the universal normalization N has been defined using with N fNL=1 loc being exactly the same quantity except defined for the f NL = 1 local model. The quantity F NL provides a much more uniform variance between different models, though we also present results with the usual equilateral normalisation f equil NL .  Table I: Results for the variance from an implementation of the isotropic modal estimator in a WMAP realistic context. These results were obtained from 1000 Gaussian maps at lmax = 500 using the KQ75 mask and the WMAP inhomogeneous noise model. The mask and noise degrades the variance relative to the Fisher matrix forecast, but this is significantly improved using modal inverse covariance weighting, even for this isotropic estimator.
Results for the variance ∆F NL from local and equilateral models are listed in Table I. This work conforms exactly to the WMAP-realistic methods described in ref. [3], using a KQ75 mask, WMAP beams and inhomogeneous noise model (see below) and with a resolution limit of l max = 500. The results were obtained using 144000 Gaussian maps from a modal bispectrum pipeline previously confirmed to be unbiased and which reproduces the correct f NL for simulated non-Gaussian maps [2]. Previous results for ideal full-sky Gaussian simulations (without inhomogenous noise or mask) demonstrate that the estimator is correctly normalised, confirming that optimality variance in the absence of experimental effects. For the local model in a WMAP-realistic context, the variance from the modal estimator is suboptimal by 30% dropping from the Fisher matrix forecast ∆F NL = 22.9 to a realised value ∆F NL = 29.5. This is consistent with the value obtained from KSW estimators previously at the same l max = 500 (e.g. ∆F NL ≈ 26 in ref. [19]). This strong degradation for the local model occurs because the noise/mask vector (38)) correlates strongly with the local shape, increasing the variance as discussed previously (see (40)). The addition of inverse covariance weighting for the local model, improves the variance to ∆F NL = 27.6 which takes us 30% closer to the projected optimal variance, demonstrating the value of incorporating this weighting even in the purely isotropic case. In contrast, nearly optimal results can always be obtained for the equilateral model and the orthogonal model because these are much less correlated with the noise/mask vectorβ. For the equilateral model, in particular, an isotropic estimator will be capable of achieving near-optimal results.

IV. OTHER POLYSPECTRA AND 3D APPLICATIONS
This methodology can be extended to other correlators. We begin with the trispectrum estimator presented in [20], for which WMAP constraints were obtained for several models in [4]. We shall restrict ourselves to the diagonal free trispectrum for simplicity, although obtaining expressions for the general case is straight forward. We build the optimal modal estimator our of the quantities A, B, C in (4) defined by where v l is a suitable weight function and the isotropising factor G l1l2l3l4 m1m2m3m4 is Using efficient modal methods calculation of the trispectrum estimator with inverse covariance weighting is no more challenging than in the bispectrum. Inverse covariance weighting for power spectrum estimation for high resolution experiments like Planck. In this case we would be constraining small deviations from a canonical model with a given C l . Considering an isotropic subspace, the building blocks for the model estimator are whereC l is understood to be the canonical model including the predicted deviation. Also we now consider C l1m1,l2m2 = a l1m1 a * l2m2 . The corresponding isotropic basis functions are which must be supplemented to incorporate anisotropic experimental effects such as cut-sky masking. We can exploit the tractability of the modal inverse covariance weighting in several ways. First we could consider constraining specific models with specific signatures in the power spectrum with a particular α n , such as oscillatory models. Secondly we could use the method for full power spectrum estimation using the recovered β to iteratively improve the best fit theoretical model. This would entail linking to Monte Carlo simulations which explore cosmological parameter space. Finally, we need not restrict ourselves to isotropic models and instead using additional anisotropic basis functions, we can to constrain anisotropy in the CMB. For large-scale structure and other 3D polyspectra estimation, this approach has an entirely analogous implementation discussed in ref. [21]. We simply replace our multipoles with the density perturbation so now we have where now the vector ℘ represents ℘ = (k 1 , k 2 , ..., k n ) which cover the full 3n-dimensional domain of possible polyspectra. Here momentum conservation and isotropy imply we can make the reduction to the subspace ℘ = (k 1 , k 2 , ..k n ) and the two are related by The rest of the modal method then follows analogously with two main caveats. The first being that the growth of structure is non-linear so even with gaussian initial conditions we produce measurable polyspectra. The form of these can be calculated from N -body simulations and we can seek deviations from these due to primordial and other effects. Secondly there is a non-linear relationship between the late time and primordial polyspectra and so a Monte Carlo approach is required to find the parameter values which provide the best fit to the data. Error bars can then be obtained from simulations with the chosen primordial spectrum. However, the inverse covariance matrix can still be calculated from Gaussian simulations, so it does not need to be recalculated for differing levels of primordial non-Gaussianity.

V. CONCLUSION
We have shown that full inverse covariance weighting can be naturally incorporated in the modal polyspectra estimator to achieve optimality. The resulting modal covariance matrix can be efficiently calculated and inverted. This greatly reduces the computational cost and memory requirements of polyspectra estimation, making optimal analysis of high resolution experiments like Planck possible. We emphasise that the primary efficiencies discussed here arise because we have effectively projected out all but the modal degrees of freedom needed to characterise the model (or models) under investigation, while retaining an adequate description of the noise and systematics with which it is correlated.
We have only briefly discussed the many applications of this work. These include improved WMAP7 constraints on a wide variety of theoretical models, together with a model-independent constraint on the overall magnitude of the bispectrum in ref. [3]. A principal component analysis can be used to identify the noise, mask and other shapes for preconditioning the inverse covariance weighting. A forthcoming paper will use eigenmodes to characterise and marginalise over foregrounds and other contaminants in preparation for the CMB analysis of Planck data. The analysis of higher order correlators in three dimensions represents a particular challenge as many billion galaxy surveys become available. The 3D modal estimator [21] allows for efficient polyspectra estimation and for the setting of arbitrary non-Gaussian initial conditions. Combined with inverse covariance weighting, this methodology may help unravel the many systematic and evolutionary effects present in large-scale structure observations. optimal estimation. WMAP non-Gaussian analysis to date has used a simple instrument noise model which assumes no correlations between pixels. If n lm is the harmonic transform of the pixel-variance map, then the covariance of the noise contribution to the CMB temperature is given by a N l1m1 a N l2m2 = 4π n pix lm (−1) m G l1 l2 l m1m2−m n lm . (A1) (See, for example, the discussion of covariance under rotations in [22].) The covariance of the cut-sky mask also can be represented in this form by taking the masked pixels to have a very large noise so that the covariance weighting neglects that section of the map. We shall now consider the contributions to the variance of the estimator both with and without the linear term. We assume CMB signal plus instrument noise can be represented as a lm + a N lm with respective covariances C S l1m1,l2m2 = a l1m1 a l2m2 = (−1) m1 C l1 δ l1l2 δ m1−m2 and C N l1m1,l2m2 = a N l1m1 a N l2m2 . Again restricting attention to the isotropic subspace by expanding using isotropic modes (27), the modal covariance matrix with the inhomogeneous noise takes the form