LSSTDESC
Expansion of PhotoZ uncertainty model
New feature requests and enhancements
I want to add a new photometric redshift uncertainty model that represents the uncertainty on the n(z) as an expansion in terms of a provided family of functions.
Problem description
If we are to reproduced the measured variability in the galaxy distributions inside
Firecrown better than order 1, we will have to include more flexible photometric redshifts uncertainty models in the library such as Gaussian processes (GP) and Principal Components (PCA).
Proposed solution
Function Expansion & PCA
Let's start writing a redshift galaxy distribution as $$n(z) = \bar{n}(z) + \delta n(z) .$$ where $ \delta n(z) $ is the uncertainty we want to capture. Then we can write $\delta n(z)$ as a expansion like: $$\delta n(z) = \sum^N_i w_i \phi_i , . $$ Plugging this back into the previous expression we can write a generative model for the $n(z)$: $$n(z) = \bar{n}(z) + \sum^N_i w_i \phi_i = \bar{n}(z) + \left< w_i \phi_i \right> .$$
Within this framework, we $w_i$ are the free parameters of the inference and $\phi_i$ are set of vectors provided by the user and stored in the class. In the case of PCA $\phi_i$ are the eigen-vectors of ensemble of $n(z)$ residuals but they can technically be any family of functions. The parameters $w_i$ might have a non-trivial prior but we can always sample them from standard Gaussian. Let us consider the individual weights of the functions as a vector $\boldsymbol{w} \in \mathbb{R}^N$. Then, a sample of $\boldsymbol{w}_j$ can be generated as follows. Consider a set of latent variables $\boldsymbol{\alpha} \in \mathbb{R}^N$ distributed as $\boldsymbol{\alpha} \sim \mathcal{N}(0, \mathbb{I})$. Then a sample of $\boldsymbol{w}_j$ can be generated as:
$$\boldsymbol{w}_j = \bar{\boldsymbol{w}} + \mathbb{L} \boldsymbol{\alpha}_j ,$$
where $\bar{\boldsymbol{w}}$ is the mean value of the weights and $\mathbb{L}$ is the lower Cholesky triangle of the covariance matrix of $\boldsymbol{w}$ which are fixed quantities stored within the class. Plugging this back into our generative model, a sample of the $n(z)$ as:
$$n(z)_j = \bar{n}(z) + \sum^N_i (\bar{\boldsymbol{w}} + \mathbb{L} \boldsymbol{\alpha}_j)_i \phi_i .$$
Since we are performing an expansion on the residues of the $n(z)$, $\bar{w} \approx 0$ such that
$$n(z)_j = \bar{n}(z) + \sum^N_i (\mathbb{L} \boldsymbol{\alpha}_j)_i \phi_i , .$$
Therefore we can write $\psi = \mathbb{L}^T\phi$ such that:
$$n(z)_j = \bar{n}(z) + \sum^N_i \psi_i (\alpha_i)_j = \bar{n}(z) + \left< \psi, \alpha_j \right> , $$
where $\psi$ are a set of matrices that we can save inside the class.
GP
Gaussian processes (GP's) can appear as a different approach from function space expansions but they are mathematically identical. We can write a GP for $n(z)$ as: $$n(z) = \bar{n}(z) + \left< \mathbb{L} , \alpha \right> , $$ where $\mathbb{L}$ is the lower Cholesky triangle of the covariance of the $n(z)$ residuals. However, this approach requires $\alpha$ to have the same dimensions as $n(z)$ which is undesirable. In order to solve this we can define $\alpha$ on a latent space $\mathbb{R}^N$ where N is the number of desired free parameters. Therefore, we can write $\delta n(q) = \left< \mathbb{L} \alpha \right>$ where $q \in \mathbb{R}^N$ is a redshift array with lower resolution than $z$. In order to map back to $z$ we simply have to compute the Wiener filter from $q$ to $z$ and apply it to $\delta n(q)$. Thus: $$\delta n(z) = \mathbb{W} \delta n(q)$$ where $$\mathbb{W} = \frac{\mathbb{C}_qz}{\mathbb{C}_qq}$$ .
plugging this back in our original expression we can obtain a generative model for the $n(z)$:
$$n(z)_j = \bar{n}(z) + \mathbb{W} \left< \mathbb{L} , \alpha_j \right> = \bar{n}(z) + \left< \psi , \alpha_j \right>$$ ,
where we have grouped $\psi = \mathbb{W} \mathbb{L}$. Therefore we can see that one can implement both GP's and function space uncertainty models by storing the $psi$ matrices inside the class .
Questions:
- Can
FireCrownstore vector variables of arbitrary length? For example a vector variable that only gets an specific dimension once the $\psi$ matrices are provided?
@marcpaterno @vitenti any chance you could have a look at how to go about implementing this?
@JaimeRZP I believe the implementation is fairly straightforward. Would you be able to join us next week to start drafting it together?
