likelihood functions for imaging

Mario Juric (UW) showed up for the day and we spoke for hours about many things. One category was image likelihood functions (for things like The Tractor or weak lensing). We came up with a very dumb (read: good!) idea for testing out ideas around likelihood functions: Take two image data sets from different telescopes that overlap on the sky. Build a catalog of sources (with positions and colors and so on) from a joint analysis of both data sets. Then do the same, but in a world in which your only interface to each data set is a callable API to a likelihood function! That is, something that takes as input a parameterized high-resolution image model and returns a likelihood value, given what it knows about its data and calibration, PSF, and so on. This would force us to figure out what would be needed in such an API. I think we would learn a lot, and it would help us think about how to construct next-next-generation data products. We also talked about image differencing, Dun Wang's Causal Pixel Model, and other matters of mutual interest.

code seminar

The research highlight for the day was a seminar led by Jeremy Magland (SCDA), about MountainView, his deployed code that sorts and presents data on neural recordings (cuts the data into spikes and clusters the spikes into neurons and presents the results visually). The great thing about the seminar is that it wasn't really about neuroscience, it was about code: How to structure and build and maintain a project of this scale, which works on both local and remote data. The audience was asked to weigh in on design issues and react to design choices. It was a very productive discussion and is a model for events we should do next year jointly between the SCDA and the new SCCA; both institutions will be building and supporting non-trivial software projects.

Gaussian processes, black-hole dark matter

At group meeting, Dan Cervone (NYU), who works on spatial statistics for sports and climate, went over the basics of Gaussian Processes for interpolation. This is all related to Boris Leistedt's project to determine galaxy photometric redshifts with a flexible spectral energy distribution model.

I had lunch with Kat Deck (Caltech) and Yacine Ali-Haimoud (JHU) and Kyle Cranmer (NYU). We talked about various things, but especially whether the dark matter could be massive (LIGO-detected!) black holes. One idea that would be easy to look at, that Scott Tremaine (IAS) mentioned to me a few years ago, is whether the dark matter granularity could be limited by looking at the dispersal of cold tidal streams of stars.

Late in the day, I worked on my few-photon image-reconstruction toy problem, tuning the stochastic gradient and going to smaller numbers of photons per image. It works well and I am excited about the implications for diffraction microscopy.

stochastic gradient FTW!

So I have been working on this toy problem: Imagine you have a scene, and each "image" you have of the scene contains only a few photons. Each image, furthermore, is taken at a different, unknown angle. That is, you don't know either the scene nor any image's angle, and each image contains only a few photons. This is a toy problem because (a) it is technically trivial, and (b) to my knowledge, no-one has this problem to solve! (Correct me if I am wrong!). I am using this problem as a technology development platform for the diffraction microscopy projects I want to do.

One piece of infrastructure I built today is a stochastic-gradient optimizer, that reads in just one image at a time (just a few photons at a time) and takes a gradient step based only on those new photons. This is standard practice these days in machine learning, but I was skeptical. Well it rocked. The image below (believe it or not) is my reconstruction of my scene. I also show four of the “images” that were used as data. The images are just scatterings of photons, of course. (As my loyal reader knows, the objective of the optimizer is a marginalized likelihood.)

I also had a great lunch with Kathryn Johnston, in which we discussed our past and future projects in Milky Way astrophysics and why we are doing them.

AAAC meeting, hoggsumexp

In big marginalized likelihood calculations, there are usually some "logsumexp" calculations floating around. Today I wrote a custom logsumexp that returns the answer (in a numerically stable way) and also derivatives with respect to parameters, for use in my diffraction microscopy (and later galaxy-pose) projects.

I spent most of my research day on the Astronomy and Astrophysics Advisory Committee meeting (that advises NSF, NASA, and DOE on areas of mutual overlap). The most interesting thing I learned on the call, from a research perspective, is that NASA now has very long-duration balloon technology that can fly around the globe many times, with no natural time limit! The balloons can carry 5000 lbs. This may create a lot of new low-cost opportunities for astronomy outside the visible.

toy problem

Friday got messed up, and only a tiny bit of research got done, but I made up for it on the weekend, building a toy version of the hard problem I presented to the mathematicians on Thursday. The toy problem is to reconstruct an image with few-photon examples, each taken at a different rotation angle. The form of the likelihood is similar, but the computational cost is much lower. It should be an intuition-building example.

talking stats at the mathematicians

I spent a very large part of the day today at the whiteboard in front of Charlie Epstein (Penn), Leslie Greengard, Jeremy Magland (SCDA), and Marina Spivak (SCDA). I presented my proposed solution to the problem of diffraction imaging of molecules in the limit of very few photons per exposure. We had a brief discussion of the physics, a very long discussion of the idea of solving this problem by optimizing a marginalized likelihood, and a brief discussion of its derivatives with respect to parameters in a representation. It was an incredibly useful session: The crowd found some mistakes on my part, and it forced me to clearly articulate how I think probabilistic inference works in these cases.

I think my proposed solution is palatable to both Bayesians and frequentists: In principle the frequentists should object to my marginalization over angles, but this is close to unassailable, because when the angles are generated by an isotropic process, they really do have a well-defined distribution that does not have to be described as a prior. That is, this integral is frequentist-safe. In principle the Bayesians should object to the fact that I am going to optimize the (marginalized) likelihood rather than fully sample the posterior, but even hard-core Bayesians recognize that sometimes you can't afford to do more than get the peak and its width!

Amusing notes from the chat: My definition of “principled” is very different from Leslie Greengard's! And the crowd was not as confident as I am that I can solve a problem where the intermediate steps to the answer require as much disk storage space as exists in all of facebook (note the reference to “10¹¹ to 10¹⁷ numbers” on the board).

Dr Ana Bonaca!

Today I had the great pleasure to sit on the PhD defense meeting of Ana Bonaca at Yale. She defended her thesis work on measuring the mass and gravitational acceleration field of the Milky Way using cold tidal streams. Her thesis is comprehensive: She develops a new method for performing the inference (a good likelihood function!); she compares the results obtained for (artificial) streams made in simple potentials and realistic, time-dependent potentials; she finds (in real data) a new, cold stream (incorrectly named Triangulum Stream) that is of great value for this work; and she performs the first ever measurement of the Milky Way using two streams simultaneously (GD-1 and Palomar 5). A great set of projects and a great seminar and discussion.

Bonaca, Marla Geha (Bonaca's advisor), and I spent lunch afterwards fighting with cosmologists Nikhil Padmanabhan and Frank van den Bosch about why this work is important, with the cosmologists taking the (exaggerated, I think) position that the Milky Way doesn't matter to cosmology! That was fun, too!

yet more derivatives

I spent the day working out derivatives of my diffraction imaging model with respect to parameters. The crazy thing is that there is a sum over components in the density model (for the molecule) and then a sum over photons within the image (or instance) and then a sum over orientation-angle samples, and finally a sum over images (or instances). And various of these sums have logs and exps inside them. So it is a mess! I wrote them out with a pen on paper, and then typed them up in a nascent paper. The short-term plan is to get to stochastic gradient.

testing works

OMG, testing works! (Duh.) I wrote proper tests for my celestial mechanics code, and sure enough I instantly located my bugs. Now I have working code, and tests to prove that it is working! Aaahh.

computational science, microlensing, and etc

In the morning I gave an informal talk at the Simons SCDA, about my work on The Cannon, and my issues (both good things and bad things) with machine learning. I discussed the point (which inspires our ABC research, and which I also discussed with Kravtsov at Chicago) that quantitative natural science is now almost entirely computational—meaning that the theory is a simulation that makes artificial data—and this leads to changes in how we do inference and speak about realities.

In the afternoon, I spoke with Megan Bedell (Chicago) about the echelle spectroscopy radial-velocity data she has; she has done some dimensionality reduction and there are promising opportunities I think to improve the end-to-end radial-velocity precision. I also worked on my celestial mechanics code for Price-Whelan; it is not working and I don't know why! My only option is to write proper tests. Tomorrow!

While all this was going on, in the background, Dun Wang has been steadily finding cool stuff in the K2C9 data, including, for example, this (previously known from the ground) baby!

new discoveries in the K2 data!

At group meeting, Dun Wang showed new discoveries in the brand-new K2 Campaign 9 data. The K2 team released the data the moment it came off the s/c, in fact before they even ran their own data-reformatting code. This was no problem for us, because K2 god Geert Barentsen (Ames) has (clean-roomed? and) open-sourced some of the core code. Wang ran our image-prediction / image-differencing version of our CPM code to predict the data and identify sources with strange excursions. He found ones that looked like possible microlensing events (known and unknown) and published them to the K2C9 group. I asked him to re-run CPM with other parameters to make it less (and more) aggressive and thereby address the (probable) over-fitting. The next step will be to incorporate a microlensing model and fit the K2 systematics (the point of the CPM code) and the microlensing parameters simultaneously. We discussed next steps.

Later in the day I started actually writing the celestial mechanics code my new eclipsing-binary team (Price-Whelan, Ness, and Foreman-Mackey). It is probably all wrong, but it is core technology, so it needs to get instrumented with tests.

celestial mechanics

As part of my project with Adrian Price-Whelan (and also Melissa Ness and Dan Foreman-Mackey), I spent my research time today figuring out time derivatives of Kepler's equations. These so we can do simultaneous fitting of the eclipsing binary light curve and the radial velocities revealed in the double-line spectrum. This was actual pen-on-paper calculus! It's been a while, although as I took these derivatives, it reminded me that I have taken them many times in the past.

In the afternoon I had a great conversation with Duane Lee (Vanderbilt) about chemical tagging and nucleosynthesis. He is close to being able to fit our data in the Milky Way halo with a mixture of dwarf-galaxy stellar populations. That would be awesome! We talked about low-hanging fruit with our APOGEE chemical abundance data.

a new project: eclipsing binaries

I have a dream! If we could get enough long-period eclipsing binaries with multi-epoch spectroscopy, we could go a long way towards building a truly data-driven model of stellar spectra. It would be truly data-driven, because we would use the gravitational model of the eclipsing binary to get stellar masses and radii, and thus give good label (truth) inputs to a model like The Cannon for the stellar spectra. (Yes, if you have an eclipsing binary and spectroscopy for radial velocities, you get everything.) And then we could get densities, masses, and radii of stars for the interpretation of transit and radial-velocity results on exoplanets, without relying on stellar models. There are lots of other things to do too, like build population models for binary stars, and exploit the stellar models for Milky Way science. And etc.

Today, because of a meeting cancellation, both Adrian Price-Whelan and I got the full day off from responsibilities, so we decided to use it very irresponsibly. We searched the (very incomplete and under-studied) Kepler eclipsing binary list for binaries with long periods, deep eclipse depths, and APOGEE spectroscopy. It turns out there are lots! We started with the system KIC 9246715, which is a red-giant pair.

In the APOGEE spectrum, the pair of velocities (double line) is clearly visible, and it clearly changes from epoch to epoch. We found the velocities at each epoch first by auto-correlation and then by modeling the spectrum as a sum of two single stars. A project is born!

extreme precision radial velocities

I continued working on my document about release of data and code. Twitter (tm) continues awesome.

Research highlight of the day was a long discussion with Megan Bedell (Chicago) about the consistency of exoplanet-host radial-velocity measurements order-by-order in a many-order high-resolution echelle spectrograph. The measurements are made by cross-correlation with a binary (ugh) template, and some orders are consistently high and some are consistently low, and we very much hope there are other more subtle regularities to exploit. Why are there these discrepancies? Probably because the model is inflexible and wrong. Unfortunately we don't have access to it directly (yet) so we have to live with the cross-correlation functions. We discussed simple methods to discover regularities in the order-by-order offsets and results, and sent Bedell off with a long to-do list.

I ended the day with a long conversation with Kat Deck (Caltech). Among other things, we discussed what we would do with our lives if exoplanet research evolves into nothing other than atmosphere transmission spectroscopy and modeling. Of course neither of us considers this outcome likely!