text:
\section{Introduction}
\label{sec:introduction}

The state of the inhomogeneous cosmic matter distribution (and the time-dependent
dark energy density) can be measured by
the gravitational lensing shearing and magnification of background luminous
sources. Using galaxies as the back lights for gravitational lensing (``cosmic shear'')
allows tomographic reconstruction of the matter distribution along the line of sight.
But, inference of galaxy lensing shears is clouded by the unknown distribution
(and possibly correlation) of galaxy intrinsic morphologies as well as image systematics
that induce spurious ellipticity correlations.

Common inference algorithms for gravitational lensing shear involve cross-correlating
the ellipticities of galaxies in a two-point function estimator under the assumption
of calibrated intrinsic ellipticity distributions and alignments~\citep[e.g.,][]{jee2015}.
This is a lossy procedure because the angular phase information in the lensing
shear and magnification fields on the sky is discarded. Traditional algorithms
are also necessarily biased because of the need to calibrate the unknown galaxy
ellipticity distributions. Said another way, the two-point function summary statistics of
cosmological large-scale structure do not capture the full statistical information in the multivariate
distribution of the lensing observable
\citep[e.g.,][]{pan2005,hamilton2006,takada2013,carron2014,2016arXiv160501100P}

In this paper we revisit the problem of inferring the lensing convergence posterior distribution
from a catalog of galaxy ellipticities to potentially capture more cosmological information from
cosmic shear than is available in two-point function estimators of the shear. We pursue a
probabilistic approach to lensing convergence inference as a means to propagate the large
uncertainties that can come in the presence of galaxy shape noise, finite survey areas, masking,
and sample selection. Moreover, we aim to build an inference framework that can fit within our previous
work on probabilistic cosmic shear~\citep{mbi-theory}.

\citet{1990ApJ...349L...1T} first presented a method to reconstruct the projected mass distribution
from weak gravitational lensing shears of source galaxy images.
\citet{1993ApJ...404..441K} derived the lensing convergence in terms
of the shear (a nonlocal relation) and applied this theory to estimators of the convergence given a
measured galaxy ellipticity catalog. The method of \citet{1993ApJ...404..441K}, while theoretically sound,
requires unbroken sky coverage and low-pass filtering to yield finite noise in the lensing
convergence estimator.
The \citet{1993ApJ...404..441K} method remains useful, however, for visualization
purposes~\citep{2015PhRvL.115e1301C,2015PhRvD..92b2006V}.
Subsequent papers have extended, applied, and explored the limitations of filtering algorithms for
cluster mass mapping~\citep{1995A&A...303..643B,1995A&A...294..411S,1998A&A...335....1L,2002ApJ...568..141G,2005A&A...440..453D,2007Natur.445..286M,2012A&A...540A..34D,2012MNRAS.424..553A,2013MNRAS.433.3373V}. Remaining data analysis challenges
using such methods include finite survey boundaries and masks,
separation of E and B modes in the shear
field, noise or significance characterization for shear or convergence extrema, and the requirement
in many algorithms to smooth or average the ellipticities of galaxies before the convergence inference
process.

Maximum likelihood (ML) estimators for the lensing
convergence~\citep[e.g.,][]{1996ApJ...464L.115B,1998A&A...337..325S,2008MNRAS.385.1431H,2009ApJ...702..980K}
can help mitigate biases arising from survey masks and admit mathematically consistent noise
characterization~\citep{2000MNRAS.313..524V}. However such estimators are biased. And many
approaches often still require a preliminary smoothing of the observed galaxy ellipticity field.

\citet{bridle1998} and \citet{bridle2000} introduced a `maximum entropy'
Bayesian prior for the lensing convergence from information theoretic and Bayesian analysis perspectives to
derive an estimator for the projected mass distribution of galaxy clusters with desirable noise
properties in a finite field.
\citet{Marshall:531310} and \citet{2001A&A...368..730S} refined the maximum entropy method to
specify the optimal smoothing
length scale for galaxy cluster mass inference via the Bayesian evidence of the observed
ellipticities. \citet{2013ascl.soft08004M} released a code\footnote{{\sc LensEnt2}\xspace, \url{http://www.slac.stanford.edu/~pjm/lensent/version2/index.html}}
implementing the algorithm from \citet{Marshall:531310}.
The algorithm of \citet{Marshall:531310} is close to meeting all the requirements for the current
analysis, except the choice of smoothing scale and application to field rather than cluster lensing
is not demonstrated in the literature. We will show further benefits of the algorithm developed
in this paper below.
\citet{2006A&A...451.1139S} applied a modification of the maximum entropy method to a cosmic
shear analysis using $N$-body cosmological simulations to create mock observations.
The work of \citet{2006A&A...451.1139S} and also \citet{jiao2011} include maximum entropy algorithm
modifications to better handle masking, shape noise, and degrees of smoothing.

For constraining cosmological parameters, several groups have considered the abundance of peaks
in lensing convergence maps~\citep{2010PhRvD..81d3519K,2011ApJ...735..119S,2014MNRAS.442.2534S,2015PhRvD..91f3507L,2016MNRAS.456..641R,2016PhRvD..94d3533L,PhysRevLett.117.051101,2016arXiv160501100P}. These studies have a common
approach in direct calculation or estimation of lensing peaks without an attempt to infer
associated peaks in the 3D cosmological mass density. The measurement process is thus direct as
lensing by multiple structures or voids along any given line of sight can easily confuse the
2D to 3D mass inference. Just as the abundance of galaxy clusters provides tight constraints on the
cosmological model, so too can the abundance of lensing peaks.
However, the estimators for lensing peaks often involve averages or line integrals over contiguous
sky areas, which are confounded by survey masks~\citep{0004-637X-784-1-31}.
The bias introduced by masking may be overcome with forward simulations of the nonlinear lensing
convergence field, which has been achieved with $N$-body cosmological simulations of large-scale
structure combined with ray-tracing predictions of the lensing statistics~\citep[e.g.,][]{2016ApJ...819..158B,2016arXiv160501100P}.

In a complimentary approach to cosmological parameter inference, other groups have considered a
linear theory approximation to (suitably smoothed) lensing mass maps.
Specifically, if (i) the shear field can be approximated as Gaussian distributed, (ii) the noise
distribution in the galaxy ellipticity measurements (and intrinsic shape distribution) can also
be approximated as Gaussian distributed, and (iii) the weak shear approximation is valid, then
mass maps can be obtained using the Wiener Filter (WF)~\citep{wiener1949extrapolation}.
The WF is a well-studied signal-inference algorithm that has been applied to cosmic microwave
background (CMB) analysis with great success in a Bayesian context~\citep{2004PhRvD..70h3511W}.
\citet{Alsing2015} and \citet{2016arXiv160700008A} have recently extended the WF technology from
the CMB literature to the problem of weak lensing shear field and power spectrum inference. However,
these WF approaches are limited not only by the WF assumptions above but also by requirements
to smooth the measured galaxy ellipticity field to a uniform grid on the sky and to work solely
within a linear theory approximation for the distribution of cosmological mass density perturbations.
Working within the stated assumptions, the primary challenge for WF approaches is the computation
of large matrix solve operations. Novel and effective algorithms for sampling from the WF
distribution have been demonstrated for the CMB~\citep{2013A&A...549A.111E,2016ApJ...820...31R}
that are also effective for smoothed, weak-shear, and linear theory cosmic shear
inferences~\citep{2016arXiv160700008A}.

This paper is structured as follows. In \autoref{sec:method} we describe the statistical
framework for lensing convergence and shear inference given galaxy images or a galaxy
ellipticity catalog. In \autoref{sec:framework} we describe the Gaussian Process (GP)
prior for the lensing potential and how this informs the correlated inferences of shear
and convergence. We give a prescription for GP parameter optimization in
\autoref{sub:optimizing_interpolation_parameters}.
We apply our method to infer lensing convergence maps of simulated ellipticity catalogs
in \autoref{sub:simulation_study} and of an observed galaxy cluster in \autoref{sub:abell}.
We describe our main conclusions in \autoref{sec:conclusions}.
We provide details of the GP covariance derivation for lensing shear and convergence
fields in \autoref{sec:Gaussian process covariances} and provide
lensing map inferences under an analogous covariance
for a cosmological model with a linear theory approximation in
\autoref{sec:cosmology_dependent_covariance_model}.

\section{Method}
\label{sec:method}

We describe a joint probabilistic model for the lensing shear and convergence given galaxy imaging
data. We previously presented the complete statistical framework for cosmic shear
inference~\citep{mbi-theory} and here focus on the specifics of lensing convergence inference
as a function of sky coordinates in both linear and nonlinear regimes of the cosmological
mass density perturbations. In \autoref{fig:pgm} (left panel) we show the relationships in
the probabilistic model for CCD galaxy imaging data for multiple observation epochs, multiple
galaxies, and multiple galaxy samples (e.g., samples selected in different photometric redshift bins).
\begin{figure*}[!htb]
  \centerline{
    \includegraphics[width=0.48\textwidth]{shear_general_pgm.png}
    \includegraphics[width=0.4\textwidth]{shear_gp_pgm.png}
  }
  \caption{Two successive levels of approximation for our statistical model for sampling
  probabilistic lensing shear and convergence $\Upsilon$ fields.
  The unshaded circles indicate sampling parameters.
  Shaded circles indicate observed parameters while dots indicate parameters with fixed values
  rather than being sampled.
  %
  Left panel: The model for galaxy image pixel data ${\mathbf{d}}_{n,i}$ for each galaxy $n$ and observation
  epoch $i$ requires specification of the pixel noise $\sigma_{n,i}$, the intrinsic (i.e., unlensed)
  galaxy ellipticity $e_n$, the lensing fields $\Upsilon$, and the PSF at each galaxy location in
  each epoch ${\Pi}_{n,i}$. The distribution of galaxy intrinsic ellipticities is described by the
  parameters $\alpha_n$, which can specify distinct distributions for each galaxy $n$.
  In \citet{mbi-theory} we infer marginal constraints on $\alpha_n$ under a Dirichlet Process (DP)
  prior. For this initial study, we assert a Gaussian Process (GP) prior on the lens fields
  $\Upsilon$ in this paper, with the assumption that the posterior inferences of $\Upsilon$
  will be related to a cosmological model in a separate analysis pipeline as we describe in the
  text. Here we assert a known PSF ${\Pi}_{n,i}$ at every galaxy location, again deferring the
  inference and marginalization of PSFs to a separate paper. The inference of the lens fields also
  depends on our assertion of the survey window function $W$.
  %
  Right panel: The approximate statistical model when the galaxy imaging pixel data is summarized as
  a galaxy ellipticity catalog ($e_{n}, \sigma_{e,n}$ for $n=1,\dots,n_{\rm gal}$) requires specification of
  the intrinsic ellipticity distribution parameters $\alpha$ (now assumed the same for all galaxies),
  the ellipticity measurement uncertainties $\sigma_{\rm ms,n}$ per galaxy, and the lens fields
  $\Upsilon$.
  We now also assert the variance, $\alpha$ of the assumed zero-mean
  intrinsic ellipticity distribution.
  }
  \label{fig:pgm}
\end{figure*}

\begin{table*
\begin{center}
\caption{Parameters for the statistical model. }
\label{tab:sampling_parameters}
\begin{tabular}{cl}
\hline
Parameter & Description \\
\hline
$i$ & index over epochs, or different exposures of a galaxy\\
$n$ & index over galaxies\\
${\Pi}_{n,i}$ & Point Spread Function (PSF) for galaxy $n_s$ in epoch $i$ \\
$e_{n}$ & Intrinsic (pre-lensing) ellipticity of a galaxy\\
$\Upsilon_{s}$ & lensing shear and convergence that modify a galaxy image\\
$\psi_{s}$ & lensing potential \\
$\theta$ & cosmological parameters\\
$W$ & survey window function \\
$\sigma_{n,i}$ & noise properties in a galaxy image\\
$\alpha_{n}$ & parameters specifying the intrinsic ellipticity distribution for a galaxy\\
$a_{\rm DP}$ & parameters specifying the distribution over $\alpha_{n_s}$\\
${\bf d_{n}}$ & data vector (measured $e_{1,2}$ for each galaxy $n$)\\
$\sigma_{{\rm ms}; n}$ & ellipticity measurement error for galaxy $n$ \\
$\mathbf{x}$ & vector of 2D spatial locations of galaxies \\
$\epsilon^\mathrm{int}_{n}$ & intrinsic galaxy ellipticity for galaxy $n$ \\
$\alpha\equiv\sigma_{e}^2$ & parameters of the distribution of galaxy parameters \\
W & window function for the survey footprint \\
$a_{\rm GP}$ & parameters of the GP kernel\\
$\lambda_{\rm GP}$ & precision of the GP kernel (element of $a_{\rm GP}$) \\
$\ell^{2}_{\rm GP}$ & squared GP correlation length (element of $a_{\rm GP}$) \\
\hline
\end{tabular}
\end{center}
\end{table*}

The model for observed galaxy images requires (at least) specification of the point-spread
function (PSF), ${\Pi}_{n_s,i}$, at all $n_s$ galaxy locations in all epochs $i$, the intrinsic
or pre-lensing shape or morphology of the galaxy,
$e_{n_s}$, the applied lensing shear and magnification, $\Upsilon_{s}$, and the noise properties
of the pixelated image, $\sigma_{n_s,i}$ (e.g., Gaussian with an asserted r.m.s.).

To simplify the discussion for this paper focused on shear inference, we assume the PSF ${\Pi}_{n_s,i}$ is
known at the locations of all galaxies in all epochs. We do not further consider errors in the
PSF inference in this work. Following \citet{mbi-theory} we allow for the intrinsic ellipticity
(and potentially size, fluxes, etc.) of each galaxy to be drawn from galaxy-specific distributions
with parameters $\alpha_{n_s}$. The parameters $\alpha_{n_s}$ in turn are hierarchically
distributed under a distribution with with parameters $a_{\rm DP}$, allowing
inference of the effective `shape noise' in the
shear inference. We use the label $\Upsilon$ for
the model for the lensing shear $\gamma$ and convergence $\kappa$, which uniquely specify the
components of the trace-free linear distortion matrix  under a weak lensing ($\kappa\ll1$)
approximation of the lens equation~\citep[e.g.,][]{bartelmann01}.

As illustrated by the `plates' in the left panel of \autoref{fig:pgm}, the PSF is unique to each
observation epoch
but the intrinsic properties of a galaxy image are common across epochs (for each unique galaxy).
A key feature in the graphical model of the left panel in \autoref{fig:pgm} is that the lens fields
$\Upsilon_{s}$ are common across all galaxy images in all epochs, indicating that the we require
a spatially correlated model for the coherent lensing shear and magnification patterns on the sky.
Specifying an appropriate correlated model for the lens fields is the main focus of this paper.
We allow the lens field models to be distinct for different galaxy samples $s$ in
\autoref{fig:pgm} because, e.g., galaxies in different photometric redshift bins will be lensed
by partially different foreground mass distributions. Only at the top level of the graphical
model do the lens field inferences from different galaxy samples become connected under a
cosmological model with parameters $\theta$.

We therefore explore an interim probabilistic model for the lens fields inferred from distinct
galaxy samples $s$ such that inferences of $\Upsilon_s$ are statistically independent for
different $s$. This will allow us to separate computationally expensive components of a cosmic
shear inference pipeline, explore multiple cosmological and systematics models for the data,
and eventually perform rigorous uncertainty propagation and marginalization of image and intrinsic
galaxy nuisance parameters as outlined in \citet{mbi-theory}.

In measuring lensing shear of galaxy images, we must marginalize over the
intrinsic ellipticities, $\epsilon^\mathrm{int}$, of the images. This is often done by averaging the ellipticities
of galaxies in neighboring regions of the sky and redshift, where the weights can include the
measurement and shape noise models. But, here we perform a
marginalization over an explicit intrinsic ellipticity distribution,
\begin{multline}\label{eq:marg_like}
  \prf{{\mathbf{d}} | \Upsilon(\mathbf{x}), \alpha} = \prod_{n=1}^{n_{\rm gal}}
  \int d^{m}\epsilon^\mathrm{int}_n\,
  \prf{{\mathbf{d}}_n | \epsilon^\mathrm{int}_n, \Upsilon(\mathbf{x}_n)}
  \\ \times
  \prf{\epsilon^\mathrm{int}_n|\alpha},
\end{multline}
where we assume the likelihood functions for each galaxy image $n=1,\dots,n_{\rm gal}$ are
statistically independent~\citep[see][for more discussion of this assumption]{mbi-theory}.
The data vector ${\mathbf{d}}_{n}$ is composed of either the pixel values contributing to the image of
galaxy $n$ or a summary statistic of those pixel values.

As a pedagogical step in the development of our probabilistic lens field model,
we will consider an approximate likelihood function for summary statistics of the pixel data;
namely estimators for the
ellipticity of each galaxy image $e\equiv e_1 + i e_2$ along with an associated measurement
error per ellipticity component $\sigma_{\rm ms}$. See the right panel of \autoref{fig:pgm}.
We will develop the approximate model for ellipticity measurements as a data vector
in this paper but advocate for the more complete algorithm of \citet{mbi-theory}
for any data analysis because of the known large biases in using ellipticity estimators for
cosmic shear~\citep[e.g.,][]{Refregier++2012,Kacprzak++2012}.

Assuming Gaussian distributed ellipticity measurement errors the likelihood
function is then,
\begin{equation}\label{eq:likelihood}
  \prf{{\mathbf{d}}_n|e_n, \Upsilon(\mathbf{x}_n)} = \mathcal{N}_{\hat{e}_n}
  \left(e_n, \sigma^2_{{\rm ms}; n}\mathbb{1}} % identity matrix: \usepackage{bbold_{2}\right),
\end{equation}
where we explicitly label the data ${\mathbf{d}}_n\equiv\hat{e}_{n}$ as ellipticity estimators,
and the distribution is bivariate given the two ellipticity components.

A general likelihood function depends on the observable galaxy properties
such as ellipticity, size, and flux, which are modified by lensing from the intrinsic
properties described by $\epsilon^\mathrm{int}$.
So, to evaluate \autoref{eq:likelihood}, we define the lensed galaxy parameters,
\begin{equation}\label{eq:lens_transform}
  e_n(\Upsilon(\mathbf{x}_n)) \equiv f(\epsilon^\mathrm{int}_n, \Upsilon(\mathbf{x}_n)),
\end{equation}
where $f(\cdot)$ denotes the function that transforms the intrinsic, unlensed, galaxy
ellipticities $\epsilon^\mathrm{int}_n$ to those in the lensed model $e_n$
under the action of the lensing convergence and shear specified by $\Upsilon(\mathbf{x}_n)$
at the galaxy sky location $\mathbf{x}_n$.
In the weak shear limit defined by $\kappa \ll 1$, \autoref{eq:lens_transform} reduces to
\begin{equation}\label{eq:weak_shear}
  \tilde{e}_{n}^{\rm weak-shear} = \epsilon^\mathrm{int}_{n} + g(\mathbf{x}_n),
\end{equation}
where $g\equiv\gamma/(1-\kappa)\approx\gamma$ is the reduced shear.

For our pedagogy, we specify a Gaussian distribution for the \emph{unlensed}
galaxy properties $\epsilon^\mathrm{int}$ to use in evaluating \autoref{eq:marg_like},
\begin{equation}
  \prf{\epsilon^\mathrm{int}_{n}|\alpha,\Upsilon(\mathbf{x}_n)} =
  \mathcal{N}_{\epsilon^\mathrm{int}_n} \left(0, \sigma_{e}^2\mathbb{1}} % identity matrix: \usepackage{bbold_{2}\right).
\end{equation}

Using the weak shear approximation in \autoref{eq:weak_shear}, we can perform
the marginalization integral in \autoref{eq:marg_like} analytically,
\begin{align}
  \prf{{\mathbf{d}} | \Upsilon(\mathbf{x}), \alpha} &= \prod_{n=1}^{n_{\rm gal}}
    \int d^{2}\epsilon^\mathrm{int}_n\,
    \mathcal{N}_{\epsilon^\mathrm{int}}\left(\hat{e}_n - g(\mathbf{x}_n), \sigma_{{\rm ms};n}^2\mathbb{1}} % identity matrix: \usepackage{bbold_2\right)
    \notag\\
    &\quad\times
    \mathcal{N}_{\epsilon^\mathrm{int}}\left(0, \sigma_{e}^{2}\mathbb{1}} % identity matrix: \usepackage{bbold_{2}\right)
    \\
    &= \prod_{n=1}^{n_{\rm gal}}
    \mathcal{N}_{g(\mathbf{x}_n)}\left(\hat{e}_n, \left(\sigma_{{\rm ms};n}^2 + \sigma_{e}^2\right)\mathbb{1}} % identity matrix: \usepackage{bbold_{2}\right)
    \\
    &= \mathcal{N}_{\Matrix{g}}\left(\hat{\Matrix{e}}, \mathsf{N}\right),
    \label{eq:marg_like_analytic}
\end{align}
where in the final line we defined the $2n_{\rm gal}$-length vector of observed ellipticities $\hat{\Matrix{e}}$,
the same-length vector of reduced shear components at each galaxy locations $\Matrix{g}$, and
the diagonal $2n_{\rm gal}\times 2n_{\rm gal}$ dimensional covariance matrix $\mathsf{N}$ with the $n$th diagonal
entry equal to $\sigma^2_{{\rm ms};n} + \sigma^2_{e}$.
Note the weak shear likelihood does not depend on the lensing convergence $\kappa$ (because the
reduced shear is approximated as equal to the non-reduced shear and we ignore lensing magnification
effects on the galaxy sizes and fluxes).

\subsection{A maximum entropy prior for lensing fields}
\label{sec:framework}

We want an interim prior on the lensing convergence and shear that is not only
independent of cosmology but also broadly encompassing of the possible
cosmological interpretations and spatially varying systematics contributions.
That is, by choosing a functional form of a prior
for interpolating shear over the sky and marginalizing statistical uncertainties,
we should not restrict the class of physical models that might explain the data.
A mathematical version of this sentiment is that we want a
maximum entropy prior~\citep{PhysRev.106.620,shore1980} on the lensing fields.

For an assumed mean and (co-)variance, the Gaussian distribution is the maximum entropy
distribution~\citep{10.2307/2984828}. Therefore, because it is principally the second moment of the
lensing convergence that we use to constrain cosmological models, we choose a
Gaussian Process (GP) interim prior for sampling convergence (and shear) fields on
the sky given measurements of galaxy image moments.

Gravitational lensing shear is a spin-2 field and is non-locally
related to the lensing convergence, which both present modeling challenges.
However, both the lensing convergence and shear can be derived as second
derivatives of a scalar valued lensing potential $\psi$. So we impose the
GP prior on the potential $\psi$ and derive the related priors on the
convergence $\kappa$ and shear $\gamma$ from this starting assertion.

Each of the lensed observables, $\Upsilon(\mathbf{x}) \equiv \left[\kappa, \gamma_1, \gamma_2\right]$
is related to the lensing potential $\psi$ via derivatives in the form of:
\begin{align}
  \kappa &=\frac{1}{2}\left(\frac{\partial^2 \psi}{\partial x_1^2} +
    \frac{\partial^2 \psi}{\partial x_2^2 }\right) = \frac{1}{2}(\psi_{,11} +
    \psi_{,22})
  \label{eq:kappa_from_psi}
  \\
  \gamma_1 &=\frac{1}{2}\left(\frac{\partial^2 \psi}{\partial x_1^2} -
    \frac{\partial^2 \psi}{\partial x_2^2}\right) = \frac{1}{2}(\psi_{,11} -
    \psi_{,22})
  \label{eq:gamma1_from_psi}
  \\
  \gamma_2 &=\frac{1}{2}\left(\frac{\partial^2 \psi}{\partial x_1 \partial
    x_2} +\frac{\partial^2 \psi}{\partial x_2 \partial x_1}\right) =
    \frac{1}{2}(\psi_{,12} + \psi_{,21}).
  \label{eq:gamma2_from_psi}
\end{align}
It is straightforward to derive via integration by parts on the moments of the field
that if $\psi(\mathbf{x})$ is Gaussian
distributed for given $\mathbf{x}$, then so too is $\Upsilon(\mathbf{x})$, but with a
modified covariance.
By specifying a GP prior on the lens potential, we therefore can derive
a GP prior on the combination of lensing convergence and shear fields that preserves the
physical correlations between these fields. We will show that we can infer the lensing
convergence and shear via correlated draws from a GP distribution given a galaxy
ellipticity catalog. We do not then compute lensing convergence and shear from
spatial derivatives of a GP-distributed lens potential. The latter operation is only
precisely defined when we have (in principle) knowledge of the lens potential at all
sky locations. However, galaxies provide a non-contiguous background for measuring the lens
potential. The relations in \autoref{eq:kappa_from_psi}, \autoref{eq:gamma1_from_psi}, and
\autoref{eq:gamma2_from_psi} are then imposed only by the GP covariance structure during
sampling and marginalization.

The derivations of using a GP to represent the lensing convergence and shear fields were
first presented in \citet{ng2016}. We refer readers to \citet{ng2016} for an introduction to
the basics of a GP. We show the derivations from \citet{ng2016} in
\autoref{sec:Gaussian process covariances} for the convenience of the reader.
In particular, the GP prior for $\Upsilon$
derived from that for $\psi$ should not mix E and B modes in the two-point function
because we only allow for GP realizations that preserve the combinations of fields that satisfy
Equations~\ref{eq:kappa_from_psi}--\ref{eq:gamma2_from_psi} in the two-point correlations.
It is also straightforward
to extend the scalar valued potential $\psi$ to a complex-valued potential
$\psi\rightarrow\psi_{E} + i\psi_{B}$
to model or infer both E and B mode contributions to a measured shear signal.

We choose a squared exponential kernel for the GP model of the lensing potential $\psi$,
\begin{align}\label{eq:gp_cov}
  \mathsf{S}_{\psi}(\mathbf{x}, \mathbf{y}; \lambda_{\rm GP}, \ell^{2}_{\rm GP}) =

  \lambda_{\rm GP}^{-1} \exp\left(-\frac{1}{2} \frac{s^2(\mathbf{x},\mathbf{y})}{\ell^{2}_{\rm GP}}\right),
\end{align}
where $\lambda_{\rm GP}$ is a precision parameter that sets the amplitude of fluctuations of the GP,
$\ell^{2}_{\rm GP}$ is a squared distance defining the correlation length of the GP kernel,
and the squared distance $s^2$ between pairs of galaxy locations is,
\begin{equation}\label{eq:euclid_dist}
  s^2 \equiv (\mathbf{x} - \mathbf{y})^T (\mathbf{x} - \mathbf{y}).
\end{equation}
The squared exponential kernel is useful in interpolating smooth response functions between
the observed galaxy locations, as we expect for low resolution or low signal-to-noise ratio (SNR)
weak lensing mass reconstructions. Exploration of other kernel choices for different mass
reconstruction resolutions and SNRs is an interesting question that we leave for later work.
A practical motivation for our kernel choice is that the two parameters of the squared exponential
kernel can be optimized for different data sets without particular numerical or computational
challenges. A further justification for our smooth kernel choice comes from our intention to use the
GP as merely an interim prior for sampling lens fields. We will later describe how these interim lens
field realizations may be re-weighted under cosmologically informed priors.
We list the derivatives of \autoref{eq:gp_cov} in Appendix \ref{sec:Gaussian process covariances}
that are required to build the joint covariance for the lensing shear and convergence.
We will denote the covariance constructed from second derivatives of \autoref{eq:gp_cov} as
$\mathsf{S}_{\rm GP}$.

We now return to the pedagogical derivation of the likelihood function for galaxy ellipticity
estimators adding the GP distribution for the shear.
Because both the marginal likelihood in \autoref{eq:marg_like_analytic} and the shear prior are
Gaussian distributions in the shear, we can specify the shear (and convergence) joint posterior
for all galaxy and grid locations as a multivariate Gaussian distribution (otherwise known
as the Wiener filter),
\begin{multline}\label{eq:shear_posterior}
  \prf{\Upsilon(\mathbf{x},\mathbf{x}') | {\mathbf{d}}, \alpha, a_{\rm GP}}
  =
  \mathcal{N}_{\Upsilon(\mathbf{x},\mathbf{x}')}
  \left({\bm \mu}_{\Upsilon}, \mathsf{S}_{\Upsilon} \right)
  \\ \times
  \mathcal{N}_{\hat{\Matrix{e}}}\left(0, \mathsf{N} + \mathsf{S}_{\rm GP}\right),
\end{multline}
where,
\begin{align}
  \mathsf{S}_{\Upsilon} &\equiv \mathsf{S}_{\rm GP}\left(\mathsf{S}_{\rm GP} + \mathsf{N}\right)^{-1}\mathsf{N}
  \label{eq:posterior_cov}
  \\
  {\bm \mu}_{\Upsilon} &\equiv \mathsf{S}_{\rm GP} \left(\mathsf{S}_{\rm GP} + \mathsf{N}\right)^{-1} \hat{\Matrix{e}}.
  \label{eq:posterior_mean}
\end{align}
Calculation of ${\bm \mu}_{\Upsilon}$ in \autoref{eq:posterior_mean} yields a posterior mean
estimate of the lensing shear and convergence at every galaxy location. However, this estimator
requires inversion of an $2n_{\rm gal}\times2n_{\rm gal}$ covariance matrix, which can be computationally expensive.

Other works have attempted to reduce the dimensionality of the covariances
in \autoref{eq:shear_posterior} by interpolating and averaging
the measured galaxy ellipticities onto a grid of coarser resolution than that
sampled by the galaxy angular distribution~\citep[e.g.,][]{Alsing2015}.
However, such interpolations not only restrict the measured dynamic range,
but also can be expected to introduce artefacts in the inferred shear field on the grid based on
the shape of the smoothing kernel, to propagate
shape measurement systematics to a broad range of angular scales,
and to ignore error propagation from individual galaxy shape measurements to the shear inference.
Our method defines an explicit interpolation from galaxy to other sky locations, with error
propagation included.

We can marginalize the lensing fields at the galaxy locations to obtain the marginal
posterior for the lensing fields at just the smaller number of grid locations $\mathbf{x}'$,
\begin{multline}\label{eq:lens_posterior}
  \prf{\Upsilon(\mathbf{x}') | {\mathbf{d}}, \alpha, a_{\rm GP}} \propto
  \int d\Upsilon(\mathbf{x})\,
  \prf{{\mathbf{d}} | \Upsilon(\mathbf{x}), \alpha}
  \\ \times
  \prf{\Upsilon(\mathbf{x}') | \Upsilon(\mathbf{x}), a_{\rm GP}}
  \prf{\Upsilon(\mathbf{x}) | a_{\rm GP}}.
\end{multline}
When all distributions are Gaussians, \autoref{eq:lens_posterior} reduces to evaluating
\autoref{eq:shear_posterior} at only the locations of the grid $\mathbf{x}'$ given the input
ellipticities $\hat{\Matrix{e}}$. Note however that we still require in \autoref{eq:posterior_mean}
the evaluation of the shear model at every galaxy location, which is then interpolated by the
WF to arbitrary sky locations. Our approach is therefore different from algorithms that require
an initial averaging of galaxy ellipticity components over local sky regions.
\autoref{eq:lens_posterior} thus specifies the interim distribution of lensing shear and convergence
that can be later propagated to cosmological model analyses.

\subsubsection{Optimizing interpolation parameters}
\label{sub:optimizing_interpolation_parameters}

We can further marginalize the lens fields $\Upsilon$ at all locations $\mathbf{x}$ and $\mathbf{x}'$ to obtain
the marginal likelihood for the GP parameters,
\begin{equation}\label{eq:gp_marg_like}
  \prf{{\mathbf{d}} | \alpha, a_{\rm GP}} = \mathcal{N}_{{\mathbf{d}}}
  \left(0, \mathsf{S}_{\rm GP} + \mathsf{N}\right).
\end{equation}
We maximize the density in \autoref{eq:gp_marg_like} to determine suitable values of the
GP parameters for interpolation of the lens fields for a given ellipticity
catalog~\citep[see a similar approach in][]{Marshall:531310}.

\autoref{eq:gp_marg_like} is informative on the GP parameters if the data vector ${\mathbf{d}}$ includes
sufficient numbers of galaxies to beat down the shape noise. In cases with smaller
signal to noise ratios we can also add a cosmologically informed prior to help constrain
the GP parameters. In this case we can replace the data vector in \autoref{eq:gp_marg_like}
with a simulation of the data vector derived from a cosmological model. Then, by marginalizing
over realizations of the simulated data vector we get a marginal prior for the GP parameters,
\begin{align}\label{eq:marg_dist_cosmo}
  \prf{a_{\rm GP}| \theta, \alpha} &=
  \int d{\mathbf{d}}^{\rm sim}\,
  \prf{a_{\rm GP} | {\mathbf{d}}^{\rm sim}, \alpha} \prf{{\mathbf{d}}^{\rm sim} | \theta},
  \notag\\
  %
  &=
  \prf{a_{\rm GP}}
  \\
  &\times
  \int d{\mathbf{d}}^{\rm sim}\, \frac{\prf{{\mathbf{d}}^{\rm sim} | a_{\rm GP}, \alpha}}{\prf{{\mathbf{d}}^{\rm sim} | \alpha}}
  \prf{{\mathbf{d}}^{\rm sim} | \theta}.
  \notag
\end{align}
In \autoref{eq:marg_dist_cosmo} we specify a model for drawing simulated data realizations
${\mathbf{d}}^{\rm sim}$ given cosmological parameters $\theta$. A standard approach to such a model
is to first draw a realization of the 3D mass density perturbations from a Gaussian distribution
with a cosmologically determined initial power spectrum and then to evolve the initial
mass density perturbations under a numerical (e.g., $N$-body) simulation of gravitational
evolution.

Modeling only linear density perturbation evolution, the late-time model for the
projected lensing mass density, and hence the shear, is also a Gaussian distribution,
\begin{equation}\label{eq:linear_cosmo_model}
  \prf{{\mathbf{d}}^{\rm sim}|\theta}=\mathcal{N}_{{\mathbf{d}}^{\rm sim}}(0, \mathsf{\Sigma}(\theta)),
\end{equation}
for some cosmological covariance $\mathsf{\Sigma}(\theta)$.
To evaluate \autoref{eq:marg_dist_cosmo}  we also must determine the evidence
\begin{equation}\label{eq:evidence}
  \prf{{\mathbf{d}}^{\rm sim}|\alpha, I} \equiv
  \int da_{\rm GP}\, \prf{{\mathbf{d}}^{\rm sim} | \alpha, a_{\rm GP}}
  \prf{a_{\rm GP} | I},
\end{equation}
where $I$ denotes prior information on the GP parameters.
Even under a linear density perturbation approximation the evidence
cannot be calculated analytically.

For illustration of the GP parameter prior informed by cosmology, we use \autoref{eq:linear_cosmo_model}
to evaluate \autoref{eq:marg_dist_cosmo} via Monte Carlo integration. We draw $N$ samples
of ${\mathbf{d}}^{\rm sim}_{i}$ from $\prf{{\mathbf{d}}^{\rm sim}|\theta}$ and evaluate
$\prf{{\mathbf{d}}^{\rm sim}_{i}|a_{\rm GP}, \alpha}$ for each $i=1,\dots,N$. We approximate
the evidence in \autoref{eq:evidence} via 2D numerical integration and then calculate
\begin{equation}
  \prf{a_{\rm GP}| \theta, \alpha} \approx \prf{a_{\rm GP}}
  \frac{1}{N} \sum_{i=1}^{N} \frac{\prf{{\mathbf{d}}^{\rm sim}_{i}|a_{\rm GP}, \alpha}}
  {\prf{{\mathbf{d}}^{\rm sim}_{i}|\alpha, I}}.
\end{equation}
We show the resulting posterior constraints on the GP parameters in \autoref{fig:gp_cosmo_prior}.
\begin{figure}[!htb]
    \centerline{
        \includegraphics[width=0.4\textwidth]{gp_cosmo_prior.png}
    }
    \caption{Posterior constraints on the GP parameters after marginalizing realizations of
    simulated lens fields under a cosmological model. The correlations and variance in the simulated
    data serve to communicate the cosmological covariance structure into constraints on the
    GP parameters. The tight degeneracy shows the tradeoff between low correlations and large
    GP precision (equivalent to white noise as the GP precision becomes large) and larger correlations
    and smaller precision as the simulated maps are fit with smoother models. The lines show the 3-$\sigma$
    confidence interval.}
    \label{fig:gp_cosmo_prior}
\end{figure}
The cosmological simulation realizations impose a tight constraint in a linear combination
of the logarithm of the two GP parameters. The constraints in \autoref{fig:gp_cosmo_prior}
indicate the cosmological lens field simulations can be fit either with a large GP precision
and small correlation length (i.e., essentially as white noise), or with a smaller precision
and larger correlation lengths (i.e., as a smooth correlated map).
We thus obtain tight constraints on the GP parameters if we can include some prior knowledge
of the correlation length scale via $\prf{a_{\rm GP}}$.

See \autoref{sec:cosmology_dependent_covariance_model} for
a description of our linear theory cosmological model for the lensing convergence.

\subsection{Cosmological parameter inference}
\label{sub:cosmological_parameter_inference}

Under a standard cosmological model the lensing potential $\psi$ is related to the
3D cosmological late-time gravitational potential $\Psi^{\rm LT}$ by a projection along
the line of sight weighted by the lensing efficiency $K(a;A_s)$,
\begin{equation}\label{eq:lenspot_proj}
  \bar{\psi}_{s}(\mathbf{x}) \equiv

  W(\mathbf{x})
  \int dz\,
  \Psi^{\rm LT}(\mathbf{x},z)
  K(z ; Z_s),
\end{equation}
where $W(\mathbf{x})$ is the survey window function, $z$ is the cosmological redshift,
and $Z_s$ defines the redshift distribution of source galaxies for a galaxy sample $s$.
The late-time gravitational potential is in turn related to the potential of the
cosmological initial conditions $\Psi^{\rm IC}$ via a deterministic function $G$
defining gravitational evolution~\citep{2013MNRAS.432..894J},
\begin{equation}\label{eq:gravpot_ic}
  \Psi^{\rm LT}(\mathbf{x}, z) = G\left(\Psi^{\rm IC}, \theta, z\right).
\end{equation}
We assume the initial conditions $\Psi^{\rm IC}$ are Gaussian distributed with mean
zero and covariance $\mathsf{\Sigma}^{\rm IC}(\theta)$, which is entirely determined by the
potential power spectrum from inflation.

A common method to simulate $\Psi^{\rm LT}$ according to \autoref{eq:gravpot_ic} is
to run a cosmological $N$-body simulation with initial conditions $\Psi^{\rm IC}$
drawn from a Gaussian distribution with a specified power spectrum.
The numerical $N$-body solver allows evaluation of the function $G$, which does not
have a known analytic form beyond low-order perturbation theory.
The standard approach to evaluate \autoref{eq:lenspot_proj} for the lens potential
in a numerical simulation is to trace light rays through the observer's light cone
given the simulated $\Psi^{\rm LT}(\mathbf{x}, a)$. The bundle of light rays is
evaluated at a discrete set of sky locations to predict the lensing shear and
convergence. The predicted lens fields must then be interpolated over the sky to
galaxy locations to complete the numerical cosmological model prediction.

We have already defined a probabilistic interpolation of lens fields over the sky
via the GP. We therefore define the conditional distribution of the lens potential
at galaxy locations $\mathbf{x}$ given the simulations evaluated at positions $\mathbf{x}'$ via the
model of \autoref{eq:gravpot_ic} and \autoref{eq:lenspot_proj}
followed by GP interpolation,
\begin{equation}\label{eq:cond_model_comparison}
  \prf{\psi(\mathbf{x}) | \bar{\psi}(\mathbf{x}', \theta),a_{\rm GP}}
  = \mathcal{N}_{\psi} \left(
    \mu_{\psi}, \mathsf{\Sigma}_{\psi}
  \right),
\end{equation}
where,
\begin{align}
  \mu_{\psi} &\equiv \mathsf{S}(\mathbf{x},\mathbf{x}')\mathsf{S}^{-1}(\mathbf{x}',\mathbf{x}')
  \bar{\psi}_{s}(\mathbf{x}'; \Psi^{\rm IC}, \theta, A_s, W)
  \label{eq:mu_psi}
  \\
  \mathsf{\Sigma}_{\psi} &\equiv
  \mathsf{S}(\mathbf{x}, \mathbf{x}) - \mathsf{S}(\mathbf{x},\mathbf{x}')\mathsf{S}^{-1}(\mathbf{x}',\mathbf{x}')\mathsf{S}(\mathbf{x}',\mathbf{x}),
  \label{eq:sigma_psi}
\end{align}
define the mean and covariance for the conditional multivariate Gaussian.
With \autoref{eq:cond_model_comparison} we have thus derived a conditional probability
distribution to compare theoretical predictions of the lens potential with interim
samples of the potential drawn under the GP prior.
Said another way, our comparison of conditional posterior
samples of the lens fields with the cosmological models is mediated by the interpolation
over the sky using the GP.

By marginalizing the initial conditions realizations and the lens fields realizations
we can now derive the posterior distribution for the cosmological parameters,
\begin{multline}\label{eq:marg_post_cosmo}
  \prf{\theta | {\mathbf{d}}, a_{\rm GP}, \alpha} \propto
  \prf{\theta}
  \\ \times
  \int d\Upsilon\, \int d\Upsilon'\, \int d\Psi^{\rm IC}\,
  \delta_{D}\left(\Upsilon' - h(\Psi^{\rm IC})(\mathbf{x}')\right)
  \\ \times
  \frac{\prf{\Psi^{\rm IC}(\mathbf{x}') | \mathsf{\Sigma}^{\rm IC}(\theta)}}
    {\prf{\Upsilon' | a_{\rm GP}}}
  \\ \times
  \prf{{\mathbf{d}} | \Upsilon, \alpha}
  \prf{\Upsilon | \Upsilon', a_{\rm GP}} \prf{\Upsilon' | a_{\rm GP}},
\end{multline}
where $h(\Psi^{\rm IC})$ indicates the deterministic gravitational evolution of the
initial conditions potential to late times where the lensing is observed.

The final line of \autoref{eq:marg_post_cosmo} is the interim sampling distribution
we defined in the previous section. We can thus perform the integrals in
\autoref{eq:marg_post_cosmo} via Monte Carlo with the same interim lens field
samples that are generated in the map making algorithm.
\begin{figure}[!htb]
  \centerline{
    \includegraphics[width=0.5\textwidth]{shear_gp_pgm_cosmo_split.png}
  }
  \caption{

    In our approximate cosmological parameter inference pipeline we assume
    the lens fields at the galaxy locations $\Upsilon$ are entirely determined
    by GP interpolation from the lens fields at a regular set of sky locations
    $\Upsilon'$ at which we compute theory predictions.
    The lens fields on the regular set are informed by the cosmological
    model with parameters $\theta$.

  }
  \label{fig:pgm_cosmo}
\end{figure}

We propose two sampling algorithms depending on whether or not we can directly evaluate the
probability density for the lens fields under the cosmological model.
\begin{description}
  \item[Sampling algorithm 1] This is an algorithm to use when it is possible to analytically
  evaluate the cosmological probability density function.
  With $K$ samples of $\Upsilon(\mathbf{x}')$ from the interim sampling distribution,
  we use importance sampling to approximate the marginalizations in
  \autoref{eq:marg_post_cosmo},
  \begin{multline}\label{eq:marglikecosmo}
    \prf{\theta | {\mathbf{d}}, a_{\rm GP}, \alpha} \approx
    \prf{\theta}

    \frac{1}{K}\sum_{k=1}^{K}
    \frac{\prf{\Upsilon_{k}(\mathbf{x}') | \theta}}
      {\prf{\Upsilon_{k}(\mathbf{x}') | a_{\rm GP}}},
  \end{multline}
  where from \autoref{eq:marg_post_cosmo},
  \begin{multline}
    \prf{\Upsilon(\mathbf{x}') | \theta} \equiv
    \int d\Psi^{\rm IC}\,
    \prf{\Psi^{\rm IC}(\mathbf{x}') | \mathsf{\Sigma}^{\rm IC}(\theta)}
    \\ \times
    \delta_{D}\left(\Upsilon' - h(\Psi^{\rm IC}(\mathbf{x}'))\right).
  \end{multline}

  \item[Sampling algorithm 2] This is an algorithm to use when cosmological modeling of the lens
  fields is only possible via forward simulation.
  For most cosmological models of interest, we can only predict the model for the
  lens fields via forward simulation of the cosmological mass density perturbation evolution.
  In this case, we have no direct mechanism to evaluate the density
  $\prf{\Upsilon(\mathbf{x}')|\theta}$.
  Instead, we follow the steps,
  \begin{enumerate}
    \item Draw $\Psi^{\rm IC}(\mathbf{x}')$ from $\prf{\Psi^{\rm IC}(\mathbf{x}') | \mathsf{\Sigma}^{\rm IC}(\theta)}$.
    \item Compute the predicted lens fields given the drawn initial conditions
    $h(\Psi^{\rm IC})$ (e.g., via $N$-body simulation and ray-tracing prediction of the lensing
    shear and convergence).
    \item Select $K$ samples of $\Upsilon(\mathbf{x})$ from the set of interim samples under the
    GP prior.
    \item Evaluate the density, $\prf{\Upsilon|\Upsilon', a_{\rm GP}}$, as in
    \autoref{eq:cond_model_comparison} for each sample $K$ and compute,
    \begin{equation}
      \prf{\theta | {\mathbf{d}}, a_{\rm GP}, \alpha} \approx \prf{\theta}
      \frac{1}{K} \sum_{k=1}^{K}
      \prf{\Upsilon_{k}(\mathbf{x})|\Upsilon(\mathbf{x}',\theta), a_{\rm GP}}.
    \end{equation}
    Note that $\Upsilon'$ here is that from the cosmological forward model, not that from the
    interim sampling.
  \end{enumerate}
\end{description}

For Sampling Algorithm 1, if we assume linear cosmological perturbations we can approximate,
  \begin{equation}
    \prf{\Upsilon(\mathbf{x}')| \theta} \approx
    \mathcal{N}_{\Upsilon(\mathbf{x}')}
    \left(\mathbf{0}, \mathsf{\Sigma}(P_{\Upsilon}(\theta))\right).
  \end{equation}
  Under this approximation we can combine the terms from \autoref{eq:marg_post_cosmo},
  \begin{equation}\label{eq:lens_field_cond_factor}
    \prf{\Upsilon | \Upsilon',a_{\rm GP}} \prf{\Upsilon'|\theta} =
    \mathcal{N}_{\Upsilon(\mathbf{x},\mathbf{x}')} \left(
      0, \mathsf{M}(\mathbf{x},\mathbf{x}'; a_{\rm GP}, \theta)
    \right),
  \end{equation}
  where $\mathsf{M}(\mathbf{x},\mathbf{x}'; a_{\rm GP}, \theta)$ is defined in \autoref{eq:signal_cov_cosmo}.

  If we further assume, as above, that the likelihood function is a Gaussian distribution in the
  galaxy ellipticities with a linear weak shear applied, then \autoref{eq:marg_post_cosmo} becomes,
  \begin{equation}\label{eq:marg_post_cosmo_linear}
    \prf{\theta | {\mathbf{d}}, a_{\rm GP}, \alpha} \propto \prf{\theta}
    \mathcal{N}_{{\mathbf{d}}}\left(
      0, \mathsf{M}(a_{\rm GP}, \theta) + \mathsf{N}
    \right),
  \end{equation}
  with the signal covariance as defined in \autoref{eq:signal_cov_cosmo}.

  To summarize, the approximations required for \autoref{eq:marg_post_cosmo_linear} are,
  \begin{enumerate}
    \item a likelihood function that is Gaussian in the galaxy ellipticities,
    \item galaxy ellipticity measurements that are unbiased estimators of the reduced shear,
      \item weak shear ($\kappa \ll 1$),
    \item linear cosmological perturbations.
  \end{enumerate}
  We can drop assumptions 1, 2, and 3 and still use Sampling Algorithm 1. But to avoid assumption
  4 we must resort to Sampling Algorithm 2.

\section{Results}
\label{sec:results}

We evaluate the mean and covariance of the marginal posterior for the lensing shear and convergence
from \autoref{eq:shear_posterior} for a simulated data set, where we can compare with a known
truth (\autoref{sub:simulation_study}). We also analyze
the ellipticity catalog~\citep{jee2013} in the
Deep Lens Survey (DLS)\footnote{\url{http://dls.physics.ucdavis.edu}}~\citep{2002SPIE.4836...73W}
in the vicinity of a massive galaxy cluster where the
amplitudes of the lensing fields are large (\autoref{sub:abell}).

\subsection{Simulation study}
\label{sub:simulation_study}

To validate that \autoref{eq:posterior_mean} can recover the correct shear and convergence fields,
we create simulated galaxy ellipticity catalogs with artificially small shape noise so we can
measure the lensing shear and convergence to a tuneable precision.

Our procedure for simulating galaxy ellipticity catalogs is,
\begin{enumerate}
  \item Calculate a lensing convergence angular power spectrum using the cosmology theory
  code {\sc CHOMP}\xspace\footnote{\url{https://github.com/karenyyng/chomp}}~\citep{2013JCAP...11..009M}.
    We assume a standard
  $\Lambda$CDM cosmology with $\Omega_m=0.3$, $\sigma_8=0.8$ and a source redshift distribution
  with a narrow peak at $z=1$.
  \item Simulate Gaussian-distributed lensing shear and convergence maps on a grid using the
  code {\sc GalSim}\xspace\footnote{\url{https://github.com/GalSim-developers/GalSim}}~\citep{galsim}.
  \item Place one galaxy in each grid cell of the simulated lensing shear maps. These galaxies are simply
  sources of illumination for measuring the lensing fields, not cosmologically clustered galaxies.
  \item For each galaxy, draw intrinsic ellipticity components from a 2D Gaussian distribution
  with mean 0 and a specified variance, $\sigma_e^2$.
  \item Calculate lensed ellipticities for each galaxy by adding the lensing shear to the intrinsic
  ellipticities, assuming a weak shear approximation.
  \item Save the galaxy angular sky locations and ellipticity components to a catalog file.
\end{enumerate}
Given a simulated ellipticity catalog, we find the GP parameters
$\lambda_{\rm GP}, \ell^{2}_{\rm GP}$ that maximize \autoref{eq:gp_marg_like}.
We then evaluate \autoref{eq:posterior_mean} and \autoref{eq:posterior_cov} using the optimized
GP parameters to obtain the
marginal posterior distribution of the lensing fields $\kappa, \gamma$ at all galaxy locations
as well as on a regular grid of locations.
We expect that in practical applications to large data sets the shear will only need to be
evaluated at the galaxy locations and the converngence (or lens potential) will only need to be
evaluated at a smaller number of sky grid locations, thus reducing the overall dimensionality of
the linear system to be solved.

We show an example of the output of this procedure in
\autoref{fig:mass_map_comparison_galsim}
compared to the input shear fields used to generate the mock ellipticity catalog. In the
left column in \autoref{fig:mass_map_comparison_galsim} we show the posterior mean fields from
\autoref{eq:posterior_mean}.
In the adjacent panel we show the input shear fields.
We show the `true' convergence that we calculated at the same time as the input shear, but we
do not use the convergence at any point in our calculation. The ``estimated convergence'' comes
from interpolating the measured galaxy ellipticities with the GP kernel.
The right panel in \autoref{fig:mass_map_comparison_galsim} shows the
signal-to-noise ratio (SNR) of the lens field maps, defined as,
\begin{equation}\label{eq:snr_def}
  {\rm SNR} \equiv \Upsilon / \sqrt{{\rm diag}\left(\mathsf{S}_{\Upsilon}\right)},
\end{equation}
with $\mathsf{S}_{\Upsilon}$ defined in \autoref{eq:posterior_cov}.
In the limit that the variance of the lens fields dominates the intrinsic shape
and ellipticity measurement variances, \autoref{eq:posterior_cov} reduces to
$\mathsf{N}$. Then \autoref{eq:snr_def} becomes
${\rm SNR} \rightarrow \Upsilon / \sqrt{\sigma_e^2 + \sigma_{\rm ms}^{2}}$,
which is similar to SNR definitions in other weak lensing mass mapping
analyses~\citep{2000MNRAS.313..524V,2014MNRAS.442.2534S}.

The simulation in
\autoref{fig:mass_map_comparison_galsim} has an artificially low intrinsic ellipticity r.m.s.\
of $\sigma_e=0.0026$ compared to $\sigma_e=0.26$ in the Deep Lens Survey that
we analyze in \autoref{sub:abell}.
We choose a small shape noise r.m.s.\ to
allow us to validate our convergence inference with a small number of only 1600 galaxies.
We place the simulated galaxies on a $40\times40$ grid, so that we do not have to interpolate
the simulated shear fields to build the mock ellipticity catalog.

Because the Gaussian model for the shape noise r.m.s.\ scales with the number of galaxies
$n_{\rm gal}$ as $n_{\rm gal}^{-1/2}$, the shape noise r.m.s.\ in
our simulation is equivalent to a DLS-like galaxy sample with 16 million galaxies.
This is about 64 times the number of galaxies we have in a single four square degree field
of the DLS. We show the effect of increasing the shape noise r.m.s.\ by a factor
of $\sqrt{64}$ later in this section.
\begin{figure*}[!htb]
  \centerline{
    \includegraphics[width=0.54\textwidth]{mass_map_comparison_galsim.png}
    \includegraphics[width=0.35\textwidth]{mass_map_snr_galsim.png}
  }
  \caption{Left: Comparison of the convergence and shear maps in our simulation study between that
  used to generate our mock galaxy ellipticity catalog (right column) and the output of our
  GP interpolation (left column). The rows show the maps for the two shear components $\gamma_{1,2}$
  and the convergence $\kappa$. These maps cover a $2\times2$ square degree field. The simulated
  intrinsic ellipticity r.m.s.\ is set to an artificially small value of $\sigma_e=0.0026$,
  which is 100 times smaller than that observed for the complete Deep Lens Survey catalog.
  Right: signal-to-noise ratio (SNR) maps for the same simulations. We calculate SNR as the
  ratio of the map to the square root of the diagonal of the covariance
  in \autoref{eq:posterior_cov}.
  These simulations use 1600 galaxies and grid sizes for the GP interpolated shear and convergence
  with 24 grid cells per dimension.
  }
  \label{fig:mass_map_comparison_galsim}
\end{figure*}

The reconstructed shear and convergence maps in \autoref{fig:mass_map_comparison_galsim} are a
good match to the simulation inputs, but are noticeably smoothed. This is for two reasons;
(1) we evaluate the interpolated lens fields on a smaller $24\times24$ grid than the
$40\times40$ grid on which the inputs are evaluated,
(2) we use a value of $\ell^{2}_{\rm GP}=0.0123$ that is larger than the Nyquist scale in the maps, even
for the $24\times24$ interpolated grid because of the way the GP parameters are optimized.
We will discuss the GP parameter optimization below.

In \autoref{fig:power_spectra_galsim} we show the E and B mode power spectrum estimators obtained
from the posterior mean shear fields shown in \autoref{fig:mass_map_comparison_galsim}.
We compare the E-mode power spectrum estimator from the mean posterior shear maps with that
using the higher-resolution shear maps that were used to generate the mock ellipticity
catalog (which we label `simulation truth'). We also show in \autoref{fig:power_spectra_galsim}
the `theory' power spectrum that we used to generate the `simulation truth' shear maps.
The `theory' and `simulation truth' spectra agree on scales below the Nyquist frequency,
shown in \autoref{fig:power_spectra_galsim} as the vertical blue line. The mean posterior
power spectrum agrees with the `simulation truth' spectrum on scales below the effective
smoothing frequency derived from the value of $\ell^{2}_{\rm GP}$ and shown by the vertical
dot-dashed line in \autoref{fig:power_spectra_galsim}.
\begin{figure}[!htb]
  \centerline{
    \includegraphics[width=0.5\textwidth]{power_galsim_lownoise.png}
  }
  \caption{Comparison of power spectrum estimators in our simulation study to the `truth' input
  spectrum.
  The vertical solid lines show the Nyquist frequencies of the grids.
  The vertical dot-dashed line shows the multipole corresponding to the length scale
  $\sqrt{\ell_{\rm GP}^{2}}$ set by the GP kernel parameter $\ell_{\rm GP}^{2}=0.0123$ for
  field coordinates are normalized to the unit square.
  }
  \label{fig:power_spectra_galsim}
\end{figure}

The B-mode power spectrum estimated
from the mean posterior shear maps is shown by the dashed red line in \autoref{fig:power_spectra_galsim}.
We expect the B-mode power to be consistent with zero because we simulated only E-mode power.
The nonzero B-mode power spectrum in \autoref{fig:power_spectra_galsim} is explained by
examining \autoref{fig:EB_maps}, which shows the E and B mode mean posterior maps from which the
E and B mode mean posterior power spectra were derived. The B-mode map in \autoref{fig:EB_maps}
is near zero throughout the field except near the boundaries. These edge effects in the B-mode map
can be explained by mathematical ambiguities in the definition of E and B mode separation in a
finite field~\citep{bunn2003}. The small value of the B-mode map in \autoref{fig:EB_maps} away from
the field boundaries indicates our GP interpolation method does not create spurious B-modes.
\begin{figure*}
  \centerline{
    \includegraphics[width=0.6\textwidth]{EB_maps_galsim.png}
  }
  \caption{E and B mode maps derived from the posterior mean shear maps shown in the left column
  of \autoref{fig:mass_map_comparison_galsim}.
  The E-mode map closely matches the interpolated convergence map in
  \autoref{fig:mass_map_comparison_galsim} as expected.
  The B-mode map is near zero throughout the center of the field, but shows non-zero
  values around the field edge because of mathematical ambiguities in the E/B mode separation
  at the field boundaries~\citep{bunn2003}. Note the value of the plotted B-mode map
  has been multiplied by 10 for easier visualization.
  }
  \label{fig:EB_maps}
\end{figure*}

To illustrate the GP parameter optimization procedure we show in \autoref{fig:sigma_e_maps} the
mean posterior convergence, the associated SNR maps, and the marginal likelihood surface for
the GP parameters (left to right columns) for increasing intrinsic ellipticity variance
$\sigma_e^2$ (top to bottom rows). The top row of panels in \autoref{fig:sigma_e_maps} show a repeat
of the $\sigma_e^2$ and $a_{\rm GP}$ values used in \autoref{fig:mass_map_comparison_galsim}.
For the very small shape noise in this case, there is a narrow peak in the log-likelihood surface
for the two GP kernel parameters. We thus select the maximum likelihood (ML) values for $a_{\rm GP}$
and obtain the convergence and SNR maps that closely resemble the `true' convergence
as shown in \autoref{fig:mass_map_comparison_galsim}. However, the ML value of $\ell^{2}_{\rm GP}$ is somewhat
larger than the Nyquist frequency of the grid to which we interpolate as shown in
\autoref{fig:power_spectra_galsim} (compare dot-dashed to red solid vertical lines). Because the
noise is sub-dominant in this example, we would expect a value of $\ell^{2}_{\rm GP}$ matching the grid
Nyquist frequency to yield a more accurate convergence map reconstruction. We are likely to
obtain more accurate results, therefore, if we impose a prior on $a_{\rm GP}$ that encodes this
expectation.
\begin{figure*}
  \centerline{
    \subfigimg[width=0.3\textwidth,pos='LOWER LEFT',hsep=15pt,font=\footnotesize]{\textcolor{white}{$\sigma_{e}=0.0026$}}{kappa_galsim_sigmae1.png}
    \includegraphics[width=0.3\textwidth]{snr_galsim_sigmae1.png}
    \subfigimg[width=0.3\textwidth,pos='LOWER LEFT',hsep=15pt,font=\footnotesize]{\textcolor{white}{$\lambda_{\rm GP}=2.2\times 10^{8}, \ell^2=0.012$}}{lnp_galsim_sigmae1.png}
  }
  \centerline{
    \subfigimg[width=0.3\textwidth,pos='LOWER LEFT',hsep=15pt,font=\footnotesize]{\textcolor{white}{$\sigma_{e}=0.0077$}}{kappa_galsim_sigmae2.png}
    \includegraphics[width=0.3\textwidth]{snr_galsim_sigmae2.png}
    \subfigimg[width=0.3\textwidth,pos='LOWER LEFT',hsep=15pt,font=\footnotesize]{\textcolor{white}{$\lambda_{\rm GP}=2.3\times 10^{7}, \ell^2=0.061$}}{lnp_galsim_sigmae2.png}
  }
  \centerline{
    \subfigimg[width=0.3\textwidth,pos='LOWER LEFT',hsep=15pt,font=\footnotesize]{\textcolor{white}{$\sigma_{e}=0.013$}}{kappa_galsim_sigmae3.png}
    \includegraphics[width=0.3\textwidth]{snr_galsim_sigmae3.png}
    \subfigimg[width=0.3\textwidth,pos='LOWER LEFT',hsep=15pt,font=\footnotesize]{\textcolor{white}{$\lambda_{\rm GP}=2.3\times 10^{7}, \ell^2=0.061$}}{lnp_galsim_sigmae3.png}
  }
  \centerline{
    \subfigimg[width=0.3\textwidth,pos='LOWER LEFT',hsep=15pt,font=\footnotesize]{\textcolor{white}{$\sigma_{e}=0.023$}}{kappa_galsim_sigmae5.png}
    \includegraphics[width=0.3\textwidth]{snr_galsim_sigmae5.png}
    \subfigimg[width=0.3\textwidth,pos='LOWER LEFT',hsep=15pt,font=\footnotesize]{\textcolor{white}{$\lambda_{\rm GP}=2.3\times 10^{6}, \ell^2=0.30$}}{lnp_galsim_sigmae5.png}
  }
  \caption{Posterior mean convergence maps (left), signal-to-noise ratio (SNR) maps (middle),
  and marginal likelihood contours (right) for simulated galaxy catalogs with varying intrinsic
  ellipticity r.m.s.\ $\sigma_e$.
  From top to bottom, $\sigma_e = \left(0.00258, 0.00774, 0.0129, 0.0232\right)$ (these values are
  also annotated in the bottom left corner of each convergence map).
  The top row matches the simulation shown in \autoref{fig:mass_map_comparison_galsim}.
  Following the right column from top to bottom shows how the maximum-likelihood estimates for the
  GP parameters shifts to longer correlation lengths and smaller precision parameters with increasing
  shape noise. That is, as the data becomes more noise dominated, the marginal GP parameter likelihood
  changes shape to prefer smoothing, or effectively averaging, more galaxies to retain a more
  significant shear and convergence signal. The ML values for the GP parameters are listed
  in each panel showing the log-likelihood contours.
  }
  \label{fig:sigma_e_maps}
\end{figure*}

As the shape noise increases (for the same input signal and measurement uncertainties), the peak
in the marginal log-likelihood surface for $a_{\rm GP}$ becomes broader and eventually disappears
as shown by the rows of panels from top to bottom in \autoref{fig:sigma_e_maps}. The ML value
for $\ell^{2}_{\rm GP}$, while less well defined, continues to yield maps that are more smoothed as
$\sigma_e^2$ increases. This helps to preserve large amplitudes of the peaks in the SNR maps (see
the color bar scales in the SNR maps of \autoref{fig:sigma_e_maps}), but
in compensation erases structures at all but the lowest spatial frequencies in the maps. This
procedure, with flat priors in the log of $a_{\rm GP}$, appears useful for visualizing the posterior
convergence maps, but is undesirable for subsequent cosmological analyses. We see again that
we would prefer a prior favoring smaller $\ell^{2}_{\rm GP}$ even as the shape noise becomes large.

We assert such a prior in \autoref{fig:sigma_e_maps_with_prior}, by imposing Gaussian priors
separately in $\ln(\lambda_{\rm GP})$ and $\ln(\ell^{2}_{\rm GP})$ with parameters given in
\autoref{tab:prior_params_galsim}. Our Gaussian prior is informed by the cosmological simulation study
shown in \autoref{fig:gp_cosmo_prior} combined with our prior that the GP correlation length be large enough
so that white noise does not dominate the fits to the lens fields.
For each value of $\sigma_e$ in \autoref{fig:sigma_e_maps_with_prior}
we see that the signal-to-noise ratio (SNR) is comparable to that in \autoref{fig:sigma_e_maps}
but the convergence maps include higher spatial frequency structures.
\begin{table}[!htb]
  \begin{center}
    \caption{\label{tab:prior_params_galsim}Parameters for the Gaussian prior on the
      GP parameters for our simulation study.}
    \begin{tabular}{lcc}
      \hline\hline
      Parameter & mean & std. dev. \\
      \hline
      $\ln\left(\lambda_{\rm GP}\right)$ & $18$ & 0.43 \\
      $\ln\left(\ell^{2}_{\rm GP}\right)$    & $-4.0$ & 0.1\\
      \hline\hline
    \end{tabular}
  \end{center}
\end{table}

\begin{figure*}
  \centerline{
    \subfigimg[width=0.3\textwidth,pos='LOWER LEFT',hsep=15pt,font=\footnotesize]{\textcolor{white}{$\sigma_{e}=0.0026$}}{kappa_galsim_prior_sigmae1.png}
    \includegraphics[width=0.3\textwidth]{snr_galsim_prior_sigmae1.png}
    \subfigimg[width=0.3\textwidth,pos='LOWER LEFT',hsep=15pt,font=\footnotesize]{\textcolor{white}{$\lambda_{\rm GP}=2.2\times 10^{8}, \ell^2=0.012$}}{lnp_galsim_prior_sigmae1.png}
  }
  \centerline{
    \subfigimg[width=0.3\textwidth,pos='LOWER LEFT',hsep=15pt,font=\footnotesize]{\textcolor{white}{$\sigma_{e}=0.0077$}}{kappa_galsim_prior_sigmae2.png}
    \includegraphics[width=0.3\textwidth]{snr_galsim_prior_sigmae2.png}
    \subfigimg[width=0.3\textwidth,pos='LOWER LEFT',hsep=15pt,font=\footnotesize]{\textcolor{white}{$\lambda_{\rm GP}=2.3\times 10^{7}, \ell^2=0.027$}}{lnp_galsim_prior_sigmae2.png}
  }
  \centerline{
    \subfigimg[width=0.3\textwidth,pos='LOWER LEFT',hsep=15pt,font=\footnotesize]{\textcolor{white}{$\sigma_{e}=0.013$}}{kappa_galsim_prior_sigmae3.png}
    \includegraphics[width=0.3\textwidth]{snr_galsim_prior_sigmae3.png}
    \subfigimg[width=0.3\textwidth,pos='LOWER LEFT',hsep=15pt,font=\footnotesize]{\textcolor{white}{$\lambda_{\rm GP}=2.2\times 10^{8}, \ell^2=0.027$}}{lnp_galsim_prior_sigmae3.png}
  }
  \centerline{
    \subfigimg[width=0.3\textwidth,pos='LOWER LEFT',hsep=15pt,font=\footnotesize]{\textcolor{white}{$\sigma_{e}=0.023$}}{kappa_galsim_prior_sigmae5.png}
    \includegraphics[width=0.3\textwidth]{snr_galsim_prior_sigmae5.png}
    \subfigimg[width=0.3\textwidth,pos='LOWER LEFT',hsep=15pt,font=\footnotesize]{\textcolor{white}{$\lambda_{\rm GP}=2.2\times 10^{8}, \ell^2=0.027$}}{lnp_galsim_prior_sigmae5.png}
  }
  \caption{Same as \autoref{fig:sigma_e_maps} except with a prior asserted for the GP
  parameters $a_{\rm GP}$ that limits the degree of smoothing in the posterior mean maps.
  }
  \label{fig:sigma_e_maps_with_prior}
\end{figure*}

\subsubsection{Cosmological parameter constraints}
\label{sub:cosmo_params}

Although we assert a non-cosmological GP prior on the lens fields to infer regular convergence
and shear maps, we now demonstrate how we can recover cosmological information in a manner
similar to common algorithms in the literature. That is, we compute an angular power spectrum
estimator for the lensing convergence from the lens field posterior distribution.

The posterior distribution for the lens fields given the GP parameters is a multivariate
Gaussian distribution characterized by a mean field and covariance as given in
\autoref{eq:posterior_mean} and \autoref{eq:posterior_cov}. The posterior mean is thus a
convenient and useful summary statistic (as we have shown above).
Also, because an isotropic Gaussian random field is fully described by the angular power spectrum,
it is common in the literature to reduce cosmological large-scale structure statistics
to two-point function estimators for cosmological parameter estimation.

We showed the angle-averaged (E-mode) convergence power spectrum estimator for our low-noise
Gaussian simulated maps in \autoref{fig:power_spectra_galsim}. We use this power spectrum
estimator as a summary statistic derived from the observed ellipticity catalog. We further
assert a multivariate Gaussian likelihood function for the power spectrum estimator
with covariance~\citep[e.g.,][]{2001ApJ...554...56C},
\begin{equation}
  \mathrm{Cov}\left(P_{\kappa}(\ell)\right) = \frac{2}{N_{\ell}}P_{\kappa}^{2}(\ell),
\end{equation}
where $N_{\ell}$ is the number of modes contributing to the band power estimator for a multipole
bin centered at $\ell$.

In \autoref{fig:cosmo_param_contours} we show 68\% and 95\% contours of the 2D posterior
distribution for the cosmological mass density $\Omega_m$ and density fluctuations r.m.s.\
$\sigma_8$ given the angular power spectrum likelihood just described and flat priors.
The values used to generate the mock data are $\Omega_m=0.3, \sigma_8=0.8$.
We limit the multipole range of the power spectrum in the likelihood to $100 < \ell < 1300$,
where the upper bound is set by the effective smoothing length imposed by the asserted GP correlation
parameter $\ell^{2}_{\rm GP}$.
\begin{figure}
  \centerline{
    \includegraphics[width=0.47\textwidth]{cosmo_param_contours.png}
  }
  \caption{Marginal constraints on the cosmological parameters using the E-mode power spectrum
  estimator derived from the mean posterior convergence field.}
  \label{fig:cosmo_param_contours}
\end{figure}

We see from \autoref{fig:cosmo_param_contours} that the GP interim prior used to derive the
convergence map from the galaxy ellipticity catalog has not biased the cosmological
parameter constraints (within the uncertainties) obtained from a reduced summary statistic
of the lens field posterior. We defer to later work more complete demonstrations of the
cosmological parameter inference algorithms in \autoref{sub:cosmological_parameter_inference}
based on marginalizing the lens field realizations.

\subsection{Abell cluster in the Deep Lens Survey}
\label{sub:abell}

We apply our lens field inference algorithm to the galaxy ellipticity catalog derived from the
DLS\footnote{\url{http://dls.physics.ucdavis.edu}}~\citep{2002SPIE.4836...73W}.
The DLS is a 20 square degree optical imaging survey optimized for cosmic shear measurements.
We analyze a $\sim1$~square~degree field centered on the Abell cluster 781, which was
previously analyzed using DLS lensing shear measurements in
\citet{2006ApJ...643..128W,2009ApJ...702..980K,wittman2014}.
Abell 781 consists of four massive galaxy clusters, which is a useful case study for our
algorithm because it provides a large lensing signal with a modest number of galaxies
and because the distribution of mass density perturbations is decidedly not Gaussian
distributed.

\citet{jee2013} presented galaxy ellipticity measurements with the DLS $R$-band imaging calibrated
to produce shear estimates with biases well below the statistical uncertainties for two-point
cosmic shear correlation function estimators. The DLS shear pipeline in \citet{jee2013} includes
correlated PSF size and ellipticity corrections in each DLS exposure, calibration of additive
and multiplicative shear biases with image simulations, and a set of null tests validating the PSF
shear calibration corrections. The shear estimation pipeline used in \citet{jee2013}, \textsf{sFIT},
was further validated as the winning algorithm in the blinded community shear measurement
challenge \textsf{GREAT3}~\citep{great3-paper1}.
The DLS was performed with four optical pass-bands ($BVRz$) that allow photometric
redshift (photo-$z$) estimates for all lensing source galaxies~\citep{schmidt2013}.
\citet{jee2015} extended the DLS shear analysis to a tomographic cosmic shear measurement using
the photo-$z$ estimates.

The DLS galaxy ellipticity catalog produced for \citet{jee2013} includes a catalog-level selection
based on measured galaxy magnitudes, sizes, ellipticity measurement error, photo-$z$ estimates, and
proximity to masks as listed in \citet{jee2013} Table 2.
We perform a further set of selections on this catalog as listed in
\autoref{tab:dls_selection}.
\begin{table}[!htb]
  \begin{center}
    \caption{\label{tab:dls_selection}Selection criteria applied to the DLS galaxy catalog.}
    \begin{tabular}{ccc}
    \hline\hline
    Parameter        & min         & max \\
    \hline
    $z_b$            & 0.45        & --  \\
    $de$             & --          & 0.1 \\
    $\sqrt{a^2+b^2}$ & $0.8^{''}$  & --  \\
    $R$              & 22          & 23  \\
    \hline\hline
    \end{tabular}
  \end{center}
\end{table}
We choose the lower bound on (maximum posterior) photo-$z$, $z_b$, to select source galaxies
that are likely to be at redshifts larger than that of the highest redshift sub-cluster in the field
at $z\approx0.43$~\citep{wittman2014}. We further select only those galaxies in the ellipticity
catalog with ellipticity measurement errors, $de$, less than 0.1 and sizes greater than 0.8~arcseconds to
obtain galaxies likely to have more precisely measured shapes for informing the lensing shear.
We exclude galaxies with sizes, as determined from the geometric mean of the semi-major and semi-minor axes $a,b$,
less than 0.8 arcseconds because the ellipticities tend to be less
well measured when the galaxy size is similar to that of the PSF.
We select the brighter galaxies based on $R$-band magnitude that are still likely to be faint
enough to avoid significant contamination from cluster members. After all the selections listed
in \autoref{tab:dls_selection} we measure an ellipticity r.m.s.\ of $\sigma_e=0.21$. This measurement
includes the ellipticities with lensing effects included, but because lensing is sub-dominant
to the intrinsic ellipticity dispersion we assert $\sqrt{\alpha} = \sigma_e=0.21$ for the posterior
inference for A781.

We show our mean posterior inference of the lensing convergence of A781 in
\autoref{fig:mass_map_Abell} using 6000 galaxies randomly selected from the cut sample described
in \autoref{tab:dls_selection}. We select only 6000 galaxies to limit the size of the lens field
joint covariance that we must invert.
The left column of panels shows the convergence while the right
column shows the signal-to-noise ratio for the mean posterior. Our calculation is in the
top row of panels in \autoref{fig:mass_map_Abell}, which we compare with the aperture mass
algorithm of \citet{2012ApJ...747L..42D} in the bottom row of panels.
Note we use the same galaxy sample as input to each mass mapping algorithm.
The algorithm of \citet{2012ApJ...747L..42D},
called `aperture densomitry', provides a mass estimator that is more localized on the sky. The shear
in apertures is averaged with weighting functions that account for both angular selections and
the expected line-of-sight lensing kernel with the aide of the photometric redshift information
in the DLS catalog. However, to make a more direct comparison with the algorithm in this paper,
we recomputed the aperture densomitry weights without using any photometric redshift information
for the source galaxies.

\autoref{fig:mass_map_Abell} shows that we obtain consistent results for the two main A781
sub-clusters using our mean posterior map and the method used in a previous analysis. The white
crosses in \autoref{fig:mass_map_Abell} indicate the locations of all sub-clusters detected
in \citet{wittman2014}. We do not detect all the same sub-clusters, which is likely
because we use a significantly
smaller number of galaxies (6000 versus $\sim 50000$) while we test the performance and scaling
of our codes. However one sub-cluster denotedin \autoref{fig:mass_map_Abell} is only detected
in the literature in x-ray emission~\citep{sehgal2008}.
\citet{wittman2014} also weight the source galaxies according to the expected
lensing kernel and the photo-$z$ estimates. We make no use of photo-$z$ information other than
in the sample selection.
\begin{figure*}[!htb]
	\centerline{
		\includegraphics[width=0.47\textwidth]{kappa_Abell781.png}
    \includegraphics[width=0.47\textwidth]{snr_Abell781.png}
	}
  \centerline{
    \includegraphics[width=0.47\textwidth]{kappa_Abell781_WD.png}
    \includegraphics[width=0.47\textwidth]{snr_Abell781_WD.png}
  }

	\caption{Posterior mean lens field maps (left) and SNR maps (right) for a
  field centered on the galaxy cluster Abell 781. We use 6000 galaxies selected
  to have photometric redshifts larger than the known redshifts of the two primary
  clusters in this field of view ($z=0.296$ for Abell 781 and $z=0.43$ for a chance
  alignment of a second cluster)
  The top row of panels show the posterior mean map and SNR using our GP
  prior. The bottom row of panels show the algorithm of \citet{2012ApJ...747L..42D}
  applied to the same galaxies with the same uniform per-galaxy weighting
  as in the top row, but evaluated on a finer grid. The normalization of the
  convergence in the lower left panel is arbitrary~\citep[see][]{2012ApJ...747L..42D}.
  The white crosses indicate sub-clusters identified by~\citet{sehgal2008} with
  x-ray detections. There is no associated mass for the white cross second
  from the right in any published analyses of this system.
	}
	\label{fig:mass_map_Abell}
\end{figure*}

\begin{table}[!htb]
  \begin{center}
    \caption{\label{tab:prior_params_abell}Parameters for the Gaussian prior on the
      GP parameters for A781.}
    \begin{tabular}{lcc}
      \hline\hline
      Parameter & mean & std. dev. \\
      \hline
      $\ln\left(\lambda_{\rm GP}\right)$ & $\ln\left(10^{6}\right)$ & 0.1 \\
      $\ln\left(\ell^2\right)$           & $\ln\left(10^{-4}\right)$ & 0.5\\
      \hline\hline
    \end{tabular}
  \end{center}
\end{table}
In \autoref{fig:abell_lnp} we show the marginal log-likelihood for the A781 ellipticities in
the plane of the two GP parameters. Unlike in \autoref{fig:sigma_e_maps} the noise covariance
is now significant in defining the contours in \autoref{fig:abell_lnp} such that a large
GP precision and small GP correlation length is favored (indicating a sub-dominant
signal covariance). We therefore impose a Gaussian prior in the logarithm of the GP parameters
with parameters listed in \autoref{tab:prior_params_abell}. We infer maximum posterior values
of $\lambda_{\rm GP}=2.7\times10^{6}$ and $\ell^2=0.012$, which we use in the convergence inference
in \autoref{fig:mass_map_Abell}. However, as shown in \autoref{fig:abell_lnp}, the marginal
posterior for the GP parameters is only weakly peaked for this data set, and a range of GP
parameters would be acceptable for the mass map inference.
\begin{figure}[!htb]
    \centerline{
        \includegraphics[width=0.45\textwidth]{lnp_Abell781.png}
    }
    \caption{The marginal posterior distribution for the GP parameters with
    priors imposed to favor correlation lengths smaller than the expected cluster sizes.}
    \label{fig:abell_lnp}
\end{figure}

\section{Conclusions}
\label{sec:conclusions}

We have demonstrated a probabilistic model inference of lensing shear and convergence
using a Gaussian Process (GP) prior. Our method is an extension of previous Bayesian maximum entropy
mass mapping methods that is applicable in both cosmic shear (i.e., `field') and cluster
lensing regimes. We validated our algorithm using simulated Gaussian distributed
shear maps with comparisons to the input `truth' and by reconstructing the mass distribution
in a cluster field in the Deep Lens Survey where we compare with the results of a previously
published `aperture densometry' method.

Because our GP kernel for the lensing shear and convergence is derived from a consistent lens
potential, we find excellent separation of E and B modes in the posterior lens field maps.
Recently \citet{2016arXiv161003345B} presented an algorithm to project Wiener Filter maps into
`pure' E and B modes, removing ambiguity in the E/B decomposition from survey masks. Because we compute
Wiener Filter solutions for the lensing convergence from galaxy ellipticities, the method of
\citet{2016arXiv161003345B} is a simple extension of the algorithm we present in this paper.

We have also described an algorithm for optimization of the GP kernel parameters given a
measured galaxy ellipticity catalog or a cosmologically-informed prior based on simulation
of lensing fields. An interesting future extension of this work could include both galaxy
positions and ellipticities in the reconstruction of the gravitational potential.

Our algorithm is computationally challenging in the solution of linear systems of dimension equal
to the number of galaxies. Next generation surveys will have $10^{7}$--$10^{10}$ galaxies, requiring
both parallelization of the linear solver routines and optimizations such as sparse matrix
approximations via tapering~\citep{Kaufman2008} or Fast Fourier Transform (FFT)-based matrix
multiplications~\citep[exploiting the special structure of the isotropic cosmological covariances][]{padmanabhan++2003}.
As a next step in the validation and scaling of
our code we plan to apply the algorithms in this paper to the variable shear `branches' of the
{\sc GREAT3}\xspace challenge\footnote{\url{http://great3challenge.info/}}~\citep{Mandelbaum++2013}, which
contain $2\times 10^{5}$ galaxies per simulated field. Coincidentally, each of the five
DLS four square degree fields contains a comparable number of $\sim2\times10^5$ lensing source
galaxies.

For the algorithm we validated in this paper we assumed an approximate data vector of galaxy
ellipticity measurements rather than the more informative imaging pixel data for each galaxy, as we
described previously in \citet{mbi-theory}.
If we use galaxy samples rather than an ellipticity catalog then we have to go back to
\autoref{eq:marg_like} and replace the integral over the intrinsic ellipticities with a numerical
importance sampling formula from \citet{mbi-theory}. In this case we do not need to assume
conjugacy of the intrinisc ellipticity prior and the likelihood function, which is good because
the pixel-level likelihood will not be conjugate. The algorithms we presented in
\autoref{sub:cosmo_params} for cosmological parameter estimation are applicable, however, when using
either an approximate ellipticity data vector or the interim sampling of galaxy image model
parameters.

We can still derive a Wiener Filter from the product of the Gaussian intrinsic ellipticity
distribution (or DP base distribution) and the interim variable shear GP prior assuming a weak
shear approximation. But then we need to evaluate \autoref{eq:posterior_cov} and
\autoref{eq:posterior_mean} for every interim sample.
The parameters of $\mathsf{S}+\mathsf{N}$ include the GP parameters in $\mathsf{S}$ and the $\alpha$ intrinsic
ellipticity distribution parameter in $\mathsf{N}$. If we update either parameter we need to perform a
new matrix factorization and solve operation at every step (but these operations can be done once
for all interim samples of galaxy image model parameters in a given step). Exploring optimized
linear algebra approaches and effective sampling strategies for this more complete framework will
be a focus of future work.

We have not included any redshift information about the source galaxies in our model. However, the
use of such information (via photometric redshifts) is a critical component of the weak lensing
analyses for cosmic shear surveys as well as cluster mass reconstructions with optimized
signal-to-noise ratios. A simple extension of our work to include redshifts could be to impose
separable GP priors on the lens fields for galaxy samples binned in redshift (or redshift estimator).
The physical correlations between the lens fields inferred in each source bin can be modeled in
the hierarchical inference stage when marginalizing lens fields to infer cosmological parameters
\citep[as outlined in][]{mbi-theory}. A more thorough approach would be to include probabilistic
redshifts for each source galaxy and marginalize over each source redshift distribution as
part of the lens field inference.

\section*{Acknowledgments}
We thank Chris Paciorek for discussion about the
statistical framework for performing shear inference and the use of Gaussian Processes.
We thank Alex Malz for critical reviews of draft versions of this work.
Thanks to M. James Jee for providing the Deep Lens Survey shear catalog.
Part of this work performed under the auspices of the U.S. Department of Energy
by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.
Funding for this work was provided by LLNL Laboratory Directed Research and Development
grant 16-ERD-013.
This work was also supported by the Director, Office of Science, Office of Advanced
Scientific Computing Research, Applied Mathematics program of the U.S. Department of Energy under
Contract No.DE-AC02-05CH11231.
This research used resources of the National Energy Research Scientific
Computing Center (NERSC), a DOE Office of Science User Facility supported by
the Office of Science of the U.S. Department of Energy under Contract No.
DE-AC02-05CH11231.
This work uses a modified version of the public code \texttt{George} available at
\url{https://github.com/karenyyng/george}, which was forked from
\url{https://github.com/dfm/george}. We also made use of the \texttt{daft} package
(developed by Dan Foreman-Mackey, David W. Hogg, and contributors, and available at
\url{http://daft-pgm.org/}) for plotting probabilistic graphical models.

\bibliographystyle{apj3auth}

redpajama_set_name: RedPajamaArXiv


For all you Stevie Nicks and Fleetwood Mac lovers out there, this is the perfect keychain. Opaque black with gold lettering and lobster clasp. There are 2 magical numbers to choose from - 7 or 13, pick your favorite.


This looks deliciously easy. I was about to say it is a bit out of season, but seeing how you are in Hong Kong, you can get the tropical Asian fruits most of the year.
That looks super delicious! Yum!
I'm always keen on trying new pairings with guacamole! This looks beautiful and delicious!


text:
Titanic Rising
Titanic Rising Weyes Blood
Label: Sub Pop Records
Genre: Alternative
Subgenre: Indie Rock
Artist: Weyes Blood
I`m sorry!
Dear HIGHRESAUDIO Visitor,
due to territorial constraints and also different releases dates in each country you currently can`t purchase this album. We are updating our release dates twice a week. So, please feel free to check from time-to-time, if the album is available for your country.
We suggest, that you bookmark the album and use our Short List function.
Thank you for your understanding and patience.
Yours sincerely, HIGHRESAUDIO
1A Lot's Gonna Change04:21
2Andromeda04:40
3Everyday05:07
4Something to Believe04:45
5Titanic Rising01:36
6Movies05:53
7Mirror Forever05:05
8Wild Time06:09
9Picture Me Better03:41
10Nearer to Thee01:05
Info for Titanic Rising
The phantom zone, the parallax, the upside down—there is a rich cultural history of exploring in-between places. Through her latest, Titanic Rising, Weyes Blood (a.k.a. Natalie Mering) has, too, designed her own universe to soulfully navigate life's mysteries. Maneuvering through a space-time continuum, she intriguingly plays the role of melodic, sometimes melancholic, anthropologist.
Tellingly, Mering classifies Titanic Rising as the Kinks meet WWII or Bob Seger meets Enya. The latter captures the album's willful expansiveness ("You can tell there's not a guy pulling the strings in Enya's studio," she notes, admiringly). The former relays her imperative to connect with listeners. "The clarity of Bob Seger is unmistakable. I'm a big fan of conversational songwriting," she adds. "I just try to do that in a way that uses abstract imagery as well."
"An album is like a Rubik's Cube," she says. "Sometimes you get all the dimensions—the lyrics, the melody, the production—to line up. I try to be futuristic and ancient at once, which is a difficult alchemy. It's taken a lot of different tries to get it right." As concept-album as it may sound, it's also a devoted exercise in realism, albeit occasionally magical. Here, the throwback-cinema grandeur of "A Lot's Gonna Change" gracefully coexists with the otherworldly title track, an ominous instrumental.
Titanic Rising, written and recorded during the first half of 2018, is the culmination of three albums and years of touring: stronger chops and ballsier decisions. It's an achievement in transcendent vocals and levitating arrangements—one she could reach only by flying under the radar for so many years. "I used to want to belong," says the L.A. based musician. "I realized I had to forge my own path. Nobody was going to do that for me. That was liberating. I became a Joan of Arc solo musician."
The Weyes Blood frontwoman grew up singing in gospel and madrigal choirs. "Classical and Renaissance music really influenced me," says Mering, who first picked up a guitar at age 8. (Listen closely to Titanic Rising, and you'll also hear the jazz of Hoagy Carmichael mingle with the artful mysticism of Alejandro Jodorowsky and the monomyth of scholar Joseph Campbell.) "Something to Believe," a confessional that makes judicious use of the slide guitar, touches on that cosmological upbringing. "Belief is something all humans need. Shared myths are part of our psychology and survival," she says. "Now we have a weird mishmash of capitalism and movies and science. There have been moments where I felt very existential and lost."
As a kid, she filled that void with Titanic. (Yes, the movie.) "It was engineered for little girls and had its own mythology," she explains. Mering also noticed that the blockbuster romance actually offered a story about loss born of man's hubris. "It's so symbolic that The Titanic would crash into an iceberg, and now that iceberg is melting, sinking civilization." Today, this hubris also extends to the relentless adoption of technology, at the expense of both happiness and attention spans.
The track "Movies" marks another Titanic-related epiphany, "that movies had been brainwashing people and their ideas about romantic love." To that end, Mering has become an expert at deconstructing intimacy. Sweeping and string-laden, "Andromeda" seems engineered to fibrillate hearts. "It's about losing your interest in trying to be in love," she says. "Everybody is their own galaxy, their own separate entity. There is a feeling of needing to be saved, and that's a lot to ask of people." Its companion track, "Everyday," "is about the chaos of modern dating," she says, "the idea of sailing off onto your ships to nowhere to deal with all your baggage."
But Weyes Blood isn't one to stew. Her observations play out in an ethereal saunter: far more meditative than cynical. "I experience reality on a slower, more hypnotic level," she says. "I'm a more contemplative kind of writer." To Mering, listening and thinking are concurrent experiences. "There are complicated influences mixed in with more relatable nostalgic melodies," she says. "In my mind my music feels so big, a true production. I'm not a huge, popular artist, but I feel like one when I'm in the studio. But it's never taking away from the music. I'm just making a bigger space for myself."
And In The Darkness, Hearts Aglow
The phantom zone, the parallax, the upside down—there is a rich cultural history of exploring in-between places. Through her latest, Titanic Rising (out April 5, 2019 on Sub Pop Records), Weyes Blood (a.k.a. Natalie Mering) has, too, designed her own universe to soulfully navigate life's mysteries. Maneuvering through a space-time continuum, she intriguingly plays the role of melodic, sometimes melancholic, anthropologist.
The track "Movies" marks another Titanic-related epiphany, "that movies had been brainwashing people and their ideas about romantic love." To that end, Mering has become an expert at deconstructing intimacy. Sweeping and string-laden, "Andromeda" seems engineered to fibrillate hearts. "It's about losing your interest in trying to be in love," she says. "Everybody is their own galaxy, their own separate entity. There is a feeling of needing to be saved, and that's a lot to ask of people." Its companion track, "Everyday," "is about the chaos of modern dating," she says, "the idea of sailing off onto your ships to nowhere to deal with all your baggage." But Weyes Blood isn't one to stew. Her observations play out in an ethereal saunter: far more meditative than cynical. "I experience reality on a slower, more hypnotic level," she says. "I'm a more contemplative kind of writer." To Mering, listening and thinking are concurrent experiences. "There are complicated influences mixed in with more relatable nostalgic melodies," she says. "In my mind my music feels so big, a true production. I'm not a huge, popular artist, but I feel like one when I'm in the studio. But it's never taking away from the music. I'm just making a bigger space for myself."

redpajama_set_name: RedPajamaCommonCrawl


Meaning and definitions of version, translation in Tamil language for version with similar and opposite words. Also find spoken pronunciation of version in Tamil and in English language.
What version means in Tamil, version meaning in Tamil, version definition, examples and pronunciation of version in Tamil language.


text:
We are pleased to announce that on August 1st at Gigi's Restaurant, 2080 Cato Ave., State College, Pa., the Penn State Quarterback Club will host Penn State Football Fans to interact and discuss the 2019 season.
Joining us from Penn State Athletics will be Spider Caldwell and Scott Sidwell.
Cost: $25.00 per member or $30.00 per non member
Cash Bar – Note: Happy Hour will be from 5:00 PM to 6:00 PM; Regular Prices after 6:00 PM
Food Menu:
Crudite Platter
Spinach and Artichoke Dip
Hummus with Chips
Shrimp Cocktail
Chicken and Waffle Sliders with Hot Honey and House Made Pickles
Pecan Chicken with Bourbon Sauce
Garden Sliders (Fried Green Tomato, Pimento Cheese, Pickled Peppers, Field Greens)
Stuffed Red Potatoes
Mini Mud Pie
Payments will be received at the door by either Cash or Check.
Please let us know if you will be joining us by emailing your response to
president@statecollegeqbclub.com no later than July 19th.
Posted in Events, In the News
Offseason standout, Journey Brown
Arguably one of the deepest rooms for Penn State, coach Ja'Juan Seider's running back group, will be a staple of the Nittany Lion offense this season. One guy who stood out during the offseason, according to coaches and players within the program, is Pennsylvania native Journey Brown.
READ MORE ..
Posted in In the News
UNIVERSITY PARK, PA. – Fifteen Nittany Lion lettermen worked out for scouts at Holuba Hall Tuesday morning. Pro Day participants completed testing and on-field drills for 30 NFL teams, two CFL squads and one WWE/XFL representative. The participants were: Ryan Bates, Jake Cooper, Donnell Dix, Koa Farmer, Kevin Givens, Connor McGovern, Trace McSorley, Shareef Miller, Amani Oruwariye, Miles Sanders, Nick Scott, Johnathan Thomas, DeAndre Thompkins, Kyle Vasey and Chasz Wright.
The 2019 NFL Draft takes place in Nashville, Tennessee Thursday, April 25 through Saturday, April 27. Coverage of the draft will air on the ESPN family of networks and the NFL Network.
PLAYERS MENTIONED
#52 Ryan Bates
6′ 4″
#33 Jake Cooper
#95 Donnell Dix
5′ 10″
#7 Koa Farmer
#30 Kevin Givens
#66 Connor McGovern
#9 Trace McSorley
#48 Shareef Miller
#21 Amani Oruwariye
#24 Miles Sanders
McSorley and Gillikin Named CoSIDA Academic All-District
Junior punter Blake Gillikin (Smyrna, Ga.) and senior quarterback Trace McSorley (Ashburn, Va.) were named CoSIDA Academic All-District 2 for their accomplishments in the classroom and on the field. Gillikin collects all-district accolades for the second-straight season. The pair is now on the ballot for CoSIDA Academic All-America.
Gillikin has a perfect 4.0 cumulative GPA as a kinesiology major. He received the President's Award for freshmen in 2016 and was presented with The Evan Pugh Scholar Award, given to those juniors and seniors who are in the upper 0.5 percent of their respective classes, in May 2018. Gillikin is a two-time Academic All-Big Ten honoree.
No. 14 Penn State Wins at Rutgers, 20-7
PISCATAWAY, N.J. – A pair of touchdowns by freshman tight end Pat Freiermuth and three takeaways on defense led No. 14 Penn State past Rutgers, 20-7, in Big Ten football action Saturday afternoon in HighPoint.com Stadium.
Quarterback Trace McSorley went 17-of-37 passing for 183 yards with two touchdowns and an interception and rushed for 39 yards on seven carries for Penn State (8-3, 5-3 Big Ten). It was the 30th time in his career that he led the Nittany Lions to victory, setting the school record. In addition, he became just third quarterback in Big Ten history to surpass 100 career touchdowns responsible for, joining Drew Brees (Purdue) and J.T. Barrett (Ohio State) as the only Big Ten quarterbacks to accomplish the feat.
Windsor Named Big Ten Defensive Player of the Week
Junior defensive tackle Robert Windsor (Fond du Lac, Wis.) earned his first Big Ten Defensive Player of the Week award after his performance against his home-state Wisconsin Badgers. He is the second Nittany Lion to win defensive weekly honors this season (Yetur Gross-Matos vs. Iowa).
Windsor was a key cog in the Penn State defense that held the Badgers to their fewest points (10) since 2016 and fewest passing yards (60) since 2014 in a 22-10 win Saturday at Beaver Stadium. Windsor finished second on the team with a career-high-tying six tackles, including a personal-best 2.0 sacks. Windsor's previous personal mark for sacks was 1.5 against Illinois earlier this season. After allowing 79 yards on Wisconsin's opening possession, the Nittany Lion defense buckled down to allow just 190 yards the rest of the game. The 269 yards of total offense for Wisconsin is 178.2 yards below the Badgers' season average (447.2). The Penn State defensive line accounted for five sacks, led by Windsor's 2.0 quarterback hits. Entering the game, Wisconsin had only allowed 13 sacks all season.
Penn State-Rutgers Game to Kick at noon
UNIVERSITY PARK, Pa. – The No. 20-ranked Nittany Lions travel to Rutgers for their final road game of the regular season Saturday for a noon kick on Big Ten Network. In addition, the game will also be available on the Penn State Sports Network (10:30 a.m. airtime) and GoPSUsports.com.
Penn State is 26-2 all-time against Rutgers, including a 5-0 mark in away games. The last time the teams met in Piscataway, Penn State posted a 39-0 win in tough weather conditions in 2016. Last season, the Scarlet Knights took the early 6-0 lead, but Penn State reeled off 35 unanswered points to secure the Homecoming win. Trace McSorley was 16-for-21 for 214 yards and two touchdowns to lead the offense, while Jason Cabinda had 11 tackles, including one tackle for loss, and Brandon Smith made 10 stops to pace the defense.
No. 20 Penn State Stifles Wisconsin, 22-10
UNIVERSITY PARK, Pa. – No. 20 Penn State limited Wisconsin to its lowest point total (10 points) since 2016 and Miles Sanders rushed for 159 yards and a touchdown as the Nittany Lions earned a 22-10 win over the Badgers Saturday afternoon at Beaver Stadium.
Wisconsin (6-4, 4-3 Big Ten) was first on the scoreboard with a touchdown on its first possession, but Penn State (7-3, 4-3 Big Ten) outscored the Badgers 22-3 over the remainder of the game for its 10th-straight win over a Big Ten West Division opponent and fourth-straight win over Wisconsin. Penn State clinched its 14th consecutive winning season with the win, becoming one of just four teams to accomplish the feat thus far this season.


ProMetric® Y is a series of high-performance imaging photometers designed for production volume testing of displays, smartphones, tablets, notebooks, keyboards, and lighting products.
The performance and appearance that consumers expect from products such as flat-panel displays, car interiors, and LED lighting are continuously advancing. Designers and manufacturers are pushed to improve the user experience. Test & Measurement instruments used in development and production must keep up with the latest technology.
ProMetric Y is built around scientific-grade CCD sensors that are thermoelectically cooled ranging in resolution from 2 to 29 megapixels. to provide accurate, repeatable 12-bit measurements. The ProMetric Y2 uses a scientific-grade 2-megapixel (1600 x 1200) CCD sensor that is thermoelectrically cooled to provide accurate, repeatable 12-bit measurements. The Y16 provides more resolution with a 16-megapixel (4896 x 3264) CCD sensor. The Y29 provides even higher resolution with a 29-megapixel (6576 x 4384) CCD sensor. Each ProMetric Y supports high-speed USB and Ethernet communications.