1.  Analytical problem :

equation to solve for {$ \mathbf E_\ell $} {$\mathbf C^{AA} \mathbf E_\ell \mathbf C^{BB} + \mathbf C^{AB} \mathbf E_\ell^T \mathbf C^{AB} = \mathbf P_{\ell} $}

of when expended :

{$ 2 \mathbf S \mathbf E_\ell \mathbf S + \mathbf N^{AA} \mathbf E_\ell \mathbf S + \mathbf S \mathbf E_\ell \mathbf N^{BB} + \mathbf N^{AA} \mathbf E_\ell \mathbf N^{BB} = \mathbf P_{\ell} $}

1.1  Approximation :

At this stage, we must approximate the shape of the pixel covariance matrices. We consider two extreme cases :

• Low signal-to-noise ratio : {$ \mathbf S \ll \mathbf S + \mathbf N$} such that {$\mathbf C^{AB} \ll \mathbf C^{AA} \sim \mathbf C^{BB} $}
• High signal-to-noise ratio : {$ \mathbf N \ll \mathbf S + \mathbf N$} such that {$\mathbf C^{AB} \sim \mathbf C^{AA} \sim \mathbf C^{BB} $}.

Both limits give a solution of the form

{$ \mathbf E_\ell \simeq \frac 1a \sum_{\ell' } B_{\ell\ell'} (\mathbf C^{AA})^{-1} \mathbf P_{\ell'} (\mathbf C^{BB})^{-1} $}

We will come back to this assumption in the simulations analysis. Where {$1 \leq a \leq 2$} is a normalization coefficient that depends on the signal-to-noise level as discussed above. We absorb this factor such that {$B_{\ell\ell'} = a \tilde B_{\ell\ell'}$}. Finally, imposing {$ Tr [\mathbf E_\ell \mathbf P_{\ell'} ] = \delta_{\ell\ell'}$} gives

{$ Tr[ \sum_{\ell'' } \tilde B_{\ell\ell''} (\mathbf C^{AA})^{-1} \mathbf P_{\ell''} (\mathbf C^{BB})^{-1} \mathbf P_{\ell'} ] = \delta_{\ell\ell'} $} {$ \Leftrightarrow ~ 2 \mathcal B \cdot \mathcal F = \mathbb 1 $} {$ \Rightarrow ~ \tilde {\mathcal B} = \frac1{2} \mathcal F^{-1} $}

where we define

{$ F_{\ell \ell'}^{AB} \equiv \frac 12 Tr[(\mathbf C^{AA})^{-1} \mathbf P_{\ell} (\mathbf C^{BB})^{-1} \mathbf P_{\ell'}] $}

as the 'cross' Fisher information matrix (see {\color{blue} cite}).

We get the fully written spectrum estimator

{$ \hat C_\ell^{AB} = \sum_{\ell'} [ F^{-1}] _{\ell \ell'} \mathbf d^A ( \mathbf C^{AA})^{-1} \mathbf P_{\ell'} ( \mathbf C^{BB})^{-1} \mathbf d^B. $}

1.2  Exacte solution of the Sylvester equation :

• proprety : {$ vec[ACB] = (A \otimes B) vec[X] $} where {$vec[X]$} is the vectorization of matrix {$X$} (meaning joining successively all its columns, starting from the first), and {$\otimes$} is the kronecker product such that {$ (A \otimes B)_{ij} = A_{ij} B$} (see https://en.wikipedia.org/wiki/Kronecker_product .

A matrix equation of the form {$ A E B = P$}, where {$E$} is the unknown, can be linearly transformed into the form {$ (A \otimes B) e = p$}. And {$e$} is the vectorization of the matrix {$E$} , same for {$p$}.

The expended equation above thus becomes

{$ ( 2 \mathbf S \otimes \mathbf S + \mathbf N^{AA} \otimes \mathbf S + \mathbf S\otimes \mathbf N^{BB} + \mathbf N^{AA} \otimes \mathbf N^{BB} ) e_\ell = p_{\ell} $}

and can be solved exactly. Also known as the Sylvester equation. Using a linear solver algorithm (rather than inversing the matrix).

• Note :
  • the kronecker product between a {$ m \times n $} matrix and a {$ p \times q $} is a matrix of dimension {$ mp \times nq $}. This means that {$n_{pix}$} square matrices, we have to solve a square matrices system of dimension {$n_{pix}^2$} !

2.  Results :

In the xQML code, the there is no normalization factor {$ B_{\ell\ell''} $} or {$a$} in {$E_\ell$}, since those cancel out in the end, for the final spectrum estimator and analytical variance estimate. We thus simply compute

{$E_\ell = (\mathbf C^{AA})^{-1} \mathbf P_{\ell} (\mathbf C^{BB})^{-1} $}

and

{$F_{\ell\ell'} = Tr[ E_\ell P_{\ell'} ] $}

Approximation Vs Exacte

El matrix

• nside= 4, thus we need a very high level of noise to be in the low signal-to-noise regime :
• three we to solve :
  • Approximation
  • exacte :
    • inverse solver : {$ M x = b $} => {$ x = M^{-1} b $}
    • linear numpy solver : solve {$ M x = b $} for {$ x $}
• We compared the diagonal of {$ \frac 12 E_0 $} from the approximation and the exacte solution {$ E_0 $}

muK =	0.1	1.0	5.0	10.0	20.0	50.0	100.0	200.0	500,0	1000.0
inv solving
numpy solving

• observations :
  • using the solver of numpy is safer given the matrix conditionning due to the low noise level.
  • The approximation and exacte solution results are of the same order for low noise, as expected. ( since El/2 means {$ E_\ell/2 \simeq S^{-1} P_\ell S^{-1} /2 $} )
  • The approximation and exacte solution results are of differ from a factor 1/2 for high noise, as expected.
  • Since we know that there is a factor of difference (btw 1 and 2 ) between the approximation and the exact solution, we cannot compare the {$ E_\ell $} matrices directly. But we can compare the impact on the spectrum variance.

Spectrum estimator and variance impact :

2.1  LiB :

We plot the estimator spectrum errors ratio (minus one) between the approximation and the exact solution for {$ E_\ell $}, using a wide range of noise level (between 0.1 muK to 2000.0 muK). The ratio is computed from the diagonal of the covariance matrix (EE ++ BB )

• remark : -1 is missing in the y label

nside = 2, fsky = 1.0
nside = 2, fsky = 0.5

• observation :
  • For extreme cases (dominant signal or dominant noise), the ratio is very low, meaning that the approximation is very close to the exact solution of {$ E_\ell$}, as expected.
  • for intermediate cases, around 20.0muK, the ratio peaks. It can reach around 1% for fsky=0.5, which is low.
  • The impact of the approximation seems higher for larger masked fraction of the sky.
  • For nside = 2, fsky = 1.0, the 0.1 muK case seems unstable seems the ratio is negative, invisible on the log scale. exept last bin of BB wich is positive (taille of the ratio curve).

2.2  Bicep :

Same analysis for bicep :

nside = 16, fsky = 0.01

• observation :
  • Peak at 1%, low.
• => Analysis should be done at higher resolution ? Computionnaly costly ...
• => using numerical algorithms ?