Science-driven 3D data compression
David Alonso
astro-ph.CO



^1Oxford Astrophysics, Department of Physics, Keble Road, Oxford, OX1 3RH, UK


  Photometric redshift surveys map the distribution of matter in the Universe through the positions and shapes of galaxies with poorly resolved measurements of their radial coordinates. While a tomographic analysis can be used to recover some of the large-scale radial modes present in the data, this approach suffers from a number of practical shortcomings, and the criteria to decide on a particular binning scheme are commonly blind to the ultimate science goals. We present a method designed to separate and compress the data into a small number of uncorrelated radial modes, circumventing some of the problems of standard tomographic analyses. The method is based on the transform (KL), and is connected to other 3D data compression bases advocated in the literature, such as the Fourier-Bessel decomposition. We apply this method to both weak lensing and galaxy clustering. In the case of galaxy clustering, we show that the resulting optimal basis is closely associated with the Fourier-Bessel basis, and that for certain observables, such as the effects of magnification bias or primordial non-Gaussianity, the bulk of the signal can be compressed into a small number of modes. In the case of weak lensing we show that the method is able to compress the vast majority of the signal-to-noiseÀ into a single mode, and that optimal cosmological constraints can be obtained considering only three uncorrelated KL eigenmodes, considerably simplifying the analysis with respect to a traditional tomographic approach.


Science-driven 3D data compression
    David Alonso^1
    December 30, 2023
==================================





§ INTRODUCTION

  Astronomical data is inherently three-dimensional: the main observable is the intensity of the electromagnetic emission in the sky as a function of wavelength and line-of-sight direction, determined by two angles. An idealized cosmological analysis would therefore use, as a data vector, the full cube I(λ,θ,ϕ) probed over sufficiently well resolved angular and frequency scales <cit.>. However, the operational costs of obtaining such a dataset imply that we can realistically only access a compressed version of it, where the compression method comes in different flavours:
  
    
  * We can decrease the measurement noise by integrating the sky intensity over large frequency bands. This approach has been used for instance in CMB observations <cit.>.
    
  * Angular resolution can also be sacrificed for wider sky and frequency coverage, as has been proposed for future intensity mapping experiments <cit.>.
    
  * The size of the dataset can also be reduced by collecting only the flux associated with the brightest extra-galactic objects. An accurate measurement of their spectra then allows a determination of their redshift, producing the well-known spectroscopic redshift surveys <cit.>.
    
  * Measuring individual spectra is a costly operation, however, and can usually only be done for a small subsample of all available sources. This problem can be mitigated by inferring the source's redshift from their emission in a small number of wider frequency bands, in what is known as a photometric redshift survey <cit.>. The redshifts thus determined are far less precise than their spectroscopic counterparts, and usually only an imperfect estimation of the redshift probability distribution for each galaxy is accessible.
  

  Even after this first compression stage, the size of the dataset makes a direct analysis of it as a data vector a computationally intractable problem. Typically this should not be an issue in terms of information loss, since large portions of the data are usually dominated by measurement noise, contaminated by sources of systematic uncertainty (both observational and theoretical) or contain only redundant information. An efficient data compression method will therefore identify these sections of data space and eliminate them, or collect them into summary statistics, while minimizing the loss of meaningful cosmological information. An example of this is the standard analysis of cosmological datasets in terms of their two-point statistics <cit.>. However, even in this case the resulting data vector can be large enough to present an important computational challenge in terms of likelihood evaluation and covariance estimation. The complexity of latter problem, in particular, scales with the square of the data vector size, and can become an important drain of computational resources <cit.>.
  
  In this work we will concern ourselves with the topic of 3D data compression: the problem of identifying the uncorrelated angular and radial modes of the data that optimally contain the maximum amount of information. This problem has been previously addressed in the literature <cit.>, and a number of approaches have been proposed depending on the definition of uncorrelatedness used, and on the type of information one wishes to preserve. Here we will present a method to derive a set of uncorrelated radial eigenmodes that are manifestly optimal in terms of information compression for any quantity, such as individual cosmological parameters or the amplitude of the cosmological signal over any set of known contaminants. The method is based on the well-understood transform <cit.>, and is similar in spirit to the derivation of optimal weighting schemes for the analysis of spectroscopic surveys <cit.>. Although we will focus here on the case of photometric redshift surveys, the method can be applied to any set of cosmological datasets.
  
  The article is structured as follows: Section <ref> describes the (KL hereupon) transform and its applicability in the context of 3D data compression. Section <ref> shows the performance of the method in a number of science cases, such as the derivation of optimal radial bases for galaxy clustering (Sec. <ref>) and weak lensing (Sec. <ref>) observations and the use of the KL eigenmodes to measure the effects of primordial non-Gaussianity (Sec. <ref>) and lensing magnification (Sec. <ref>) with a small number of modes. Finally Section <ref> summarizes our findings and dicusses the advantages and shortcomings of the method.



§ METHOD

  

 §.§ The Karhunen-Loeve transform

    The idea behind the transform, as developed within the field of cosmological data analysis in e.g. <cit.>, is to compress a given data vector into a small set of modes containing most of the useful information on a particular parameter (or set of parameters). Let x be a data vector of dimension N_s, and let θ be a particular parameter we want to measure. Under the assumption that x is Gaussianly distributed with mean 0 and covariance C, a set of linear combinations y_p≡ e_p^†  x can be found such that the y_p are white and uncorrelated (⟨ y_py_q^*⟩=δ_pq), and such that the first m<N_s combinations contain most of the information about θ. This is done by solving the generalized eigenvalue problem <cit.>:
    
    ∂_θ C  e_p=λ_p  C  e_p,

    where ∂_θ≡∂/∂_θ.
    
    Although the transform can be used to compress the information on any particular parameter, it has been most commonly used to separate signal-dominated and noise-dominated modes by optimizing for the amplitude of the signal, as we explore below. Before moving on, however, it is worth noting that a generalized eigenvalue problem such as Eq. <ref> can always be recast as a standard eigenvalue problem of the form A ẽ_p=λ_p ẽ_p, where
    
    A≡ C^-1/2 (∂_θ C)  C^-1/2,ẽ_p≡ C^1/2 e_p,

    and we have made use of the fact that C is positive-definite (and therefore C^1/2 is well-defined and invertible).

    
    

  §.§.§ The KL transform for the signal-to-noise

      Let us decompose the data vector x into uncorrelated signal and noise components x= s+ n where, in this context, the signal is the part of the data containing any information of cosmological interest, and the noise is any contaminant preventing us from accessing it[n could include, for instance, the contribution of foregrounds in intensity mapping experiments, which motivates the use of the KL transform as a foreground cleaning method <cit.>.]. In this particular case, the data covariance matrix can be split into their independent contributions C= S+ N.
      
      The KL transform has traditionally been used to design an eigenbasis that maximizes the overall signal-to-noise ratio (e.g <cit.>). This can be done by defining a fictitious parameter ρ multiplying the signal part of the data with fiducial value ρ=1 (i.e. x=ρ s+ n). In this case, after some trivial manipulations, the eigenvalue equation (Eq. <ref>) takes the form:
      
    ( S+ N) e_p=λ_p N e_p,

      where we have redefined 2/(2-λ_p)→λ_p. This can be cast into a standard eigenvalue equation using the Cholesky decomposition of the noise covariance matrix N= L L^†:
      
    [ L^-1 C ( L^-1)^†] ẽ_p=λ_pẽ_p,

      where ẽ_p≡ L^† e_p.
      
      At this point it is worth noting that the generalized eigenvalue problem in Eq. <ref> can be understood as the problem of diagonalizing C under a non-standard dot product ∘ given by the inverse noise covariance matrix (i.e. a∘ b≡ a^† N^-1 b). Under this dot product, an eigenbasis F≡( f_1, f_2,..., f_N_s) can be found such that F is orthonormal F∘ F= I, and the covariance of the transformed data vector y≡ F∘ x is diagonal:
      
    ⟨ y  y^†⟩= F^† N^-1 C N^-1 F=Λ≡ diag(λ_1,...,λ_N_s).

      Using the orthonormality of F (with respect to the non-standard dot product), this can be cast into the same form as Eq. <ref>, where f_p= Lẽ_p= N e_p.
      
      Finally, note that, because both S and N are positive-definite matrices, their eigenvalues will also be positive. Since the eigenvalues of N under the KL transform are, by construction, 1, the elements of Λ above will all be greater than 1, and converging to 1 for the noise-dominated modes.
      
    

  §.§.§ The KL transform with correlated contaminants

      Let us now consider a more general case in which we further split the noise into two parts n→ n+ m, where m is a contaminant with a non-zero correlation with the signal. The covariance matrix of the data is then given by:
      
    ⟨ x  x^†⟩=ρ^2 S+2ρ M_s+ M+ N,

      where M_s≡(⟨ m  s^†⟩+⟨ s  m^†⟩)/2, M≡⟨ m  m^†⟩ and we have kept the fictitious parameter ρ defined in the previous section.
      Eq. <ref> then reads:
      
    ( S+ M_s)  e_p=
            λ_p/2 C  e_p.

      
      Unfortunately, in this case the manipulation that lead us to Eq. <ref> cannot be performed. If we were to do so, the matrix remaining on the right hand side of this equation would not be positive-definite, and the corresponding generalized eigenvalue problem would be ill-defined. This is not a problem, since the solutions to Eq. <ref> still separate the modes with the highest signal. The separation of the noise-dominated modes becomes less obvious, however, since the resulting eigenvalues cannot be simply compared with 1, corresponding to noise-dominated modes in the previous section.
      
      The eigenvector solutions to the generalized eigenvalue problem in Eq. <ref> can be collected as columns of a matrix E that simultaneously satisfies the equations:
      
    E^†( S+ M_s) E=Λ,
             E^† C E= I,

      where I is the identity and Λ= diag(λ_1,...,λ_N_s).
      Since the second equation implies C E≡( E^†)^-1, the original vector x can be recovered from the coefficients y≡(y_1,...,y_N_s) as x= C  E  y. More interestingly, one can identify the principal eigenvectors of the Eq. <ref> (e.g. those with associated eigenvalues λ_p above a given threshold λ_ thr) and project out the remaining modes, which are presumably more contaminated by m. This procedure defines a filter W≡ C  E  P  E^†, where P is a projection matrix with 1s in the diagonal elements corresponding to the principal eigenmodes and zeros everywhere else. The filtered data vector is therefore x̃= W x.

  

 §.§ Application to tomographic datasets

    The standard method to draw cosmological constraints from photometric redshift surveys is to divide the galaxy sample into bins in photo-z space and use the information encoded in all the relevant auto- and cross-correlations between different bins <cit.>, making use of various calibration methods in order to estimate the true redshift distribution of each bin. Several criteria can be followed in order to select these redshift bins, such as minimizing the correlation between non-neighbouring bins or preserving a roughly constant number density on all bins. Other approaches
    (<cit.>)
    involve projecting the main observable (e.g. galaxy overdensity or shear) onto the Fourier-Bessel eigenbasis. None of these schemes are manifestly optimal from the point of view of S/N, final cosmological constraints or contaminant deprojection, however. This section presents an alternative slicing scheme addressing these shortcomings, based on the KL transform.
    
    

  §.§.§ Tomographic analyses

      Let us start by assuming that we have split the galaxy sample into N_s subsamples. As mentioned above, we will think of each of these subamples as some kind of redshift binning (e.g. binning galaxies in terms of their maximum-likelihood redshift), but the formalism applies to any set of subsamples. Let a^α() be the a field on the sphere at the angular position  and defined in terms of the properties of the sources in the α-th sample (e.g. the cosmic shear field γ^α or the galaxy overdensity δ^α), and let ϕ^α(z) be the redshift distribution of these sources. Finally, let a^α_ℓ m be the spherical harmonic coefficients of a^α[Spin-2 fields, such as the cosmic shear, will be decomposed in spin-2 spherical harmonics, however the discussion below holds for fields of arbitrary spin.]. The power spectrum for our set of subsamples is defined as the two-point correlator of a^α_ℓ m:
      
    ⟨ a_ℓ m  a^†_ℓ' m'⟩≡δ_ℓℓ'δ_mm' C_ℓ,

      where we have packaged a^α_ℓ m as a vector for each (ℓ,m): a_ℓ m≡(a^1_ℓ m,...,a^N_s_ℓ m). Usually the observed field can be decomposed into uncorrelated signal and noise components a= s+ n, with a similar decomposition in the power spectrum, C_ℓ= S_ℓ+ N_ℓ.
  
      Once the choice of subsamples α is made, the standard analysis method would proceed by performing a likelihood evaluation of the two-point statistics of these subsamples. While this procedure is relatively simple, it suffers from a number of drawbacks, an incomplete list of which is:
      
        
  * It is not clear what the optimal strategy should be to define the sub-samples. The brute-force solution to make sure one exploits all of the information present in the data would be to use a large number of very narrow redshift bins, and let the likelihood evaluation pick up the information encoded in them.
        
  * C^αβ_ℓ is a N_s× N_s× N_ℓ data vector. Thus increasing N_s will increase the computational time required for each likelihood evaluation like N_s^2 and number of elements of the covariance matrix of C^αβ_ℓ like N_s^4, with the corresponding increase in complexity needed to estimate this covariance. Although this can be partially alleviated by considering only correlations between neighbouring redshift shells, the amount of information lost by neglecting all correlations beyond a given neighbouring order is not clear a priori.
        
  * Estimating the redshift distribution for a large number of subsamples can be inaccurate, depending on the method used to do so, on the quality of the photometric redshift posterior information and on the statistics of the available spectroscopic sample.
     
  
    

  §.§.§ Optimal radial eigenbasis

      Following the description in Section <ref>, it is straightforward to derive an optimal set of radial, uncorrelated eigenmodes.
      
        
  * We start by assuming that the field a has been measured in a number of narrow redshift bins, and by defining the inverse-variance weighted field ã_ℓ m≡ N^-1_ℓ  a_ℓ m.
        
  * Let us consider a set of linear combinations of the weighted field measured on narrow redshift bins:
        
    b_ℓ m= F_ℓ^†·ã_ℓ m≡ F_ℓ∘ a,

        where F_ℓ is a yet-unspecified matrix and, as in Section <ref>, we have let N^-1_ℓ define the non-standard dot product v_ℓ∘ w_ℓ≡ v^†_ℓ· N^-1_ℓ· w_ℓ. The power spectrum for this new observable would then simply be given by:
        
    D_ℓ≡⟨ b_ℓ m  b^†_ℓ m⟩= F_ℓ^†∘ C_ℓ∘ F_ℓ.

        
  * Requiring that the new modes be uncorrelated, we can identify Eq. <ref> with the generalized eigenvalue equation <ref>, which defines the KL eigenbasis F_ℓ by additionally requiring that it be orthonormal (F_ℓ∘ F_ℓ= I). Note that, after this transformation and without any further optimization, some of the practicalities of the original problem are already simplified, since we can now focus on the diagonal elements of the new power spectrum and its covariance.
        
  * The data can be further compressed by assuming that we are interested in measuring a set of cosmological parameters Θ≡{θ_1,...}. The information regarding this set of parameters encoded in a given data vector x can be quantified in terms of its Fisher matrix (the expectation value of the Hessian of the log-likelihood with respect to Θ), which assuming ⟨ x⟩=0 reads
        
    F_ij≡⟨∂_i∂_j L⟩=1/2 Tr(∂_i C  C^-1∂_j C  C^-1),

        where C≡⟨ x  x^†⟩ is the covariance matrix of the data.




        We can thus rank the eigenvectors ( F_ℓ)^p_α in terms of their information content (in a Fisher-matrix sense). In the simplest scenario one may be interested in maximizing the overall signal-to-noise (S/N), in which case each mode contributes independently to the Fisher matrix element of the signal amplitude.
        
  * The final set of uncorrelated modes can then be truncated to the first M defined by this procedure, which will contain the bulk of the information needed to constrain Θ.
    
    Besides the elegance of this method in defining a natural set of radial basis functions for the particular dataset under study, analogous to the Fourier-Bessel basis in a translationally-invariant system (see Section <ref>), its merits are better evaluated in terms of data compression. This strategy allows one to reliably and significantly reduce the dimensionality of the data vector from N_s^2× N_ℓ to M× N_ℓ while minimizing the loss of information. This can lead, for instance, to a substantial reduction of the computational costs of likelihood sampling and covariance estimation.
    
    Note that, although the method is based on an initial thin-slicing of the galaxy distribution, the fact that the final datased comprises only a small set of samples means that the method is not penalized in terms of photometric redshift uncertainties. Once the KL eigenmodes F_ℓ are found for a fiducial cosmological model, they can be directly applied as weights to all the objects in the survey to generate the b^p modes. These modes are be characterized by their own window function:
    
    ϕ̃_ℓ^p(z)=∑_α( F_ℓ)^p_α ϕ^α(z)/N^αα_ℓ,

    where we have assumed a diagonal noise power spectrum for simplicity.
    The same methods used to calibrate photo-z uncertainties in the standard tomographic analysis can be applied on b^p to calibrate ϕ̃^p with minor modifications (e.g. weighed and ℓ-dependent stacking of photo-z pdfs, or cross-correlations of the b^p maps with a spectroscopic survey in the case of clustering redshifts). Furthermore, using F_ℓ for the fiducial cosmology as model-agnostic weights and inserting them in Eq. <ref>, the theoretical prediction for the power spectrum of each mode D^p_ℓ can be computed in a model-independent way.
    
    Finally, the method outlined in this section is based on the KL decomposition that maximizes the amplitude of the signal under study. This is the main application advocated in this article, since it is plausible that the set of modes containing the bulk of the cosmological signal will also drive the constraints on any comprehensive set of cosmological parameters (we explore this in more detail in Section <ref>). However, we must note that, for individual parameters, the optimal degree of data compression is achieved by solving the general KL eigenvalue problem (<ref>), which can lead to substantial improvements with respect to the S/N-optimal basis. We explore one particular example of this in Section <ref>.



§ PERFORMANCE AND PARTICULAR EXAMPLES

  This Section explores the performance of the KL decomposition in a number of specific science cases.
  
  

 §.§ Special case: the harmonic-Bessel basis

    Let us start by considering a simplified case where the field a is the overdensity field of a non-evolving galaxy population for which we neglect the effects of redshift-space distortions. Let us further assume that we have perfect redshift information, such that we can split the sample into thin radial slices of equal width δχ, which we label by their comoving radius χ. The noise in the measurement of a is given purely by shot noise, and since (as per our initial assumptions) the number density of sources does not change with χ, the noise power spectrum is diagonal and scales like N_ℓ(χ,χ')∝δ_χ,χ' χ^-2. Thus, the dot product is just given by:
    
    b^†∘ c∝∫ dχ χ^2 b(χ)^* c(χ).

   
    In this case, the cross-shell signal power spectrum is given by,
    
    S_ℓ^χχ'=2/π∫_0^∞ dk k^2 P_k j_ℓ(kχ)j_ℓ(kχ'),

    and it is trivial to show that the KL eigenmodes are simply given by the spherical Bessel functions: ( F_ℓ)^k_χ∝√(2/π)j_ℓ(kχ):
    
    
    D_ℓ^kk'   ≡∑_χ,χ'(F_ℓ)^k_χ(F_ℓ)^k'_χ'S_ℓ^χχ'
       ∝2/π∫ dχ χ^2∫ dχ' χ'^2 j_ℓ(kχ)j_ℓ(k'χ') S_ℓ^χχ'
       =∫ dq q^2P_q[2/π∫ dχ χ^2 j_ℓ(qχ)j_ℓ(kχ)][2/π∫ dχ' χ'^2 j_ℓ(qχ')j_ℓ(k'χ')]
       =∫ dq q^2P_qδ(k-q)/q^2δ(k'-q)/q^2=P_kδ(k-k')/k^2=P_k/k^2Δ kδ_k,k'

    
   
    This choice of basis defines the so-called harmonic-Bessel (or Fourier-Bessel) decomposition, and has been used as a data-compression method for the analysis of photometric redshift datasets (e.g. <cit.>). In any realistic scenario – e.g. in the presence of redshift uncertainties, redshift-space distortions or in the analysis of weak lensing data – this basis is non-optimal (among other things different k-modes will be correlated), as opposed to the KL basis described in the previous section.
    
  

 §.§ Galaxy clustering - Bessel-like eigenfunctions

    
    The assumptions used in the previous section are an ideal limit of the data collected by a photometric survey. In a more realistic (although still idealized) scenario, the information about the radial position of a given source is encoded in its posterior photo-z distribution p(z|α), where α is a continuous variable determining the properties of the photo-z (e.g. the mean of the posterior). The cross-power spectrum of two samples with photo-z properties α and β is given by
    
    C_ℓ^αβ=S^αβ_ℓ+N^αβ_ℓ,
       S_ℓ^αβ=2/π∫_0^∞ dk k^2 Δ_ℓ^α(k) Δ_ℓ^β(k),
       N_ℓ^αβ=δ(α-β)/n_t p(α),

    where n_t is the total angular number density of sources and
    
    Δ^α_ℓ(k)≡∫ dz p(z|α) Ψ_ℓ(k,z) √(P(k,z)),
       Ψ_ℓ(k,z)=b^α(z)j_ℓ(k χ(z))-f(z)j_ℓ”(k χ(z)).

    Here b^α(z) is the linear galaxy bias, f(z)=dlogδ/dlog a is the growth rate of structure, P(k,z) is the matter power spectrum at redshift z, p(α) is the probability that a source has photo-z properties α, and p(z|α) is the conditional redshift distribution of these sources (we have labelled this quantity ϕ^α(z) in previous sections). Note that, for simplicity, we have kept the contribution of redshift-space distortions at linear order and neglected the effect of magnification (this will be studied in Section <ref>).
    
    For a continuous variable α, the generalized eigenvalue problem in Eq. <ref> becomes a homogeneous Fredholm integral equation of the second kind:
    
    ∫ dβ C^αβ_ℓ e^p_ℓ(β)=λ_p∫ dβ N^αβ_ℓ e^p_ℓ(β)⇒
    ⇒∫ dβ n_t p(α) S^αβ_ℓ e^p_ℓ(β)=(λ_p-1)e^p_ℓ(α).

    In the limit of perfect photo-zs (p(z|α)=δ(z-α)), and in the absence of redshift-space distortions, the solution to this equation are the spherical Bessel functions, as proven in the previous section. For general kernels, however, no analytical solution to the homogeneous Fredholm equation can usually be found, and the standard procedure to solve it is through discretization, which is equivalent to taking finite bins in α. We will use this method here to find the KL eigenmodes that maximize the signal content for galaxy clustering.

    To do so, we have considered a specific strawman photometric survey targeting a sample of red galaxies, characterized by their higher bias and better photo-z uncertainties than their blue counterparts (making them better suited for clustering analyses). The sample we consider is compatible with what could be observed by the Large Synoptic Survey Telescope <cit.>, characterized by the redshift distribution shown in the left panel of Fig. <ref> (full details can be found in <cit.>). We assume a photo-z uncertainty of σ_z=0.02 (1+z) and split the sample into redshift bins in photo-z space with z_ ph>0.5 and a width given by the photo-z uncertainty at the bin centre. The auto-power spectra for our set of 23 bins are shown in the right panel of Fig. <ref>. The large overlap between bins implies that a choice of thinner slices is unlikely to unveil significantly more information, and we have verified that the results shown below do not change after doubling the number of bins. All power spectra were computed using a modified version of the code presented in <cit.>.
    
        

    Using the prescription described in Section <ref>, we find the KL eigenmodes and associated power spectra, and rank them according to their contribution to the total signal-to-noise ratio (defined here as the Fisher matrix element of the signal amplitude). The power spectra of the resulting KL modes are shown in Figure <ref>. Unlike the case of weak lensing, explored in Section <ref>, the information encoded in the galaxy overdensity is local in redshift, and thus the correlation between different bins decays rapidly with redshift separation. The signal-to-noise is therefore spread over ∼15 signal-dominated modes, and the noise-dominated modes can be thought of as the radial scales filtered out by the finite photo-z uncertainty (as we mentioned in Section <ref>, the noise power spectrum gets mapped into 1 under the KL transform). The relative contribution of each mode to the total signal-to-noise is shown in the top panel of Fig. <ref>. 90% of the total constraining power can be achieved by considering the first 13 eigenvectors. The form of the first 7 of these eigenvectors for ℓ=30 are shown in the right panel of Fig. <ref>. The eigenmodes are sinusoids with increasing frequencies, in agreement with the expectation that, in the limit of σ_z→0 and no background redshift dependence, the KL decomposition is achieved by the spherical Bessel functions. A Fourier-Bessel decomposition is therefore probably a near-optimal analysis method for galaxy clustering, although the KL decomposition allows a more precise determination of the truly orthogonal radial modes.
    
  

 §.§ Galaxy clustering - optimal basis for f_ NL

    It is expected that future large-scale photometric surveys will make the search for primordial non-Gaussianity one of their main science cases. This can be achieved by measuring the excess power on large scales caused by a non-zero value of f_ NL[The reader is referred to <cit.> for a thorough review of non-Gaussianity and a definition of f_ NL.] generates in the two-point statistics of biased tracers of the matter distribution <cit.>. Since the signal is most relevant on large scales, we can expect the bulk of it to be concentrated in a small number of radial modes, which makes the general KL decomposition outlined in Section <ref> an ideal analysis method. Similar approaches have been explored in the literature to devise optimal weights for spectroscopic galaxy surveys <cit.>.
    
    We again consider the red galaxy sample used in the previous section, but now estimate the KL basis of eigenmodes that optimize the information content on f_ NL instead of the overall signal amplitude. I.e. we solve the generalized eigenvalue problem in Eq. <ref> where θ=f_ NL. We compare the performance of this basis with other choices of radial modes as follows: for a given number of modes, we estimate the associated uncertainty on f_ NL, σ(f_ NL), by summing the contributions to the corresponding Fisher matrix element of those modes, and compute the excess of σ(f_ NL) with respect to the best achievable constraint σ_ best(f_ NL). The results are shown in Fig. <ref> for three choices of radial functions:
    
      
  * The KL eigenbasis resulting from optimizing the information content on f_ NL discussed in this section (in red).
      
  * The KL eigenbasis resulting from optimizing the overall signal-to-noise of the galaxy clustering signal, as discussed in the previous section (in gray).
      
  * Photo-z tomography: the result of dividing the galaxy sample into a number of top-hat photo-z bins of equal width (in blue).
    
    As demonstrated by this figure, for a fixed number of modes the optimal KL basis always outperforms any other data compression prescription. In particular, the constraints on f_ NL are only degraded by ∼30% when considering only the first principal eigenmode, and almost 90% of the total constraining power is contained in the first three. Interestingly, a naive tomographic approach achieves the same uncertainty on f_ NL with a smaller number of modes (redshift bins) than the KL eigenbasis for the S/N. However, since the tomographic bins are not orthogonal, unlike the KL modes, for a fixed σ(f_ NL) both KL bases typically outperform the tomographic approach in terms of the size of the associated power spectrum. In any case, this example serves to stress the fact that the optimal radial basis in terms of overal S/N is not necessary optimal in terms of final constraints for cosmological parameters that depend on specific features of the power spectrum.
    
    
  

 §.§ Galaxy clustering - magnification bias
  
    
    Gravitational lensing of the observed galaxy positions alters their clustering pattern. This appears as an extra term in the galaxy clustering transfer function (Eq. <ref>):
    
    Δ^M,α_ℓ(k)=-2ℓ(ℓ+1)∫ dχ W^M,α(χ)j_ℓ(kχ)/k^2a(χ)√(P(k,z(χ))),
       W^M,α(χ)=3H_0^2Ω_M/2∫_z(χ)^∞ dz' ϕ^α(z')2-5 s/2 χ(z')-χ/χ(z')χ,

    where s is the tilt in the number counts of sources as a function of magnitude limit. This effect, commonly labeled “magnification bias” <cit.>, can be used as an alternative measurement of gravitational lensing, through galaxy positions instead of shapes. The contribution of the magnification term is, however, weak in comparison with the density and RSD terms (Eq. <ref>), and therefore its measurement can be hampered by the cosmic variance contribution of these terms.
    
    
    
    One can therefore think of the density and RSD terms as correlated contaminants of the magnification signal, and use the KL formalism described in Section <ref> to devise an optimal basis of radial eigenmodes containing the bulk of its signal-to-noise.
    
    To test this approach we consider, as in the previous section, an LSST-like survey.  Since lensing magnification is an integrated effect, it is less hampered by poor photo-z uncertainties, and it is most easily measured by cross-correlating high-redshift and low-redshift data <cit.>. For this reason, in this case we consider a sample of blue galaxies, with inferior photo-z errors but wider redshift support. Full details can be found in <cit.>. In summary, we consider a sample with ∼40 objects per arcmin^2 with the redshift distribution shown in Figure <ref>. We also approximate the photo-z distributions as Gaussians with a scatter σ_z=0.05(1+z), and divide the sample into 16 top-hat bins in photo-z space with z_ ph<0.5 and widths given by the value of σ_z at the bin center (again, we verified that our conclusions did not change after decreasing the width by a factor 2).
    
    A key property of the magnification bias effect is the fact that, since gravitational lensing is caused by the integrated matter distribution between source and observer, the magnification signals in widely separated redshift bins can be tightly correlated. This is shown explicitly in Figure <ref>. The figure shows the correlation coefficients between the 16 redshift bins, defined as R^αβ_ℓ=C^αβ_ℓ/√(C^αα_ℓ C^ββ_ℓ), at ℓ=400, with (bottom panel) and without (top panel) the magnification bias effect. Although the contribution of lensing magnification to the correlation between neighbouring bins is subdominant, it produces noticeable correlations between distant ones.
    
    This property is particularly interesting in the context of the KL decomposition: a signal that is tightly correlated across samples will contribute significantly only to a small set of eigenmodes. To explore this possibility, we follow the prescription outlined in Section <ref> for correlated contaminants. The contribution of each eigenmode to the total signal-to-noise of the magnification bias (in a Fisher-matrix sense) is shown in Figure <ref>. As expected, most of the signal (>80%) is contained in the first eigenvalue, with the practical totality of it concetrated in the first three modes.
    
    We finish this section by noting that this approach is similar to the “nulling” method of <cit.>, and that an analogous treatment could be carried out on the cosmic shear field to separate the lensing and intrinsic alignment contributions <cit.>.


  

 §.§ Weak lensing

    
    

    The effects of gravitational lensing can be measured directly by studying the correlation it induces on the shapes and orientation of galaxy images. This effect, labeled “cosmic shear” is arguably the most promising observational probe for photometric redshift surveys, and therefore we will discuss the KL analysis of this signal in particular detail.
    
    As in the case of lensing magnification, and unlike the dominant galaxy clustering terms, the cosmic shear signal is correlated between widely separated redshift bins due to the integrated nature of gravitational lensing. Thus we can expect that a KL transform should be able to compress most of the signal to noise into a small set of radial eigenmodes. To quantify this we consider the same survey configuration used in Section <ref>. The signal part of the cross-power spectrum between the cosmic shear measurements made in two different redshift shells is given again by Eq. <ref>, where now the transfer functions Δ^α_ℓ take the form:
    
    Δ^α_ℓ(k)≡√((ℓ+2)!/(ℓ-2)!)∫ dχ W^α(χ)j_ℓ(kχ)/k^2a(χ)√(P(k,z(χ))),
    
             W^α(χ)≡3H_0^2Ω_M/2∫_z(χ)^∞ dz ϕ^α(z')χ(z')-χ/χ(z')χ.

    The noise power spectrum is white and simply given by the intrinsic ellipticity scatter weighed by the angular number density of sources in each redshift bin n̅^α:
    
    N^αβ_ℓ=δ_αβσ_γ^2/n̅^α,

    with n̅^α in units of srad^-1 and σ_γ=0.28 <cit.>. The lensing auto-power spectra (both signal and noise) for these bins are shown in Figure <ref>.
    
    We compute the KL modes for this setup and rank them according to their contribution to the total lensing signal (in a Fisher matrix sense). The power spectra of the resulting set of modes are shown in Figure <ref>. Comparing against Fig. <ref> we can see that the KL decomposition effectively separates the signal-dominated and noise-dominated modes, with all modes p>3 dominated by noise. The fractional contribution of each mode to the total signal-to-noise is shown in the top panel of Figure <ref>. Most of the signal (∼95%) is contained within a single mode, and the first two modes are able to recover more than 99% of the total. The eigenvectors corresponding to the first three principal modes for different values of ℓ are shown in the right panel of the same figure. We observe that the eigenvectors preserve roughly the same shape for all ℓ, and converge to the same shape at large ℓ. The first eigenvector upweights the parts of the redshift range with the highest signal-to-noise, penalizing the low-z regime due to its poor lensing signal and the high-z bins due to their high shot noise. The second and third eigenmodes then recover part of this information by marginally upweighting these regions. The dashed black line in the same figure shows the weighting scheme associated with the measurement of the lensing signal integrated over a single bin encompassing the redshift range covered by the 16 bins in Fig. <ref>. These weights are similar to the principal KL eigenmode, and thus the KL decomposition determines, broadly speaking, that the bulk of the signal-to-noise is mostly concentrated in the redshift-integrated signal, and extra information regarding the growth of structure can be picked up by up- or down-weighting the contributions in different sections of the redshift range.
    
    

    As we have discussed in the previous sections, the principal KL eigenmodes that optimize the recovery of the cosmological signal are not necessarily optimal in terms of encoding cosmological information, although it is plausible to expect so in general. In order study this further we have performed a Fisher-matrix forecast of the final constraints on cosmological parameters achievable by collecting the information encoded in the first M principal eigenmodes, and compared them with the best possible constraints coming from the use of the full set of 16 redshift bins (or, equivalently, all of the KL eigenmodes). We do so following the approach described in Section 3 of <cit.> and using, as observables, the corresponding set of KL modes b^p_ℓ m. For these forecasts we considered a set of 9 parameters: the relative density of cold dark matter ω_c, the relative contribution of baryons ω_b, the normalized local expansion rate h, the amplitude A_s and spectral index n_s of primordial scalar perturbations, the sum of neutrino masses Σ m_ν, the equation of state of dark energy w and two parameters, log_10M_c and η_b, parametrizing the contribution of baryonic effects in the matter power spectrum as described in <cit.>.
    
    Figure <ref> shows the results of this analysis in terms of 1σ contours in the Σ m_ν-w plane marginalized over all other parameters. The results are shown for the set of 1, 2 and 3 principal KL eigenmodes in red, blue and yellow respectively, while the best achievable constraints using all of the modes are shown as a solid black ellipse. We have removed the axis labels to focus the reader's attention on the relative improvement of the constraints with the number of modes. Even though the first eigenmode contains the vast majority of the lensing signal, as evidenced by the top panel of Fig. <ref>, the extra information contained in the second and third eigenmodes is necessary in order to break the degeneracies between cosmological parameters. Three modes are however sufficient to recover the full constraining power with negligible loss of information.
    
    To finalize, we would like to emphasize the fact that, as shown in the bottom panel of Fig. <ref>, the three principal eigenmodes preserve roughly the same shape as a function of multipole order, converging to the same curve for large ℓ. An ℓ-independent basis of radial functions would be advantageous from the point of view of data analysis since, for instance, the radial window functions associated with each mode (see Eq. <ref>) would only have to be calibrated once (independently of ℓ). It is therefore interesting to explore the cosmological constraints achievable by using the radial functions associated with the KL eigenmodes at high ℓ for all ℓ, even though, for a fixed multipole order, the corresponding set of modes will not be exactly orthogonal. We have verified that, doing so for the first three KL eigenmodes, the final constraints on either w or Σ m_ν degrade by less than 0.5%. This is a reasonable result, given the larger statistical weight of the small-scale (large-ℓ) fluctuations.
    


§ DISCUSSION

  Next-generation cosmological observations will gather their constraining power from a variety of observables, and will therefore have to deal with enormous data vectors. This will present a computational challenge, both from the point of view of likelihood evaluation and in the estimation of the covariance matrix. An efficient data compression scheme would be able to not only alleviate these problems, but also to separate the most significant and less contaminated modes in the data.
  
  In this article we have studied the problem of 3D data compression in the context of photometric redshift surveys, and presented a method, based on the KL transform, to derive a basis of orthogonal radial functions that optimally separate the data into uncorrelated modes, where optimality can be defined in terms of overall signal to noise ratio, information content on a particular cosmological parameter or separability between clean signal and contaminants. This basis is a general and natural extension of the well-know harmonic-Bessel (or Fourier-Bessel) decomposition of spherically-symmetric and translationally-invariant systems, adapted to the particular properties of the dataset under study. Even though the definition of this basis requires prior knowledge of some of these properties, including uncertain ones such as a model for the photo-z distributions, once the radial eigenfunctions are selected, the analysis of the resulting data eigenmodes can proceed as usual, including any calibration of these properties. Thus, a suficiently well-educated model of the survey parameters should preserve the near-optimality of the associated eigenbasis, while not hampering the robustness of the analysis.

  We have shown that, for the study of galaxy clustering in an idealized spectroscopic survey, the optimal set of eigenmodes corresponds to the standard harmonic-Bessel basis, and that this would not be the case in the presence of redshift uncertainties, RSDs or in the analysis of weak lensing observables. For the study of galaxy clustering in a photometric redshift survey, we have shown that the KL basis that maximizes the recovery of the cosmological signal is Bessel-like, although more optimal compression schemes can be derived to optimize the measurement of individual cosmological parameters. In particular, in the case of f_ NL we have shown that the bulk of the constraining power is concentrated in ∼3 radial modes. We have also extended the method to maximize the recovery of a particular signal in the presence of correlated contaminants, and shown that it could be used to simplify the measurement of the effect of lensing magnification as a subdominant contribution to the statistics of the galaxy distribution.
  
  In the case of cosmic shear measurements we have shown that, due to the integrated nature of the gravitational lensing effect, the bulk of the signal (∼95%) is concentrated in a single radial mode, qualitatively equivalent to the measurement of the weak lensing effect over the full redshift range of the survey. The next subdominant modes are however needed in order to break degeneracies between different parameters, and we have shown that three modes are enough to recover the best achievable cosmological constraints.
  
  Further work is needed in order to alleviate some of the practical shortcomings of the method: the KL decomposition is arguably less connected to real-space, directly observable quantities (although not less so than standard Fourier-space methods). Some of the usual methods for systematics calibration thus need to be adapted for a KL-based analysis, and this is particularly relevant for photo-z systematics. In the case of weak lensing, however, we have shown that the shape of the radial eigenfunctions converges to the same curve on large multipole orders, and that the use of ℓ-independent eigenfunctions would have a negligible impact on the final cosmological constraints. In this case, photo-z calibration methods could be used in exactly the same manner as in the standard tomographic analysis.
  
  It is also worth emphasizing that, as is the case for the standard harmonic decomposition of fields defined on the sphere, the KL radial eigenmodes are no longer uncorrelated in the presence of an incomplete sky coverage, and a standard pseudo-C_ℓ analysis reveals non-zero coupling between different multipole orders (ℓ,ℓ')as well as different KL indices (p,p') (see Appendix <ref> and <cit.>). The impact of these correlations on the performance of the KL decomposition should be studied in more detail, and well-understood contaminant-deprojection techniques, implemented in standard power spectrum methods <cit.>, should be adapted for this analysis.
  
  Finally, although we have explored the applicability of this method to independent galaxy clustering and weak lensing measurements, current and upcoming photometric redshift surveys will draw cosmological constraints from a joint analysis of both observables <cit.>. The direct application of this method to the joint data vector would in general produce eigenmodes that mix both signals. Alternatively a joint analysis of the KL modes of each observable, taken individually, could be performed, and the merits and drawbacks of each approach should be studied in detail.  
  


§ ACKNOWLEDGEMENTS

  The author would like to thank Justin Alsing, Pedro Ferreira, Alan Heavens, Boris Leistedt, Jason McEwen, Anže Slosar and David Spergel for useful comments and discussions, and the Center for Computational Astrophysics, part of the Simons Foundation, for their hospitality. He also aknowledges support from the Science and Technology Facilities Council and the Leverhume and Beecroft Trusts. 
  




§ PSEUDO-C_ℓ ESTIMATION OF THE KL MODES

  One of the standard methods to estimate the angular power spectrum of any two quantities in the cut sky is the so-called pseudo-C_ℓ estimator <cit.>. This method can be directly applied to the two-point statistics of the KL eigenmodes, and reveals the correlation between radial modes generated by an incomplete sky coverage <cit.>.
  
  The standard pseudo-C_ℓ method is based on computing the spherical harmonic coefficients of the masked field:
  
    â^α_ℓ m=∫ d a^α() w^α(),

  where w^α is the weights map characterizing the mask of the field a^α. One then estimates the power spectrum of this object by averaging over m for each ℓ:
  
    Ĉ^αβ_ℓ≡∑_mâ^α_ℓ mâ^β *_ℓ m/2ℓ+1.

  This object is then related to the true underlying power spectrum through a mode-coupling matrix M^αβ_ℓℓ' such that
  
    ⟨Ĉ^αβ_ℓ⟩=∑_ℓ'M^αβ_ℓℓ'C^αβ_ℓ',
        M^αβ_ℓℓ'≡∑_ℓ”(2ℓ'+1)(2ℓ”+1)/4πW^αβ_ℓ”(
        [  ℓ ℓ' ℓ”;  0  0  0 ])^2

  where the coupling matrix M depends solely on the power spectrum of the masks W^αβ_ℓ≡(2ℓ+1)^-1∑_mw^α_ℓ mw^β *_ℓ m.
  
  The extension of this estimator to the power spectrum of the KL modes is straightforward: we project the masked harmonic coefficients â^α over the KL eigenvectors F_ℓ (i.e. b̂_ℓ m≡ E_ℓ∘â_ℓ m) and compute their power spectra by averaging over m. The resulting estimator takes the form D̂^p_ℓ=∑_ℓ'M_ℓℓ'^pp'D^p'_ℓ', where the new mode-coupling matrix is given by:
  
    M^pp'_ℓℓ'≡∑_αβM^αβ_ℓℓ'[∑_α'( F_ℓ)^p_α( N^-1)_αα'( F_ℓ')^p'_α'][∑_β'( F_ℓ)^p_β( N^-1)_ββ'( F_ℓ')^p'_β']
        =M_ℓℓ'[∑_αβ( F_ℓ)^p_α( N^-1_ℓ)_αβ( F_ℓ')^p'_β]^2

  where the second equality holds only if all the maps a^α_ℓ share the same mask w. Note that, for full-sky coverage M_ℓℓ'=δ_ℓℓ', and using the orthonormality of F we recover M^pp'_ℓℓ'=δ_ℓℓ'δ_pp'.

