rogerblandford / Music

Speculative project to map the entire universe over time with low spatial resolution using CMB, 21 cm and local surveys
MIT License
0 stars 5 forks source link

SVD tests #21

Open drphilmarshall opened 8 years ago

drphilmarshall commented 8 years ago

The galactic plane area in the CMB T maps has been infilled, and so seems to be leading to an ill-conditioned $a_{lm}$ covariance matrix. So, let's try inverting it with SVD, and setting the singular values below some threshold to zero (which should be equivalent to assigning infinite uncertainty to the infilled pixel values, at least approximately.

1) We decided to start with mock data, and try three different reconstruction tests:

a. Make $\phi$, generate noisy $T$ with no mask, reconstruct $\phi$ using SVD, compare with input $\phi$ - this is a functional test b. Make $\phi$, generate noisy $T$ with no mask, then mask and zero $T$, reconstruct $\phi$ using SVD, compare with input $\phi$. Does SVD enable us to work with infilled maps? Is this independent of the mask chosen? c. Make $\phi$, generate noisy $T$ with no mask, then mask $T$ and put very large variance on the masked pixels (instead of zeros), reconstruct $\phi$, compare with input $\phi$. This will double-check that this is effectively what SVD on the masked. zero'd data is doing.

2) Working with real data, we want to:

a) Start from a masked and zeroed $T$ map, reconstruct $\phi$, make a T replica map.

Notes:

LaurencePeanuts commented 8 years ago

Hi @drphilmarshall !

I think 1 a) works pretty well now, check it out here. It includes multipoles from $l_min =1$ to $l_max=10$, which corresponds to $n_min=1$ to $nmax=3$ in the 3d box. I needed to use SVD both for inverting the covariance matrix $C{yy}$ and the matrix $A=(R^T C^{-1}_{yy} R + C_f^{-1})$ , but with this it seems to be behaving very well.

I can check how it behaves as I change $l_max$, and I can calculate the figure of merit from eq. 10 in the allegro draft, but I haven't checked what sort of of numbers I should expect, although from looking at it I'm guessing it should be $\mathcal{O}(1)$.

Let me know when you want to chat!

Cheers, Laurence

drphilmarshall commented 8 years ago

Excellent! Nice work Laurence :-) I'm at SLAC tomorrow if you want to chat. It'd be good to walk through your nice demo, I have a few questions and some ideas for making it clearer.

I see the residuals are about 10% of the map - is that to be expected, given the amount of noise you put in? Or are we seeing some model leakage? Either way, it could be nice to show the T maps all on the same fixed scale, -0.8 to +0.8 perhaps - to give a visual impression of the goodness of fit. What we really ought to do is plot the residuals scaled by the uncertainty map - which we could get by drawing 100 sample potential models from the posterior PDF for the fn, making the T map for each one, and building up map pixel statistics. It won't be perfect, as the standard deviation map doesn't capture the covariance between predicted map pixels, but it's something. The animated GIF of the 100 samples probably gives a better impression of whether our uncertainties are accurate. Is it very expensive to compute the covariance matrix of the Gaussian posterior PDF for the fn?

On Thu, Jan 21, 2016 at 3:45 PM, Laurence Perreault Levasseur < notifications@github.com> wrote:

Hi @drphilmarshall https://github.com/drphilmarshall !

I think 1 a) works pretty well now, check it out here https://github.com/rogerblandford/Music/blob/master/Reconstruction_Demo.ipynb. It includes multipoles from $l_min =1$ to $l_max=10$, which corresponds to $n_min=1$ to $nmax=3$ in the 3d box. I needed to use SVD both for inverting the covariance matrix $C{yy}$ and the matrix $A=(R^T C^{-1}_{yy} R + C_f^{-1})$ , but with this it seems to be behaving very well.

I can check how it behaves as I change $l_max$, and I can calculate the figure of merit from eq. 10 in the allegro draft https://github.com/rogerblandford/Music/blob/master/doc/music-allegro.pdf, but I haven't checked what sort of of numbers I should expect, although from looking at it I'm guessing it should be $\mathcal{O}(1)$.

Let me know when you want to chat!

Cheers, Laurence

— Reply to this email directly or view it on GitHub https://github.com/rogerblandford/Music/issues/21#issuecomment-173752051 .

drphilmarshall commented 8 years ago

PS. Can you make 12:30pm tomorrow?

On Thu, Jan 21, 2016 at 4:11 PM, Phil Marshall dr.phil.marshall@gmail.com wrote:

Excellent! Nice work Laurence :-) I'm at SLAC tomorrow if you want to chat. It'd be good to walk through your nice demo, I have a few questions and some ideas for making it clearer.

I see the residuals are about 10% of the map - is that to be expected, given the amount of noise you put in? Or are we seeing some model leakage? Either way, it could be nice to show the T maps all on the same fixed scale, -0.8 to +0.8 perhaps - to give a visual impression of the goodness of fit. What we really ought to do is plot the residuals scaled by the uncertainty map - which we could get by drawing 100 sample potential models from the posterior PDF for the fn, making the T map for each one, and building up map pixel statistics. It won't be perfect, as the standard deviation map doesn't capture the covariance between predicted map pixels, but it's something. The animated GIF of the 100 samples probably gives a better impression of whether our uncertainties are accurate. Is it very expensive to compute the covariance matrix of the Gaussian posterior PDF for the fn?

On Thu, Jan 21, 2016 at 3:45 PM, Laurence Perreault Levasseur < notifications@github.com> wrote:

Hi @drphilmarshall https://github.com/drphilmarshall !

I think 1 a) works pretty well now, check it out here https://github.com/rogerblandford/Music/blob/master/Reconstruction_Demo.ipynb. It includes multipoles from $l_min =1$ to $l_max=10$, which corresponds to $n_min=1$ to $nmax=3$ in the 3d box. I needed to use SVD both for inverting the covariance matrix $C{yy}$ and the matrix $A=(R^T C^{-1}_{yy} R + C_f^{-1})$ , but with this it seems to be behaving very well.

I can check how it behaves as I change $l_max$, and I can calculate the figure of merit from eq. 10 in the allegro draft https://github.com/rogerblandford/Music/blob/master/doc/music-allegro.pdf, but I haven't checked what sort of of numbers I should expect, although from looking at it I'm guessing it should be $\mathcal{O}(1)$.

Let me know when you want to chat!

Cheers, Laurence

— Reply to this email directly or view it on GitHub https://github.com/rogerblandford/Music/issues/21#issuecomment-173752051 .

LaurencePeanuts commented 8 years ago

12:30 works for me, and Roger says he'd like to join in later in the afternoon, around 1:30 ish. From my realizations of the noise, it looks like it's pretty big, i.e. 1% to 100% of the individual a_y values . There might be some units problems here (i.e. the covariance matrix for the noise is possibly in temperature perturbation, but i'm doing everything in potential perturbations, so that gives a factor of 3^2 right there). Without noise I get the map with lmax=10 back within 6-7% lmax=30 back within 0.2% . We should talk more tomorrow :)

Best, Laurence

drphilmarshall commented 8 years ago

Fantastic! Sounds like we have a functioning system :-)

On Thu, Jan 21, 2016 at 5:39 PM, Laurence Perreault Levasseur < notifications@github.com> wrote:

12:30 works for me, and Roger says he'd like to join in later in the afternoon, around 1:30 ish. From my realizations of the noise, it looks like it's pretty big, i.e. 1% to 100% of the individual a_y values . There might be some units problems here (i.e. the covariance matrix for the noise is possibly in temperature perturbation, but i'm doing everything in potential perturbations, so that gives a factor of 3^2 right there). Without noise I get the map with lmax=10 back within 6-7% lmax=30 back within 0.2% . We should talk more tomorrow :)

Best, Laurence

— Reply to this email directly or view it on GitHub https://github.com/rogerblandford/Music/issues/21#issuecomment-173772276 .