National Energy Research Scientific Computing Center (NERSC) Data Seminars Series, 19 Apr 2022

Previous Meeting Next Meeting

⏯

youtube image

►

From YouTube: 2022-04-19 - Huang, Gu - Discovering and Modeling Strong Gravitational Lenses at NERSC

Description

NERSC Data Seminars Series: https://github.com/NERSC/data-seminars

Title:
Discovering and Modeling Strong Gravitational Lenses with Cori and Perlmutter at NERSC

Speakers:
Xiaosheng Huang, USF
Andi Gu, UCB

Abstract:
We have discovered over 1500 new strong lens candidates in the Dark Energy Spectroscopic Instrument (DESI) Legacy Imaging Surveys with residual neural networks using NERSC resources. Follow-up observations are underway. Our Hubble Space Telescope program has confirmed all 51 observed candides. DESI observations have confirmed more systems spectroscopically. Preliminary results from our latest search will increase the number of lens candidates to over 3000. We have also developed GIGA-Lens: a gradient-informed, GPU-accelerated Bayesian framework, implemented in TensorFlow and JAX. All components of this framework (optimization, variational inference, HMC) take advantage of gradient information through autodiff and parallelization on GPUs. Running on one Perlmutter A100 GPU node, we achieve 1-2 orders of magnitude speedup compared to existing codes. The robustness, speed, and scalability offered by this framework make it possible to model the large number of strong lenses found in DESI, and O(10^5) lenses expected to be discovered in upcoming large-scale surveys, such as the LSST.

Bios:
Xiaosheng Huang received his PhD from UC Berkeley and has been a faculty member in the Physics & Astronomy Department at the University of San Francisco since 2012. He works on problems in observational cosmology with collaborators in the Supernova Cosmology Project, the Nearby Supernova Factory, and the Dark Energy Spectroscopic Instrument experiment, and of course, with students.

Andi Gu is a current senior at UC Berkeley. He has been working in the Supernova Cosmology Project and DESI since 2019, applying his computer science and physics background to gravitational lens detection and modeling.

Host of Seminar:
Steven Farrell
Data & Analytics Services Group
National Energy Research Scientific Computing Center (NERSC)
Lawrence Berkeley National Laboratory

A

Welcome everybody. Thank you for joining today for the nurse data seminar.

A

We are uh very pleased and excited to have xiao xiang, huang and andy gu, talking to us today about their work on discovering and modeling strong gravitational lenses, not only because they're using uh nurse resources, but that's a cool extra extra feature of this work um and they're doing cool stuff using gigalands, which is, uh I guess, what I would call it: a sort of differentiable simulator and uh some machine learning techniques. So xiaoshan got his phd from uc berkeley he's now.

A

A faculty member at university of san francisco works on supernova cosmology project nearby, supernova factory and uh desi and andy is a senior now at uc berkeley and has also been working on the supernova cosmology project project in desi.

A

So they will both be giving the talk, and at this point I think I'll hand it off to xiaoshan. You can get started.

B

Thanks steve for inviting us, I'm excited to talk about strong gravitational lensing, and so I I think, given that the audience may not all have physics background I'll try to go through what lensing is at the beginning, but I also want to make sure that I leave enough time for andy so towards the end.

B

I may leave out a few things, but you know hopefully there there will be a q, a time so feel free to ask more about those aspects later and as steve and andy, and I were talking about earlier, feel free to interrupt any time all right. So let me share my screen.

B

Okay, hopefully you can see my screen.

B

And so today I will be talking to you about discovering and modeling uh strong gravitational lenses with corey and pro motor uh and nursk, and I will be talking more about the search aspect uh and in will talk more about uh the lens modeling side of things all right. um So, first of all, let me acknowledge my collaborators.

B

uh Davis lego was really the person that got this project started and has been a person. That's been sustaining this project going forward and there are quite a few physics division, people uh as well as people from other institutions, um but I want to give a special thanks to steve farrell, uh laurie, um steffi and roland thomas uh four years ago, when we started this project. There were so many things that we didn't know about.

B

Questions about nurse facilities, questions about machine learning and those three really were extremely helpful in getting the project off the ground. So thanks- and this is a project that has a number of undergraduate students. In fact, this is almost exclusively an undergraduate student research project.

B

I've been amazed astounded by the energy, the imagination, the expertise that they've offered and contributed to this project and finally, on the lens modeling side, we've added a few more collaborators in recent years, all right! So um what is gravitational lensing? Well, uh it is a chance alignment. It's due to a chance alignment between the observer, a foreground, massive galaxy and a faraway background galaxy, and this happens for about one in ten thousand massive galaxies and the mass of the foreground.

B

Galaxy uh will warp space time in this vicinity and this will cause photons the the path of photons to bend, and this will create highly distorted images and sometimes multiple images of the background object. So this is a simulation that is generated by one of the team members with them shu. So these two images correspond to the same object in the background and.

B

Here are a few examples of real lensing systems, so this is a doubly lens system, just like the simulation, and if the alignment is really good and if the lens mass distribution is reasonably spherically symmetric, you will get what's called an einstein ring and what that means is that photons can reach you from the background galaxy in all directions around the lens.

B

Now, if the lens has uh it has some ellipticity. Even if the alignment is perfect, you will not get a ring. Instead, you will get four images of the same background object and this this is called the einstein cross and you will see all three image configurations in the discoveries that I will show you later.

B

Lensing strong lensing, it turns out, is a powerful tool for cosmology.

B

It is the most direct test for co-dark matter model beyond the local universe, and so, if you want to know how that can actually how lensing can actually be used as a test, how you can ask me about that later?

B

If the background galaxy has a transient, meaning something that has a varying brightness so, for example, a quasar or supernova, then the the two copies of the same event will be offset in time, and we call that offset time delay. If you can measure that time delay and you can accurately model the mass distribution of the lensing galaxy, you will be able to measure directly the expansion rate of the universe or this quantity that astrophysicists call h, naught or the hubble constant and finally, in very rare, um there are these very rare systems.

B

They have two sources, so these images correspond to one source. The blue images correspond to another source. In those scenarios, these systems can be used to constrain dark energy, as pointed out by a recent paper by eric linder and one of his students.

B

um Okay, so then, how do we find these strong lenses? Well, you need a very wide and deep image survey, because these systems are so rare and the desi legacy. Imaging surveys are perfect for this purpose.

B

Desi is a dark energy spectroscopic instrument, it's a lbl, led doe, dark energy experiment, and the footprint covers much of the northern sky, as you can see here and terabytes upon, terabytes of data are stored and processed at nursk, and the objects are extracted and typed by this algorithm called tractor, which also runs on ersk, and let me point out that for such a large survey, if you want to finish that survey in a reasonable amount of time, it turns out that you really need more than one telescope.

B

So so so these are the three telescopes that did a survey and in part because of that, the survey is more homogeneous than most surveys than most imaging surveys. So, for example, the northern part which is shown with this green outline, is observed by these two telescopes, whereas the southern part of the survey, which is observed by the blanco four meter telescope with the red contour.

B

So the southern part has slightly better image quality, and this is something that that we, we need to think hard about how, when we try to look for lens lenses, and also you can see that the depth, the observation depth is not uniform, and this is also something that we had to think about how to overcome.

B

And so all this processing and nurse really provided a really good starting point for our search.

B

We first select galaxies that are brighter than 20th magnitude, so for astronomers, the smaller the magnitude value the brighter they are so we're looking for massive galaxies, and these galaxies tend to be among the brightest, and so we want to find these bright galaxies, which are likely to be massive galaxies.

B

uh Don't worry about these acronyms! These are basically uh uh indicating that these are elliptical galaxies, so these are the most massive galaxies in the universe and uh um most uh the vast majority of strong lenses are elliptical galaxies. So we focus on this type of galaxies. First rex are galaxies or objects that the tractor recognizes as galaxies, but couldn't tell what type it is so there's an unknown percentage of these galaxies that are elliptical galaxies.

B

So in total we wanted to search through 22 million galaxies to find strong lenses.

B

um So I think to this audience uh it's probably not a surprise that we wanted to use a resnet to do image recognition. um The architecture that we followed is from the news at all 2018.

B

We translate it into tensorflow. We do the usual image pre-processing and data augmentation.

B

We use a highly asymmetric uh training sample, meaning we use about 30 times more non-lenses than lenses in the training sample, and the reason for that is because we use observed images to train on the internet and when we started there were only several hundred known lenses.

B

That's one reason, so we are limited by the number of lenses at the time, but the other reason is that there are many different kinds of non-lensing configurations and we want the resnet to be exposed to as many of those configurations as possible.

B

um So, as I said, uh I can't remember if I, if I give you a numerical estimate, so we have about 30 times more non-lenses than lenses in our training sample.

B

Earlier I mentioned that the observational depth can vary. So what we found to be really helpful is to maintain the lens, the non-lens ratio, regardless of the depth uh that would that was very helpful in preventing the neural network to be biased one way or the other.

B

um We use the cross-entropy as the uh loss function.

B

Where y I is the label and y, I hat is the predicted probability and we achieve a uc of 0.992, and um so then we deployed on nursk, and um so we actually did two searches. The first was conducted on data release, seven of the legacy surveys and we found 335 lens candidates and and dr eight, when dr8 was completed, we did another search and we found four times more new candidates.

B

Here are the candidates that are elliptical galaxies, so we didn't include those lenses that are classified as wrecks, but those, but only those that are classified by the tractor as elliptical galaxies, um and so I mentioned earlier that image quality is a little bit better in the southern part of the survey and all of our training images are taken from this part of the survey, both lenses and the non-lenses.

B

So uh our performance, in this part, uh is, is, is really a sort of a fair metric to see how um the new net the neural network is doing.

B

um So, unlike for terrestrial um image recognitions in image recognition tasks to find lenses, we cannot rely on the resnet entirely. Human inspection is still needed.

B

It turns out that for our algorithm about 30, one in 30 of the images that are above the neural network threshold is classified by by a human as a lens candidate, so one in 30 is about as good as anyone can reach these days. This may be, uh you know one of the best, if not the best efficiency, since we did so well.

B

We thought why not apply our neural network trainer network to the northern part of the survey, so we did this without retraining the neural network and the efficiency drops a little bit, but not too badly. We still get one in 45.

B

um and finally, one of the team members- andy gool, from whom you will hear more later in this talk, modified the original resnet that we used a little bit by adding what he calls shielding layers, which are basically one by one convolution layers that compresses uh the output, uh the the number of output channels- and this reduces the number of training parameters- reduces training. Time and uh surprisingly, uh this modification mild, as it is added one-third more lens candidates, so neural network uh model.

B

Comparison, uh is something that we we thought about doing, uh but we haven't kind of done this in a systematic way. But this is just a an interesting uh early result that we found um so in total. We've discovered over 1600 lens candidates and those that are found in the second paper are broken down uh by the two models that we used by grades by the type. So dc stands for dead or comp which, as I said, are basically elliptical galaxies and by regions.

B

The decals is the southern part, and the mdos is the northern part.

A

I think we have a question from shane yep.

C

Go ahead, yeah just I wanted to make sure I was following this like 1 in 30 and 1 in 45. So are these false positives or these cases where it's leaving it's leaving out a candidate that should have been included.

B

Yeah, that's a really good question. uh No doubt some of these are false positives.

B

uh We found, for example, uh there there's a type of astronomical objects called ring galaxies, and uh uh um you know we kind of gradually got a little bit better in uh in terms of our ability to tell when um you know you know whether an object is a ring galaxy versus you know a lens einstein ring, but in other cases it's the evidence is just not enough, so we basically want to maintain you know a certain level of sort of candidate quality, so there's systems that we have a suspicion for that.

B

This may be lens, but just using ground-based images there just isn't enough um evidence for us to be confident and I'll show you examples of of of that kind of uh uh kind of uh hard cases later on. Sorry,.

C

Even a human has a hard time in those cases.

B

Yes, yes, so this yeah! This is because ground-based observations have inherent limitations.

C

Yeah and then the other question I had was you know if it's um leaving candidates out that should have been, is there any possibility that this could introduce any kind of biases that you might you know, draw conclusions from.

B

Yeah yeah, so uh if, if we want to do a population level study, then I would say that at this point the selection function is not absolutely clear.

C

B

um That's a great question and there is a way to address that and we're working towards that.

C

Okay, thank you.

B

B

um But I also want to say that lenses are so rare that each individual system is is worth studying in itself. You know, even if at the population level, we don't understand, what's what gets left out, what gets left out or what? What is kept in.

B

All right, um so, let's look at some of the discoveries, uh so these are from our first. First of all, all of our candidates are on our website and uh the vast majority are galaxy scale lenses, and you will see what I mean by that in just a minute.

B

It's not a surprise, because our training sample mostly consists of galaxy scale, lenses about 10, also our group and cluster scale lenses, and then about a dozen or so on lens quasars, which we deliberately excluded from from our training samples. So it's good to see that we can find these systems all right. So this is. This is from our first paper and the first column shows two near einstein rings.

B

As you can see, and this one shows an einstein cross, those four images and this one uh I would still call this one, a quad, because you got one two and then these this arc is the merging of two images.

B

uh This is a classic image, counter image configuration. So is this one an image and a counter image? What's interesting about the last two? Is that, as you can see, the lens really consists of two galaxies. As I said, we didn't have too many of those in our training example, so it was good to see that we can find those and don't worry about these numbers. We'll come back to this towards the end of my part of the talk.

B

This is the uh some of the systems that we found from our second search. uh Here are these big arcs and this dark? You can see a counter arc here, and this is a near einstein ring. uh These are doubly lens systems, I'm pointing at the two lensed images.

B

And then these are the quads. um This is faint, but if you look at it, there really are four images. There are four images and there are four images and then this is again a very a near perfect einstein, cross and I'll quickly go through these cluster lenses. These are very interesting, huge arcs, uh much larger, a factor of 10 larger than the galaxy scale arcs. So when I say group or cluster scale lenses, what I mean is that the lens consists of multiple galaxies in this case.

B

It's a group in this case it's a cluster, so cluster is basically a very large group and you can see these arcs there's a very thin arc here. There's an arc here, and this is one of our best systems. It's another group lending scenario, but it has four distinct sources: the labeled one.

B

So one has four images and two has this long arc kind of a broken in the middle uh three in the middle here and then four, this giant red arc out here um and so remember, multiple source uh lensing systems uh can be used to constrain, dark energy.

B

A

So sorry interrupt. Can I ask: how do you disentangle those things to know like ah these three go together? Those two definitely good use like simulations.

B

Very good question, so one uh one quick way is color, so gravitational lensing, unlike optical lens, is achromatic, so it does not change color. So, for example, these four images just visually. You can tell that they have very similar color, whereas this one's a little reddish this one's, clearly a different color and also geometry. So you can see that the aspect ratio of this and the color are more consistent compared with these four images.

B

um But I tell you that you know when a new student join uh the group, that's one of the first things that we do, which is uh kind of train them to recognize lenses.

B

um Yeah great question thanks, um so we thought. Oh, you know. We found a number of very interesting systems. We proposed to observe these systems on the hubble space telescope on the hubble space telescope and uh our proposal was accepted. uh We uh requested an infrared band to complement the optical observation from the ground.

B

uh In total, 51 systems were observed and every single one has been confirmed as an actual lensing system, and so now let me show you so these are half of the system 26 of the 51, and these are these are systems that primarily have a single galaxy as the main lens, as you can see, and some of them just by looking at it's actually very hard to tell so, for example, this one you look at it, you go.

B

Is that really a lens and we we thought we were confident enough to send this to hubble and you will see what hubble has what has hubble found. So let me now go to the hubble images. So here you can see that these arcs are really very clear.

B

Gravitationally lens, you can see the einstein cross here: another quadruply line system, another quadruply lens system and so forth. So this one was the one that you know when may legitimately ask the question: is that really lensing? It turns out that we were right, and so let me just for the fun.

A

B

Me kind of blink a little bit, um and so so the einstein radius is something that andy will talk more about. So it's basically the average distance between the arcs. Let me use that hse image. The average is the average distance between the arc and the lens um and the other half of the systems are group or cluster lenses. So here you clearly see there's a group of galaxies this one's a very interesting one again you can probably legitimately ask you know: is that really an arc?

B

Well, let's take a look at what it is. It is really arc. You can see actually other arcs, it's actually a very interesting, complex system. So is this one you can see there's basically little arcs everywhere.

B

I really like this one. If you look carefully it's a spiral galaxy with a bulge and with spiral, arms- and you can see it because of highly magnified by gravitational lensing.

B

You can see the spiral arm sticking out so again, you know I can. I can look at these for a long time, but in interesting time. Let me just blink back and forth. You can see that all of these are so you know this is almost complete einstein ring.

B

uh This is almost a full ring more like an egg shaped ring all.

C

B

uh Moving on um so we did it again because desi did it again, so desi kept applying their powerful image processing pipeline using nurse resources.

B

This time they went beyond the desi spectroscopic footprint to deck more negative than -18, so we follow suit and we did a search. So this is actually an existing survey called the dark energy survey and what's interesting about this search, is that this part of sky has been mined by other research groups for strong lenses already, but when we did it, we found thousands more galaxies, thousands more strong lenses. So, as I said, I I think that uh um you know machine learning, technique.

B

Comparisons could be interesting um and so basically, we've doubled the number of lenses compared to our first two papers and one of the team members christopher, is writing up. The results here are three of the new lenses that we found. This is a clear double. This is probably lens quasar. You can barely see uh the lensing galaxy, the red little few pixels down the middle crazy dots can be very bright. This is a giant arc and a counter arc produced by a group in the interest of time.

B

Let me just uh quickly walk you through some of the new candidates, so these are doubles. As you can see the quads. It's not always obvious. Here you see four red images here. You see four images.

B

You can ask me more if you want so here are clusters. This is arc, that's an arc, and this arc is interesting because in the middle, it's lensed again by this galaxy. So now it's got four images and then finally, the last row are all high red shift so very far away lens. Not only are the sources far away, but the lenses themselves are far away and you can tell by the deep red color of these galaxies in the middle and each system is really nice, so this is a quadrupole lens system.

B

This is another quadrupoli lens system. This one has a giant arc and a counter arc. So you you get the picture. um I I also want to say that uh we continue to work with uh desi uh dave.

B

If desi, if you don't know, is this powerful instrument. This amazing instrument that can measure the spectra of 5 000 galaxies with one pointing using fiber optics and so they've generously offered to observe our lensing systems. So why do we want to get spectra? Well, spectral, for one thing, spectra tells you spectra, tell you what's called redshift, so it's basically a measure of distance.

B

So one way to tell whether a suspected lensing system is actually a lensing system is that the lens is closer and what you suspect to be sourced is further away and typically much further away.

B

So the typical ratio of the red, the redshift for the source versus the redshift of the deflector cd or the lens, is a factor of two, and so let me just show you one example. This is very interesting, so this system already is confirmed by our hubble data, but desi has gotten the redshift for the lens. Well, we suspect to be lens, it's it's at a relatively low 0.25 and then we also targeted these two arc. Looking objects.

B

Not only are they much further away than this object, but they are the same distance, and not only that, but this emission feature are identical for these two objects, so in combination the red shifts of this guy being lower than these two and these two red shifts being exactly the same. They have the same spectra and they have a high resolution image that shows you the features that correspond to each other, for example this knot and that knot this little feature this little feature.

B

This is the most solid evidence you can get to confirm a lensing system, so it's really great and all the spectroscopic reduction is done on nurse and even the visual inspection tool called prospect is run on uh nurse jupiter hub and uh so, and does he even offer to observe system that has not yet been observed by hubble, and so these are subsets of systems that being confirmed spectroscopically as strong lenses. You see a nice arc here. You see almost einstein radius einstein ring here.

B

So with all these data coming in, we would want a reliable and scalable and fast modeling tool to extract science to reach our science goals, and so here I will hand it over to andy about lens model.

D

Perfect um all right session can I steal the screen. Okay,.

D

Perfect, can you guys see the presentation yeah.

A

D

Awesome: okay, um thank you, shaoxing and thank you everyone for having us, I'm really honored to have been invited to give this talk. So in my segment of the talk, I'm going to be talking about our new um fast differentiable bayesian lens modeling framework, um which we've called giga lens.

D

Okay. So before I start talking about our tool, um I think it's important to kind of walk through some very basic modeling theory.

D

So the primary lensing equation, or the central lensing equation, is expressing actually a very simple geometric fact that I've illustrated in the following diagram. um So here the sorry, let me see if I can get my cursor yeah. um So here the s marks the source and the disk marks a lens and o marks.

D

The observer and theta is the angle at which the observer um sees the image, and alpha is kind of this quantity known as the deflection angle, because it essentially measures how much the light from the source is bent by the lens um so for different types of galaxies. The functional form of alpha can actually differ depending on the mass distribution of the lens. um It turns out that dark matter consists of around eighty percent of the mass in a galaxy.

D

So that's why, as joshing mentioned before, modern lenses can help us understand uh dark matter better, probably the simplest archetypal lens model um is known as the singular isothermal sphere, where you basically assume that the stars inside the galaxy can be described as free atoms in a gas of uniform temperature so for systems belonging to this class of lens. The deflection angle is actually analytically expressible in a very simple form. It's actually just a constant, and this constant is known as the einstein radius.

D

So it's a parameter that depends on both the mass and the redshift of the lens, and you know, as you should mention before the reason it's called a radius is that if a source is perfectly behind the lens, um the image that is formed will have a radius exactly equal to the einstein radius, so even for more complicated lens models, there's still typically an einstein radius parameter for which the intuition is roughly the same. It basically measures the typical radius of the images produced by the ones and I've shown on this slide.

D

Here an example of a simulated lens and marked with the arrow um roughly what the einstein radius of that one corresponds to.

D

So our problem here is actually fairly typical. We do not have access to the ground truth parameters of the lens. um We only have the observable, which is typically just a noisy image of the lensing system, and to move towards the science gold that we mentioned before so, for example, dark matter testing a hubble, constant inference.

D

We need the ability to do several things um so using the observed image, we mainly want to learn about the light profile of the source, galaxy, uh the distribution of matter in the lensing galaxy and potentially even dark matter, halos sitting near the lens or along the line of sight, so line of sight, meaning um it sits between the observer and the source.

D

um So in the lamina in the language of bayesian statistics, what we have here is essentially an inference problem. Then, of course, as you might imagine, we will be making extensive use of the bayes equation, which essentially expresses the posterior distribution as a product of the likelihood and prior. um I just want to note uh you might catch me calling this inference problem, um modeling problem, but note that they refer to the same thing.

D

So I'm going to describe here basically a high level strategy for how we're going to tackle this inference problem. um The first thing that we need to do is to write down a model for the physics. um This includes a functional form for the source and the lens light, as well as the lens mass. So this mass model is always translated into a deflection angle, which is what we actually use in practice during the modeling.

D

um So the next step is to write down a model for the probability um that really means we need to define a prior for the parameters. In the likelihood. Oh sorry, we need to define a prior for the parameters as well as a likelihood function.

D

um The likelihood typically assumes that the noise at each pixel in the image is gaussian, so the log like, if it really looks like a squared error, or essentially just the chi squared up to a multiplicative factor, speaking heuristically, um all we're trying to do- is to match the simulated image to the observed image and the matching is again measured by the squared error. All the bayesian formalism does is to make this rigorous and to characterize an uncertainty for the parameters.

D

um So once we've written down this physical and probabilistic model we're basically ready to sample from the posterior um simple sounding enough.

D

um But when we actually start thinking about how to sample we're. Faced with a number of problems, the central one is that most regions of a high dimensional parameter space um has vanishing likelihood. So if we use some sort of out of the box sampler um out of the box, monte carlo sampler, like mcmc, we're likely to waste a lot of time in regions of the parameter space that we should not really be looking at because they have vanishing liquid.

D

So our solution, which is an approach which was an approach first introduced to us by an open source lens package, known as one astronomy, um is we're essentially going to look for the global maximum of the posterior density, and if we can find that, then and the posterior is unimodal in the sense that all the other modes of the posterior have vanishing likelihood.

D

We can just initialize our sampler at this point of high posterior density, and it will do much better well. The problem with this basic procedure is that it leads to another problem, which is that the parameter space is highly non-convex and therefore is also hard to search over. um So the the figure I've included on this slide shows a simple problem, which is that when you have a non-convex loss landscape, which in our case the lost landscape, is the posterior density.

D

Depending on your starting point, you may end up with different solutions um that different quality.

D

So with our goals and the problems from the last slide in mind, we had three main digital autumn uh when making a new lensing package. So first we required fast simulation. So as we head to that simulation is a key to forward modeling. um So if we want a fast inference code, we need a fast simulation code.

D

Luckily, most popular lens models are actually analytic. So more specifically, their deflection angles can be calculated through a series of very simple matrix operations. um So this is what I mean when I say the lens models are linear algebraic.

D

Since this is the case, we can actually use gpus for fast simulation, because we know that gpus are exceptionally good at elementary operations on large matrices.

D

So the next thing we need to do is we need some way to guide us through high dimensional parameter spaces. um We'll make the case that without some guide, so these methods are known as zeroth order methods. Modeling can be very inefficient, so, for instance, linstrom uses zero with order methods, and we found this to be an issue um for the primary reason that there's virtually no performance guarantees with these methods. um So for us we make our guide a local one, specifically the gradient.

D

um This is particularly fitting because, as I mentioned previously, lens models are analytic, hence differentiable. So therefore we can easily access the gradient via automatic differentiation, which is now technology, that's widely available in many many algebra libraries.

D

So, finally, we want a robust framework. um What this means is that it produces consistent results across different runs, so it does not fail often, and it doesn't require careful tuning of the initialization. So our approach to this is simple. We basically hedge our bets, that is as a general design principle.

D

We use many many candidate solutions so that the failure probability vanishes geometrically with the number of candidates we use so I'll elaborate further on this in coming slides, but to maintain numerous candidate solutions, we'll use parallel computation with gpus.

D

So we developed our framework known as gigalens, which is an abbreviation for a gradient, informed and gpu accelerated lens modeling framework. We wrote this with the three designer atom from the last slide in mind um and, as some of you may have realized, the popular library tensorflow is a good choice for most of the features that we wanted.

D

So our code is fully vectorized, fully differentiable um using tensorflow's, auto diff technology, and it also runs natively on gpus, but furthermore, we've also written a parallel branch using jaxx that has further support for distributed computing using multiple gpus. So we should note that we also make um extensive use of the tensorflow probability library as kind of a substrate for many of the bayesian inference tasks. We do in our modeling procedure.

D

Okay, so now that we've established or written down a simulation code, it basically remains to actually apply it and start modeling.

D

So, as I mentioned before, our strategy is to first identify the global maximum of the posterior density, and then we use this to inform how we initialize our sampler. um So we do this with multi-starts gradient descent, and what this means is that we initialize many candidate solutions by sampling from our prior, which ought to be broad enough to include the global maximum.

D

Next, we use a bijector to unconstrain the various parameters in our learn's model um so, for instance, uh take the einstein radius. So clearly it must always be positive, because it's a length, um however, it's undesirable to like manually, put this constraint into the gradient descent with some sort of like hard wall. um So therefore, we use bijectors that map parameters into an unconstrained space um so for enzyme radius.

D

We take the log and optimize over that um for a different parameter, known as the mass slope, um which basically measures how fast the mass distribution is changing uh at einstein, radius, roughly speaking um so from a physical argument, we know that it always sits between one and three um and since we know a sigmoid maps of the real line to a bounded interval, we can use the inverse sigmoid to unconstrain the slope gamma.

D

um So the next thing we do is locally optimize each candidate solution using gradient descent with respect to the posterior density. So I want to emphasize we're doing this optimization in parallel over every single candidate solution and we typically use atom with an annealed learning rate, as is typical for non-convex optimization problems.

D

um So, finally, after a fixed number of iterations, we stop and select the best performing solution and use this as the maps.

D

So how well does this perform in practice? um Let's look at an example system uh on the left. I've shown the trajectory of each parameter um or of two primary parameters: uh one being einstein, radius and one being the mass slope. um So specifically, I'm showing the error with respect to the ground truth of each parameter as a function of the gradient descent iteration.

D

So the solutions highlighted converge, the correct solutions, uh so the ones in bold and we've marked in red boxes. Basically, the error for the two key parameters um so note that the failure probability of each particle is here roughly 96 percent. uh In the test we actually did for our paper. um We found uh 98 upper bound on the failure probability, but for 300 particles that actually results in success, so success, meaning identifying the global maximum. um It results in a success being virtually certain um on the right.

D

We show the lost trajectory for the best performing six particles or the six uh best performing six solutions. uh We see that by the 300th iteration the particles seem to have come into the neighborhood of the optimum, because after this they exhibit a geometric convergence which is characteristic of loss, trajectories in convex regions of the loss function. um So that's. This actually informs our cutoff iteration for the gradient descent.

D

Okay, um so with the map estimate uh before sampling, we still have one more important auxiliary step. um We want to estimate the covariance, which will actually be important later, but the way that we do this is we estimate it using variational inference. So this uh procedure, variational reference- is essentially one where you fit an onsite's distribution, which we call our surrogate posterior to the true posterior.

D

So, as indicated on the diagram um say, our posterior is kind of like a funny object like the one in blue, for which it's hard to get an analytic expression.

D

If we want a rough characterization of it, we might use a simpler analytic distribution and try to fit this to the ground truth or the true posterior and the fit here is measured by the kl divergence between the two distributions. So for us, since we want the covariance we'll make the surrogate a multivariate normal and you know later I'll touch on whether or not this is a reasonable approximation.

D

um But with this we can run vi by optimizing, the kl as a function of the mean and the covariance, so the kl up to an additive constant, can be written as something known as the evidence lower bound called the elbow. This is expressed as an expectation over the surrogate posterior of a difference between the log densities, so the gradient is actually similarly expressible as an expectation with an additional gradient term just going into that expectation.

D

So to optimize the kl I divergent or our strategy will basically be as follows: we'll initialize mu at the map estimate from the previous step and then we'll set an initial small diagonal covariance, with a variance of 10 to the minus six for each parameter um and then we'll minimize the elbow by approximating the two expressions um I've shown on the slide using monte carlo estimates. So that just means we sampled many thetas from the surrogate posterior calculate gradients average and then update using gradient descent.

D

So the final optimized covariance matrix is an estimate for the true covariance um so for a frequentist interpretation of what's going on here. This can also be understood as taking the hessian of the log posterior at the mode, um so indeed this actually works pretty well for just as well as vi for small parameter spaces, but actually we found in high dimensions the hessian becomes kind of unstable and that vi is significantly more reliable. So we opt to use vi rather than taking the hashtag.

D

So, finally, we're prepared to sample, um since we have a gradient since we have gradient information uh natively using giga lens, we can use a sampling technique known as hamiltonian monte, carlo, um so this sampler was proposed to handle the inefficiency of mcmc in high dimensions. um This inefficiency essentially has its origins in the idea of a typical set. So when we're sampling, we want to spend most of our time in an area known as a typical set, which can be very non-trivial in high dimensions.

D

So it can exhibit some sort of like ring structure, as shown on the right. um So while mcmc will make a lot of rejected proposals with the structure like this and hence be inefficient, hmc is actually able to use gradient information to follow this typical set very closely and we get very efficient sampling. um Hmc actually requires a couple of hyper parameters and uh initialization configurations.

D

um The first thing is that we initialized hmc walkers, with samples drawn from the surrogate vi posterior from the previous line, and the next point is kind of critical, which is that we set the mass matrix to be the inverse of the covariance matrix.

D

So this is essentially to inform hmc of the various scales associated with each parameter, which can vary significantly for our problem. So some parameters might have a variance of 10 to the minus 3, and some have just 0.1 and hmc does not handle this anisotropy well without um a mass matrix or like a set mass matrix to basically inform.

A

D

You have these different scales, um other hyper parameters like the step size and the leapfrog length. Those are all adaptively tuned over um the course of the burnin.

D

So to give you an idea of the performance improvements, we're able to extract using hmc, we compare it to a popular ensemble sampler in astrophysics, known as emc. um So look at the so take a look at the trace plots for the einstein radius alone. We can see almost purely independent samples being drawn at every iteration with hmc, whereas there's a high degree of correlation uh in emc. So this is quantified using metrics like the effective sample size which measures the number of independent samples drawn in a sampler.

D

So it's also measured by something called the potential scale reduction factor, r hat, which basically measures what degree to what degree different rockers uh are sampling from different distributions. So, ideally we want them to look like the top marginal distribution where all the walkers show very similar distributions, um and this leads to an r hat, that's as close to one as as possible. um So we see something far less desirable for emc, where we get in our head of roughly 1.5 um and 1.1 is kind of the standard that most people typically reach.

D

For so you can also see a similar degree of performance differences in the effective sample size. There's almost a two two orders of magnitude: difference in the effective sample size.

D

So with that that actually completes our modeling pipeline uh once hmc is done, sampling we've drawn many samples from the posterior which we can now do inference with. um So here. I've shown modeling results for one typical system, which is close to an einstein ring um the quarter. Plot shows the distribution of hmc samples, as well as the bi posterior. So we see that vi actually did a fairly good job of capturing.

D

The typical scales of the posterior observed, like the scale of einstein radius, is 10 to the minus 3 and gamma is 0.07, so we have very different scales here and vi was actually um actually able to capture that. um However, vi can't be perfect because we found our posterior is not perfectly gaussian, despite how it may appear in the corner plot. So, for instance, we show a cross section of the posterior found by slicing along a two-dimensional subspace on the top, where the title says a non-gaussian parameter cross-section.

D

So we see that our posterior is not perfectly gaussian, because there's um the the true posterior is not perfectly gauging along this cross section. um So on the right, we've shown a diagnostic for the time breakdown of our modeling. um So essentially note that um vi is actually um the most expensive part of our modeling when we run on four 800 a100 gpus uh on pearl motor early access, um so I'll actually address this later on um and some ideas for what we have to address this.

D

So now, I'm going to show modeling results for four archetypal systems that were uh simulated um so observe the significant disagreements between vi and hmmc. Here um power. It actually still captures the scale roughly accurately, which is all we need. um So, despite um despite the uh in performance of the the er, despite the powerful performance of the vi, um the reconstructions are still pretty good and we've basically found an unbiased fit to the observed image.

D

So here are two more systems. um The one on the left is known as a cross, so josh, I showed you um some real observed images uh that exhibit this cross structure. This is a simulated. um The one on the right is known as a double um so notably. Actually, the vi is almost perfect for the double um so to verify to verify the correctness of our modeling code.

D

We simulated a hundred systems blinded and then modeled them and comparing our inferred parameters with the ground truth, we found no bias as well as well calibrated uncertainties, so the mean chi square for each of these is very close to one. um So this was actually um a good test for us to verify the speed, the speed and the robustness of our code.

D

So, finally, I'm actually excited to show modeling results on a real observed systems. um What we have here is a classic quad, so you can see um in the false color I've indicated the quad configuration in gold and then there's also a double which I've shown in the silver.

D

So we actually have here two sources being lens making it a very interesting system, um and we found that giga lens was actually able to successfully model the system. So take a look at the residual um you'll find that it looks pretty good.

D

Okay, awesome so to wrap up um our paper reporting. This work um has been posted on archive and has also been accepted for publication in the astrophysical journal. um What next? Well, the first thing we want to do is develop a suite of comprehensive, automated tests, as well as improve the documentation for our code, and this goes to our second objective, which is to do a public code release so that we can see um the interesting results that other people can get when using giga lens.

D

um The next thing we want to do is improve the speed of our modeling code by using two gpu nodes, um so that means eight gpus, uh and this is actually an option. We hope to use uh like the configurability option on promoter for um and a sub point of, that is to ask.

D

Can we eliminate the vi step and an idea that we have, for that is to use mass matrix adaptation for hmc, so uh we adaptively tune the mass matrix out um the mass matrix over the burning steps as well, and the final thing is to apply our code to observe systems from our hst program.

D

So there's two potential goals with this. uh One being sub galactic dark matter, halo detection and the other being using ground-based data. Can we model the systems. um Spectra spectroscopically, confirmed using desi yeah that about wraps up um what I have to talk about for modeling. um I'm sure me and joshua would be happy to take any questions.