Strong Gravitational Lens Modeling with Spatially Variant Point Spread Functions

1.7 MB·English

by  Adam Rogers

Publishedin ApJ,743:68,Dec. 10,2011. Submitted2011July21;accepted2011Sept. 30
PreprinttypesetusingLATEXstyleemulateapjv.8/13/10
STRONG GRAVITATIONAL LENS MODELING WITH SPATIALLY VARIANT POINT SPREAD FUNCTIONS
Adam Rogers & Jason D. Fiege
DepartmentofPhysicsandAstronomy,TheUniversityofManitoba,Winnipeg,Manitoba,R3T-2N2,Canada
Published inApJ, 743:68, Dec. 10, 2011. Submitted 2011 July 21; accepted 2011 Sept. 30
ABSTRACT
Astronomical instruments generally possess spatially variant point-spread functions, which deter-
2 minetheamountbywhichanimagepixelisblurredasafunctionofposition. Severaltechniqueshave
1 been devised to handle this variability in the context of the standard image deconvolution problem.
0 We have developed an iterative gravitational lens modeling code called Mirage that determines the
2 parameters of pixelated source intensity distributions for a given lens model. We are able to include
n the effects of spatially variant point-spread functions using the iterative procedures in this lensing
a code. In this paper, we discuss the methods to include spatially variant blurring effects and test the
J results of the algorithm in the context of gravitational lens modeling problems.

0 Subject headings: gravitationallensing: strong — methods: numerical
1

] 1. INTRODUCTION σ,andgivenasourceintensitydistribution,theresulting
O χ2 statistic between the data and model is
Modeling gravitational lens systems is fundamentally
C a twostep process. A successfulmodel requiresanouter 2
h. olepntsimmizoadteiol,napnrodceadunreestteoddiospcotivmeriztahteiopnarsatmepetteorsdoefttehre- χ2 = di− jfijsj (2)
p mine the structure of the lensed source object. This in- i (cid:16) Pσi2 (cid:17)
- X
o ner source optimization step must incorporate the effect
By requiring derivatives of this equation with respect to
r of the point spread function (PSF). Observed data are
t the elements of the source vector s vanish, we find that
s blurred by the PSF, and the presence of noise compli- theoptimalsourcepixelintensitiessatisfyaleast-squares
a catesthedeconvolutionprocessingeneral(Hansen et al.
equation:
[ 2o0f0a6u)t.hoTrhsiusspinrgobslpeamtiahlalysibneveanriawnetllPsStFudsi(eWdabryreanv&arDieytye FTFs=FTdˆ (3)
1
(2003); Koopmans (2005); Suyu et al. (2006)).
v where we have absorbed factors of σ into the matrix
i
7 Toincludeblurringandlensingeffects,wedescribeop- F = f /σ and data vector dˆ = d /σ . The details of
2 erationsonimagesbytheapplicationoflinearoperators. thijisderiivjatiioncanbefoundinWi arrein &i Dye(2003)and
1 We use “flattened” images to facilitate this notation, in
Koopmans (2005), where the system is solved by direct
2 which each column of the image is stacked upon the
matrix inversion with regularization. This least squares
. next, forming a vector. Our preferred scheme for solv-
1 form is commonly found in the context of large scale
ing the lens modeling problem is a modification of the
0 imagedeconvolutionproblems,whicharetypicallysolved
semilinear method of Warren & Dye (2003), as outlined
2 by iterative methods (Golub & Reinsch (1970); Bj¨orck
inRogers & Fiege(2011). Tobegin,wedefineimageand
1 (1996);Hansen(2010)). Notethatthesemilinearmethod
sourcecoordinatesystemsthatareconnectedbythethin
: providesa technique for solving the linear parametersof
v lens equation:
thelensedsystemonly(thesourceintensitydistribution)
i β =θ−α(θ), (1)
X inan“innerloop”,whilethenonlinearparametersofthe
r where θ and β are the image and source coordinates re- lensmassdistributionmustbesolvedinaseparate“outer
a spectively, and α(θ) is the deflection angle determined loop” optimization step (Rogers & Fiege 2011).
by the gravitational potential of the lens density dis- Spatial dependence of the PSF is not considered in
tribution. The lens equation is a nonlinear equation most conventional deconvolution problems. This simpli-
that maps pixels from the image to the source plane fiestheconstructionoftheblurringmatrixB,sinceonly
(Schneider et al. (1992); Petters et al. (2001)). onePSFistakenintoaccount. However,itiswellknown
The source plane pixels are represented as a vector s, that the PSF cannot always be treated as constant over
with elements that can be varied independently to ac- an image in cases of astronomicalinterest. For example,
countforthedetailsoftheunknownsourceintensitydis- spatially variant PSFs have been studied in the context
tribution. ThedetailsofthePSFareencodedintheblur- ofadaptiveoptics(Lauer(2002);Gilles et al.(2002))and
ring matrix B and gravitational lens effects, described thePSFofastronomicalinstruments,suchastheHubble
byEquation1,areincludedinthe lensing matrixL. We AdvancedCameraforSurveys,canbeextremelyposition
then define the total lens matrix f = BL and data vec- dependent (Bandara et al. 2009). Several schemes have
tor d. Denoting the standard deviation of the noise as been designed to deal with this variability (Boden et al.
(1995); Biretta(1994); Adorf(1994); Lauer(2002)). De-
[email protected] scribing a spatially variant PSF is much more compli-
2 Rogers & Fiege
cated than for the invariant case, since each row of the of the image. The image is then divided into a square
blurring matrix B will be derived from a unique PSF in grid, where the PSF is assumed constant in each region.
general. The position of a pixel in the image determines These image regions and the PSFs are then padded to
the amount by which it is blurred. match in size. The two-dimensional fast Fourier trans-
Weillustratethe effectofspatiallydependentPSFson form(FFT)isusedtocalculatetheresultantblurredim-
gravitationally lensed images in Figure 1. Consider the ageregionsindependently,resultinginaneffectivepiece-
lensing effect produced by a singular isothermal sphere wiseconvolution. Bysubstitutingefficientalgorithmsfor
(SIS), which has three parameters: velocity dispersion the explicit matrix andmatrix-transposemultiplications
σ and lens center (x,y). The deflection angle due to a in Equation 3, the least squares form of the problem is
v
SIS lens has a simple analytical form most conveniently preservedand the system can be solved efficiently.
described in standard polar coordinates (r,ω): Inprinciplespatiallyvariantblurringcanbedescribed
by a blurring matrix compatible with the semilinear
σ 2 D
v ls method. However, in practice there are several prob-
α(r)=4π , (4)
c D lems with the matrix approach. First, the size of the
os
(cid:16) (cid:17) blurring matrix is N × N , so the matrix quickly
whereD andD aretheangulardistancesbetweenlens pix pix
ls os becomes large as the image resolution is increased. Sec-
andsourceandobserverandsourcerespectively,andcis
ond, since the PSFs vary over regions of the image, it is
thespeedoflight. WemodeltheblurringinFigure1with
possiblethatB maycontainalargenumberofsmallbut
σ = 265 km s−1 and use source redshift z = 1.5 and
v s non-zeroentries,particularlyforlarge,complicatedPSFs
lens plane redshift z =0.12. This model was calculated
l that are not well approximated by Gaussians or other
using cosmological parameters H = 70 km s−1 Mpc−1,
0 simple analytical functions. This complicates the opti-
Ω0 =0.3andΛ0 =0.7whichweadoptfortheremainder mizationbecauseM=FTF mustbeinvertedinthesemi-
ofthisstudy. Thesourceiscomprisedofasetofcircular
linear scheme. It is generally required that M is sparse
disksin the sourceplane asshownin the left-handpanel
in order to store and invert this large matrix. The spar-
ofFigure1andthegravitationallylensedimageisshown
sity requirement helps to reduce computation time and
in the center panel. The lensed image is then blurredby
reduces the amplification of noise in the reconstructed
aspatiallyvariantPSFandisshownintherightpanelof
source. In practice the semilinear method requires regu-
Figure1. The distortionusedto createthis imagevaries
larization to control the amount of noise present in the
from a delta function in the lower-left corner (negligible
solutionof Equation3. The details andeffects of several
blur) to a Gaussian with standard deviation σ = 6.0
g distinct regularizationmethods used with the semilinear
pixels in the upper rightcorner. Each PSF is defined on
method were studied in detail by Suyu et al. (2006).
anarbitrary33×33gridandisnormalizedtounitysum.
Our previous work (Rogers & Fiege 2011), compared
The source and image plane size are 240×240 pixels.
the semilinear method with severaliterative methods to
Unlike constant PSFs, spatially variant PSFs cannot
solve the least-squares problem (Equation 3). Iterative
be described by a simple convolution operation. For-
schemeshave the advantagethat time is savedby avoid-
tunately, numerical methods have been devised to han-
ingtheexplicitconstructionofthelensandblurringma-
dlethem, including sectioningmethods(Trussel & Fogel
trices. This is done using direct interpolation on the
1992), which deconvolve each PSF independently and
sourceplaneunderthe effectofthe lensequation(Equa-
forms the source from the sum of the reconstructions.
tion 1).
Nagy & O’Leary (1998) devised a clever method to
Rogers & Fiege (2011) used the Qubist Optimization
model the effects of spatially variant PSFs within the
Toolbox (Fiege 2010) to map the χ2 surface over the
frameworkofthestandardimagedeconvolutionproblem.
space of the nonlinear lens parameters using the Ferret
Thisapproachdiffersfromsectioningmethodsinthatthe
GeneticAlgorithm(GA)andLocustParticleSwarmOp-
separatePSFsareusedtobuildanapproximationtothe
timizer (PSO).Since this mapping requiresa largenum-
blurred image of a given source, and a single iterative
ber of function evaluations (≈ 105) over the course of a
deconvolutionoperationisneededtosolveforthesource
run, speed is of the essence when choosing an inner loop
intensity distribution. The method was implemented in
optimization to determine the source plane parameters.
Nagy et al.(2002)andrepresentsthespatialdependence
Using the techniques introduced by Nagy et al.(2002)
of the PSF as a summation of piecewise blurring ma-
asa foundation, we haveadded the capabilityto include
trices, each of which applies over a limited area of the
spatiallyvariantPSFstoourgravitationallensmodeling
image. In this study, we use the method of Nagy et al.
codeusingpiecewiseconstantPSFs. Thisnewcapability
(2002) to incorporate spatially variant blurring into our
is the subject of the current exploration.
gravitational lens modeling code. We briefly review the
method here and discuss the procedure in detail in the
2. ASMALL-SCALETEST
Appendix.
Toincludetheeffectsofspatiallyvariantblurs,theim- In this section, we provide an example of modeling an
age of the unblurred lensed source is padded to enforce extended source under the effects of a spatially variant
a boundary condition(Hansen et al.2006). We focus on PSF.Wegeneratethelensedimageofananalyticalfunc-
theuseofreflexiveboundaryconditions,inwhichtheim- tion that describes a spiral source intensity distribution,
ageispaddedbysymmetricreflectionsofitself. Reflexive according to the equation:
boundary conditions tend to reduce ringing artifacts if a
S
significant amount of structure is located near the edges S(r,ω)= 0 exp −2sin2 ω−ω −τr2 , (5)
r2+r2 0
c
(cid:2) (cid:0) (cid:1)(cid:3)
Strong Gravitational Lens Modeling with Spatially Variant PSFs 3
where S is the maximum brightness in arbitrary units between the model solution and the true solution im-
0
and core radius r . The tightness of the arms about the proves until a minimum is reached and then begins to
c
centralbulgeiscontrolledbyτ,andω controlstheorien- increase. This is due to the properties of the local opti-
0
tation of the spiral, in standardpolar coordinates (r,ω). mizer used to determine the optimal source, and arises
This artificial “galaxy”, originally described by Bonnet in the deconvolution step due to noise in the observed
(1995), serves as a convenient test pattern. To draw image. Regularization methods are generally used to
comparisons between our results for spatially invariant control the increase of noise in the reconstructed source
PSFs (Rogers & Fiege 2011), we will again make use of foundbythe semilinearmethod(Suyu et al.2006). Sev-
a Singular Isothermal Ellipse (SIE; Keeton & Kochanek eral optimization methods have been applied to prob-
(1998)) with deflection angle components: lems with spatially variant blur including Landweber
iteration (Nocedal & Wright (1999); Fish et al. (1996);
bq x 1−q2 Trussel & Hunt(1978)), Richardson-Lucydeconvolution
αx = tan−1 (6) (Faisal et al.1995),andLanczos-Tikhonovhybridmeth-
1−q2 pψ+s ! ods (Chung et al. 2008) in the context of the standard
p imagedeconvolutionproblem. FollowingRogers & Fiege
α = bq tanh−1 y 1−q2 , (7) (2011), we focus on the conjugate gradient method for
y 1−q2 ψp+q2s ! least-squaresproblems (CGLS) and the steepest descent
method (SD). Figure 6 shows the convergencehistory of
with ψ2 = q2(sp2 +x2)+y2 and q = (1−ǫ)/(1+ǫ), theSDalgorithm. AsintheinvariantPSFcasediscussed
and b is the equivalent Einstein radius when q =1. The in Rogers & Fiege (2011), the SD solution converges
parameter b is related to the velocitypdispersion σ by more slowly than CGLS and therefore it is less sensitive
v
Equation 4. tothestoppingcriteria. Whenusinganiterativemethod
Theparametersusedinthistestarevelocitydispersion forlocaloptimization,thenumberofiterationsitselfacts
σ = 265 km s−1 with z = 0.3 and z = 1.05 giving an as a regularization parameter. The optimal stopping it-
v d s
equivalent Einstein ring of b = 1.32 arcsec, ellipticity erationoftheselocaloptimizersisattheminimumofthe
ǫ = 0.35, lens center (x,y) = (0.11,0), core size s, and relativeerrorcurveforagivensetoflensparametersand
orientation angle θ = π/4 measured counterclockwise PSF tiling. This critical iteration represents a balance
L
from the right of the image. We set s=0, resulting in a between the reduced image χ2 and the amount of regu-
singular mass distribution. larization used (Press et al. 2007). Established methods
We used Equation 1 to form the lensed image of the exist to determine this critical iteration, including the
source (Equation 5) using the SIE deflection angle for- L-Curve criterion (Hansen & O’Leary 1993) and Gener-
mulae. We generated a 20×20 grid of spatially variant alized Cross Validation (Golub et al. 1979). In previous
Gaussian PSFs where each PSF is defined by a 33×33 work(Rogers & Fiege2011)wemadeuseoftheL-Curve
pixel mesh and has an FWHM ranging from 2.35 to 4.8 criterionbut Generalized Cross Validation is also imple-
pixels,showninFigure 2. This gridofPSFs wasusedto mented in our software.
blurthe gravitationallylensedimageandadditiveGaus- We find that the execution time of the problem in-
sian white noise with standard deviation σ = 1.05 was cluding a spatially variant PSF increases approximately
g
added after the blurring operation, resulting in the arti- linearly with the number of separate PSFs used in the
ficialdatashowninFigure 3. We define the peaksignal- inversionasshowninFigure7. Thissuggeststhatsignif-
to-noise ratio (PSNR) as icantgainscouldbe madeinthe efficiencyofthe routine
byparallelizingtheimplementation,sinceeachimagere-
I
max gion is independent. By splitting up the problem over
PSNR= , (8)
σ several processors, the runtime for very large PSF grids
g
can become feasible.
giving PSNR=105.83. To illustrate the effect of vary-
ingthenumberofPSFs,wemodelthedatausingsmaller
gridsof3×3,5×5,7×7,10×10,and20×20PSFs. As 3. ALARGE-SCALETEST
shownin Figure 4, the best reconstructionwith the low- To demonstrate the code in operation on a large scale
est reduced χ2 is obtained using a grid of 20×20 PSFs, problem, we simulate the lensing effect of the mass dis-
which is the same number used to generate the data. tribution of a galaxy cluster on a portion of the Hub-
This source and corresponding model image after 20 it- ble deep field using an elliptical potential. This test is
erations are also shown in Figure 3. The 3×3 and 5×5 intended as a demonstration of the feasibility and effi-
image residuals show significant structure, which is not ciency of our method on a problem that would be diffi-
present in the finer approximations. The residuals using cult using the semilinear method while including a spa-
a grid of 20×20 PSFs appear featureless. This demon- tially dependent PSF. Problems of this size are realis-
strates the improvement in image reconstruction as we tic for a number of practical modeling situations. For
include successively more information characterizingthe example, Alard (2009) has modeled the lensed system
blur. SL2SJ021408-053532,whichproducesa setoflargearcs.
Figure 5 shows the relative error between the model Thissystemhasalensthatiscomprisedofasmallgroup
source and the true solution as a function of iteration. of six galaxies. Due to the large size of the lensed arcs,
For all PSF grid sizes, we find that the solutions display the scope of the source modeling prohibited the direct
semi-convergence behavior such that the relative error application of the semilinear method.
4 Rogers & Fiege
We form the lensed image of a portion of the Hubble terscouldbe foundusingglobaloptimizationmethods if
deep field (Williams et al. 1996) by applying the ellipti- one ofthe following strategieswereemployed: (1) a low-
calpotentialofBlandford & Kochanek(1987)whichwas resolutionapproximationtothedatacouldbeusedearly
used by Link & Pierce (1998) to simulate the lens effect during the lens parameter optimization, with successive
of the dark matter distribution of galaxy clusters. This refinement occurring later during the run; (2) a global
potential function is given by optimizercouldbe usedtoroughlyapproximatethe lens
parameters,shiftingtoafasterlocaloptimizationscheme
ψ(x,y)= b2(1−q) s2+(1+ǫc)x2+2ǫsxy+(1−ǫc)y2 q, otenrcsepsaocleu;toiorn(s3)argelolboacalloipzetdimtiozaatisomnaclolurledgiboenuosfedpaforarmthee-
2q
entire problem making use of large-scale parallelization.
(cid:2) (9(cid:3))
which results in deflection angle α(θ) = ∇ψ(θ). The 4. CONCLUSION
elliptical potential depends on seven parameters: b is
Wehavedevelopedamethodtoincludetheeffectsofa
the equivalent Einstein radius in the limit of vanishing
spatially variantPSFin gravitationallens modeling. In-
core radius s, ellipticity ǫ, and power law index q, where
0 ≤ q ≤ 0.5. The position angle of the lens φ deter- cluding these effects in the standard semilinear method
wouldbedifficultduetothecomplicatedblurringmatrix
mines the functions ǫ = ǫcosφ and ǫ = ǫsinφ. We
c s
required. Thesecomplicationscanbeovercomeeasilyby
use the Einstein radius b=9, power law index q =0.25,
incorporatingthemethodofNagy et al.(2002). Ourap-
φ = π/4, position (x,y) = (0,0) and s = 0.5. The lens
proachcanaccommodate largelensing problems like the
and source redshifts are z = 0.12 and z = 1.5, re-
d s
casestudiedbyAlard(2009),whichlimitstheapplicabil-
spectively. We used an array of 25 PSFs arranged on
a 5× 5 grid to blur the image. This set of PSFs has ity of the direct semilinear approach. Techniques to in-
cludetheeffectsofspatiallyvariantPSFsareimportant,
been used to test image restoration schemes for Hub-
as the response varies over the detector area for many
ble Space Telescope (HST) images and represents the
astronomical instruments. Our algorithm allows this ef-
spatially variant nature of the aberrations affecting the
fect to be included in lensing problems, thus improving
HST before it was repaired (Katsaggelos et al. (1994);
Nagy & O’Leary(1998)). ThesizeofeachPSFis60×60 the qualityofreconstructionswhenthe variabilityofthe
PSF is significant. The CGLS and SD algorithms allow
pixels, and the source and image plane used to generate
our lensed image are 800×800 pixel2. Gaussian white a regularizedinversionto be found quickly by truncated
iteration.
noise was added with standard deviation σ =1.37, giv-
g
ing the image PSNR=138.4.
5. ACKNOWLEDGEMENTS
Theimageafter100iterationsisshowninFigure8,and
a reducedχ2 =0.995was found. The system was solved A.R. acknowledges NSERC for funding this research,
using the CGLS algorithm with all 25 PSFs using the and J.F. acknowledges funding from an NSERC Discov-
lens parameters defined above. The model took approx- ery Grant. The authors also thank the anonymous ref-
imately 7 minutes to solve using a single 2.4 GHz CPU eree, whose comments helped to improve the flow of the
core. An approximation to the nonlinear lens parame- paper.
APPENDIX
SPATIALLY VARIANTBLURRING EFFECTS
TodescribeblurringbyaspatiallyvariantPSFwefirstpresentanefficientmethodusingtwo-dimensionalFFTs. We
then show how to treat the problem in terms of blurring matrices and flattened image vectors. See Nagy & O’Leary
(1998) for more details on the approachand Nagy et al. (2002) for a MATLAB implementation.
Consider anN×N gridof independent PSFs P and split the unknown blurredimage Y into regionsY , eachof
ij ij
size k×k:
Y Y ··· Y
11 12 1N
Y Y ··· Y
21 22 2N
Y = ... ... ... ... (A1)
Y Y ··· Y
N1 N2 NN
Each of these blocks will be affected by an independent PSF. Suppose that the size of each PSF is (r+1)×(r+1)
with r even, and let the unblurred N ×N image be represented by X.
Let us define a set of “mask” matrices w . In the case of piecewise constant PSFs, these masks are the same size
ij
as the unblurred image and are comprised of 0 entries everywhere except for the k×k block at position (i,j), where
the entries of w are set to 1.
ij
TofindthecomponentsofagivenregionweconvolveX withthecorrespondingPSFP ,followedbyanelement-wise
ij
multiplication by the mask w . The non-zero elements of this product give Y . Proceeding in this way we build up
ij ij
the blurred image block by block:
N N
Y = w ◦(P ∗X), (A2)
ij ij ij
i=1j=1
XX
Strong Gravitational Lens Modeling with Spatially Variant PSFs 5
where the symbol “◦” represents element-wise multiplication and symbol “∗” is the convolution operation. Note that
each term in the sum is determined by the convolution of the entire image X with the appropriate PSF before the
mask is applied. This is crucial to ensure that “seams” will not be visible between regions in the blurred image Y.
Ingeneral,it is possible to speedup this routine by calculatingY directly. Considersplitting the unblurredimage
ij
into regions Xk where the superscriptdenotes the size ofthe block,in this case k×k. In order to avoidartifacts and
ij
keepthe correctintensity near the edges of this block after convolution,we include a number of neighboring rowsand
columns on each side of Xk. The width of this border is set by the size of the PSF, r/2, with regions on the image
ij
boundarypaddedtoenforcethe boundaryconditionsdiscussedin Section1. Theseextendedregionsarethendenoted
X(r+k). ThePSFsarepaddedtomatchtheextendedregionsinsize,resultinginP(r+k). Theblurredextendedregion
ij ij
is found by the convolution
Y(r+k) = P(r+k)∗X(r+k) . (A3)
ij ij ij
(cid:16) (cid:17)
The central k×k block of this product is clipped out and placed in the (i,j) position of Y. The process is repeated
until the entire blurred image is filled in. Time is saved working with extended regions and padded PSFs since we
only need to calculate the convolution over the (r+k)×(r+k) block for each PSF rather than the entire image as
in Equation A2, and the construction of masks is not needed. The convolutions can be carried out efficiently with
two-dimensional FFTs.
Thebasicprocedurecanalsobedescribedbyananalogousmatrix-vectoroperation. ToexpressthesuminEquation
A2 in terms of matrix multiplication, we define the unblurredflattened image as a vector x, and the flattened blurred
image as y. We build a set of N2 blurring matrices to describe the effect of each PSF on x, which we denote as B .
ij
Themaskmatricesw areusedtoconstructanalogousweightingmatricesD . ThesematricesareofsizeN ×N ,
ij ij pix pix
where N is the number of pixels in the image, identical to the size of the blurring matrices B . The total blurring
pix ij
matrix B is then written as a weighted sum of blurring matrices B , B ,...,B .
11 12 NN
N N
B = D B . (A4)
ij ij
i=1j=1
XX
Theblurredimageisthenfoundbyamatrixmultiplicationy=Bx. TheweightingmatricesD havethemthdiagonal
ij
entryequalto1providedthatimagepixelmisinregion(i,j),andallotherelements0. Theweightingmatricessatisfy
N N D =I where I is the N ×N identity. We adopt the use of piecewise constant PSFs but in general
i=1 j=1 ij pix pix
it is possible to include higher order interpolation schemes between PSFs using the weighting matrices. The case of
P P
linear interpolationin solvingsystems with spatially variantblur has been studied by Nagy & O’Leary(1998), but its
inclusion complicates the procedure and did not provide a significant improvement to the quality of the solution and
increased computation times (Nagy et al. 2002).
REFERENCES
Adorf,H.M.,1994,inTheRestorationofHSTImages and Golub,G.H.,Heath,M.,&Wahba, G.1979,Technometrics, 21,
SpectraII,ed.R.J.Hanisch&R.L.White(Baltimore,MD: 2,215-223
SpaceTelescopeScienceInstitute), 72 Golub,G.H.,&Reinsch,C.1970,Numer.Math.,14,403
Alard,C.2009,A&A,506,609 Hansen,P.C.2010, DiscreteInverseProblems:Insightand
Bandara,K.,Crampton,D.,&Simard,L.2009,ApJ,704,1135 Algorithms(Philadelphia,PA:SIAMpublishers)
Biretta,J.1994inTheRestorationofHSTImagesandSpectra Hansen,P.C.,Nagy,J.G.,&O’Leary,D.P.2006,Deblurring
II,ed.R.J.Hanisch&R.L.White(Baltimore,MD:Space Images:Matrices,SpectraandFiltering(Philadelphia,PA:
TelescopeScienceInstitute), 72 SIAMpublishers)
Bjo¨rck,˚A.1996,NumericalMethods forLeastSquaresProblems, Hansen,P.C.,&O’Leary,D.P.1993,SIAMJ.Sci.Comput.,14,
(Philadelphia,PA:SIAMPublishers) 1487
Blandford,R.D.,&Kochanek, C.S.,1987,ApJ,321,658 Katsaggelos,A.K.,Kang,M.G.,&Banham,M.R.1994,inThe
Boden,A.F.,Redding,D.C.,Hanisch,R.J.,&Mo,J.,1995,J. RestorationofHSTImagesandSpectraII,ed.R.J.Hanisch&
Opt.Soc.Am.A,13,1537 R.L.White(Baltimore,MD:SpaceTelescopeScience
Bonnet,H.1995, PhDthesis,L’Univ.PaulSabatierdeToulouse Institute), 3.
Chung,J.,Nagy,J.G.,&O’Leary,D.P.2008,Elec.Trans. Keeton,C.R.&Kochanek, C.S.1998,ApJ,495,157
Numer.Anal.,28,149 Koopmans,L.V.E.2005,MNRAS,363,1136
Faisal,M.,Lanterman,A.D.,Snyder,D.L.,&White,R.L.1995, Lauer,T.2002,Proc.SPIE,4847, 167
J.Opt.Soc.Am.A,12,2593 Link,R.,&Pierce,M.J.,1987,ApJ,502,63
Fiege,J.D.,2010,QubistUsersGuide:Optimization,Data Nagy,J.G.&O’Leary,D.P.1998,SIAMJ.Sci.Comput.19,
Modeling,andVisualizationwiththeQubistOptimization 1063
ToolboxforMATLAB(Winnipeg:nQubeTechnical Nagy,J.G.,Palmer,K.M.,&Perrone,L.2002,Numer.
Computing) Algorithms,36,73
Fish,D.A.,Grochmalicki,J.,&Pike,E.R.1996,J.Opt.Soc. Nocedal,J.&Wright,S.J.1999,NumericalOptimization(New
Am.A,13,1 York:Springer)
Gilles,L.,Vogel,C.R.&Bardsley,J.M.2002,InverseProbl.18, Petters,A.O.,Levine,H.&Wambsganns, J.2001,Singularity
237 TheoryandGravitational Lensing(Boston, MA:Birkh¨auser)
6 Rogers & Fiege
Press,W.H.,Teukolsky, S.A.,Vetterling,W.T.,&Flannery,B.P. Trussell,H.J.&Fogel,S.1992,IEEETrans.ImageProc.1,123
2007, NumericalRecipes:TheArtofScientificComputing(3rd Trussell,H.J.&Hunt,B.R.1978, IEEETrans.Acoust.Speech,
ed.;NewYork:CambridgeUniv.Press) SignalProcessing,26,608
Rogers,A.&Fiege,J.D.2011,ApJ,727,80 Warren,S.J.&Dye,S.2003,ApJ,590,673
Schneider,P.,Ehlers,J.,&Falco,E.E.1992,Gravitational Williams,R.E.,Blacker,B.,Dickinson,M.,etal.1996,AJ,112,
Lenses (Berlin:Springer) 1335
Suyu,S.H.,Marshall,P.J.,Hobson,M.P.,&Blandford,R.D.
2006, MNRAS,371,983
Strong Gravitational Lens Modeling with Spatially Variant PSFs 7
Source Plane Unblurred Image Blurred Image
2 2
1
y 0 y 0 y 0
β θ θ
−1
−2 −2
−1 0 1 −2 0 2 −2 0 2
β θ θ
x x x
Fig.1.—Exampleofspatiallyvariantblurring. Left: asetof regulardiskswithradius0.268tilethe sourceplane. Center: thecircular
disksareseenunderthelensingeffectofaSingularIsothermalSphere(SIS)lensmodel. TheSISdistortsthebackgroundcirclesintoarcs,
andthediskatthecenteroftheSISbecomesacompletering. Right: thesamediskpatternundertheeffectoftheSISlens,withaspatially
variantPSFblurringtheobservation. TheblurisdescribedbyadeltafunctioninthelowerlefthandcornertoaGaussianwithstandard
deviationσg =6.0pixelsintheupperrightcorner,introducingasignificantblur.
8 Rogers & Fiege
Fig.2.—Gridof PSFs usedinFigure3. The PSFs varyfromaGaussian of FWHMof 2.35 pixelsinthe lower-leftcorner producing a
modestblurtoaGaussianwithFWHM4.75pixelsintheupperrightcorner.
Strong Gravitational Lens Modeling with Spatially Variant PSFs 9
Observation Model Image
1 1
y 0 y 0
θ θ
−1 −1
−1 0 1 −1 0 1
θ θ
x x
Source Model Source
0.4 0.4
y 0 y 0
β β
−0.4 −0.4
−0.4 0 0.4 −0.4 0 0.4
β β
x x
Fig.3.— Top left: artificial data on a 120×120 grid. Bottom left: artificial source on a 50×50 grid. Top right: model observation.
Bottomright: modelsource. Theresultsafter19iterationsareshown. Notethepresenceofreconstructed noiseinthesource. Themodel
hasareducedχ2=0.998.
10 Rogers & Fiege
3x3 5x5
7x7 10x10
Reduced χ2 vs N
20x20 psf
1.1
1.05
1
0 200 400
Fig.4.—Imageresidualsfora3×3,5×5,7×7,10×10,and20×20PSFgridsafter19CGLSiterations. ForasmallnumberofPSFs
thereisasignificantamountofresidualstructure,buttheseartifactsarereducedasthegridofPSFsisenlarged. Thereducedχ2 isshown
asafunctionofthenumberofPSFs(Npsf)usedintheinversion.