Terra Hunting Fall Science Meeting, day 3

Today was a delight! In a working session, Clark Baker (Cambridge) gave a beautiful, conceptual and concrete description of how an echelle spectrograph works and the blaze and the resolution and etc. My favorite moment was the aha! moment I had when he described the Littrow condition. This was followed by Alicia Anderson (Cambridge) explaining how the data reduction proceeds. Then she and Federica Rescigno (Exeter) helped us install the data-reduction software for the ESO instruments (ESPRESSO, HARPS-N, etc) and we started reducing raw echelle data.

Before all this there was a wide-ranging discussion of measuring 3-point functions of radial-velocity time series data. This was inpired by the question: Is a Gaussian process a good model for these data? I hope this turns into a project or set of projects.

Terra Hunting Fall Science Meeting, day 2

So many good things happened in the meeting today! Highlights were presentations by Niamh O'Sullivan (Oxford) Ben Lakeland (Exeter) who showed amazing results running models of stellar variability on data from the Sun. O'Sullivan can see that the sun goes through many different phases of spots, granulation, and super-granulation. She finds these by fitting Gaussian processes of certain forms. Related: Suzanne Aigrain (Oxford) showed that even in very gappy data, the GP fits are unbiased, whereas naive use of periodograms is biased!

Lakeland showed that super-granulation can in principle be modeled in the Solar time series, and maybe the tiniest hint that when he corrects for super-granulation well, the RV variability might be even lower than at times at which there is no super-granulation in play at all. Does super-granulation suppress other kinds of variability?

I'm very optimistic—between Liang yesterday, Zhou's work at Flatiron, and these presentations—that we will be able to mitigate many difficult sources of stellar variability. I was inspired to outline a conceptual paper on why or how this is all going to work.

writing an interpolator in jax

Today I pair-coded with Matt Daunt (NYU) a general kernel-based interpolator in jax. Actually, I am not certain that what we wrote can be handled by jax gracefully, because I don't understand the functional programming model. But I learned that interpolation is harder than I thought: If your interpolator is higher-order than linear, the interpolation involves solving an inverse problem. The fast interpolators (like spline) have such solutions coded with clever iteration schemes that require only one pass through the data, like the fast GPs we love. But if you want to write something general, you might have to suck up some (sparse) linear algebra. Daunt and I discussed strategies for our spectroscopic model, given these realities.

Hekker group visit

Today I visited the group of Saskia Hekker (HITS). We discussed many things asteroseismological! We discussed:

• the ESA Plato observing strategy
• is the asteroseismic signal a Gaussian process to any degree of accuracy?
• using asteroseismic information to improve and inform open-cluster membership
• synchronization of orbital periods with primary-star rotation periods
• are two distributions different?

and much, much more. I had a lovely day at HITS.

GPRV, day 1

The GPRV meeting started in Oxford today. The meeting brings together people working on data analysis in extreme precision radial-velocity projects, but united by interests in and uses of Gaussian processes. The first day ended with a very nice tutorial by Foreman-Mackey (Flatiron) on applied-math and computational tools for scalable Gaussian processes. He even live-coded and blew everyone's mind with Python jax.

Many talks (including Barragán (Oxford), Delisle (Geneva), and Tran (UT Austin) to name a few) are using Gaussian processes and their derivatives or two Gaussian processes to model the star's variability, with photometry, radial-velocity measurements, and activity indicators modeled as linear combinations of these latent processes. That's a really interesting theme, and connects somehow to my evil plan (with Bedell, Luger, Zhao, et al) of modeling the whole stellar surface. It is definitely an exciting time.

One issue that came up is how to judge or assess over-fitting. There was no consensus or answer, and most of the GP practitioners are very Bayesian. But Bayesian approaches aren't always sensitive to true statistical violations of the model; I want to see some cross-validation in this house.

In other news, Halverson (JPL) told us about publicly available solar data (and lots of it) from NASA NEID. I might want to play with that when I get home!

generative models for quasars

I spent part of the day working with Christina Eilers on her Gaussian process latent-variable model for quasar spectra and physical properties. We re-wrote our title and abstract and went through the math in the paper. It's time to finish this up! We find that we can predict quasar masses with good accuracy (based on held-out data) based on single-epoch, limited-coverage optical spectra. It's sweet. And Eilers has beautiful demonstrations that she can predict unobserved spectral regions, because the model is trained on different quasars at different redshifts with different data. The big problem with this model is that it scales poorly; we can't imagine training on thousands of objects without substantial engineering efforts (and maybe not ever).

finishing the response to referee and adjusting the paper.

As is usual with Publications of the Astronomical Society of the Pacific (great journal!), Soledad Villar and I got a constructive and useful referee report on our fitting paper. We finished our comments and adjustments to the paper today. The referee made an excellent point, which is: Since there are fast Gaussian process codes out there, why ever do interpolation or flexible fitting any other way? Good question! We answered it in the new revision (because sometimes fast GPs don't exist, and sometimes you don't want a stationary process and sometimes you are in a weird geometry or space), which we will update on arXiv soon.

more MPIA dust mapping

At Milky Way Group Meeting at MPIA, Thavisha Dharmawareda (MPIA) showed her first results on building pieces of a three-dimensional dust map from observed extinctions/attenuations to stars. As usual, the problem is to infer a three-dimensional map, preferably non-parametrically, from measurements of line-of-sight integrals through the map. She uses a Gaussian process, variational inference, and inducing points. She has some nice features in her maps (she started with star-formation regions with interesting morphologies). She sees extinction-law variations too; we discussed how those might be incorporated.

asteroseismic p-mode noise mitigation

I had a call with part of the HARPS3 team today, the sub-part working on observations of the Sun. Yes, Sun. That got us arguing about asteroseismic modes and me claiming that there are better approaches for ameliorating p-mode noise in extreme precision radial-velocity measurements than setting your exposure times carefully to null the modes. The crew asked me to get specific, so I had a call with Bedell (Flatiron) later in the day to work out what we need to assemble. The issues are about correlated noise: Asteroseismic noise is correlated; those correlations can be exploited for good, or ignored for bad. That's the argument I have to clearly make.

finishing a paper

I spent time today working through comments from Kate Storey-Fisher (NYU) on the document that Soledad Villar (JHU) and I have written about fitting flexible models. I made those changes, while Soledad put in some proofs of some of the key math points. We are so close to being done! But I don't mind being slowed down by amazingly constructive and useful comments from my students!

zeroth draft done!

Today I finished the zeroth draft of a (first-author; gasp!) paper about linear regression with large numbers of parameters. My co-author is Soledad Villar (JHU). The paper shows how—when you are fitting a flexible model like a polynomial or a Fourier series—you can have more parameters than data with no problem, and in fact you often do better in that regime, even in predictive accuracy for held-out data. It also shows that as the number of parameters goes to infinity, your linear regression becomes a Gaussian process if you choose your regularization correctly. It is designed to be like a textbook chapter so we are faced with the question: Where to publish (other than arXiv, which is a given).

writing down the linear regression GP relationship

My loyal reader knows that I have been working on the fundamentals of linear regression for a bit now. Today I did some writing on this topic. Last week, Soledad Villar (JHU) and I got the point that we could write down a specific case where the limit of infinite features in a particular, carefully designed linear regression becomes exactly a Gaussian Process with a particular, carefully chosen kernel. I understand how to generalize this result partially, but not completely: Apparently this will work in an infinity of different bases, with an infinity of different weighting functions or kernels. My goal is to write something pedagogical and useful for practitioners.

human aspects of data analysis

We had a fun data group meeting today, in which we discussed many human aspects of data analysis (like asking questions in talks and seminars, and sharing work when it is in pre-publication status). I spoke about the connection between Gaussian processes and linear fitting with enormous numbers of basis functions; there is a limit in which they become identical, which is awesome. Group meeting was followed by a conversation with Storey-Fisher about what we are going to work on next: Pulsar timing? Looking for anomalies in large-scale structure? Intensity mapping?

big, huge linear regressions

I spoke (remotely) at CCA today about linear regression (fitting linear models for the purposes of prediction), when the linear regressions have huge numbers of parameters. Yes huge: More than the number of data points! It turns out that even though you can thread the data perfectly—your chi-squared will be exactly zero—you can still make good predictions for held-out data. That surprised the crowd, which, in turn, surprised me: Many in this crowd use Gaussian processes and deep learning, both of which have these properties: More parameters than data, can fit any training data perfectly, and yet still make good, non-trivial predictions on held-out data.

My slides are here. Should I write something about all this?

forward modeling of asteroseismic modes—for exoplanets

As my loyal reader knows, I have four separate projects with Bonaca, Feeney, Casey, and Bedell to forward-model asteroseismic modes, with Bonaca concentrating on ground-based data, Feeney on principled Bayes, Casey on ESA Gaia, and Bedell on removing the modes as nuisances when we want to find planets. Today Bedell and I discussed where we are at, and came up with some things to try. In Sun-like stars, the modes have days-ish coherence times (apparently) and the modes have minutes-ish periods. So there are different regimes as your exposure cadence ranges from minutes to days to weeks. We have some qualitative predictions, and we are trying to make quantitative results that will influence survey design in the near future (especially for Terra Hunting Experiment and NASA NEID).

One funny thing about quasi-periodic (as in: finite coherence-time) oscillation processes is that they can be generated as a subset of Gaussian Processes, if you have good kernel machinery. We do! Another funny thing is that a GP can fit anything it is asked to fit! It literally has infinity free parameters (yes, literally). But the more appropriate kernels will do better (we hope) at predicting held-out data.

GPs for p-modes in RV data

I checked in with Megan Bedell (Flatiron) on our projects today. She showed really nice results in which she fits simulated radial-velocity data for a star that is oscillating in finite-coherence asteroseismic modes. My loyal reader knows that we have been working on this! But the cool thing is that she can mow fit the oscillations with a Gaussian Process with a kernel that is roughly correct, or exactly correct, even when the observations are integrated over finite exposure times. That's a breakthrough. It depends in large part on the magic of exoplanet.

Now GPs are extremely flexible, so the question is: How to validate results? After all, any GP can thread through any set of points. We came up with two schemes. The first is a N-fold cross-validation, in which we train the GP on all but 1/N-th of the data and then predict that 1/N-th, and cycle to get everything. First experiments along these lines seem to show that the more correct the kernel, the better we predict! The second is that we make fake data that includes the p-modes and a simulated planet companion. We show that our planet-companion inferences become more accurate as our kernel becomes more accurate.

We're hoping to improve on the results of this paper on p-mode mitigation. My conjecture is that when we use an accurate GP kernel, we will get exoplanet inferences at any exposure time that are better than one gets using the quasi-optimal exposure times and a “jitter” term to account for the residual p-mode noise.

interpolation duh

As my loyal reader knows, I am working with Lily Zhao (Yale) to calibrate the EXPRES spectrograph. Our approach is non-parametric: We can beat any polynomial calibration with an interpolation (we are using splines, but one could also use a Gaussian Process or any other method, I think). The funniest thing happened today, which surprised me, but shouldn't have! When Zhao plotted a histogram of the differences between our predicted line locations (from our interpolation) and the observed line locations (of held-out lines, held out from the interpolation), they were always redshifted! There was a systematic bias everywhere. We did all sorts of experiments but could find no bug. What gives? And then we had a realization which is pretty much Duh:

If you are doing linear interpolation (and we were at this point), and if your function is monotonically varying, and if your function's first derivative is also monotonically varying, the linear interpolator will always be biased to the same side! Hahaha. We switched to a cubic spline and everything went unbiased.

In detail, of course, interpolation will always be biased. After all, it does not represent your beliefs about how the data are generated, and it certainly does not represent the truth about how your data were generated. So it is always biased. It's just that once we go to a cubic spline, that bias is way below our precision and accuracy (under cross-validation). At least for now.

got it!

Adrian Price-Whelan (Flatiron) and I spent time this past week trying to factorize products of Gaussians into new products of different Gaussians. The context is Bayesian inference, where you can factor the joint probability of the data and your parameters into a likelihood times a prior or else into an evidence (what we here call the FML) times a posterior. The factorization was causing us pain this week, but I finally got it this weekend, in the woods. The trick I used (since I didn't want to expand out enormous quadratics) was to use a determinant theorem to get part of the way, and some particularly informative terms in the quadratic expansion to get the rest of the way. Paper (or note or something) forthcoming...

imaging asteroseismic modes on the stellar surface

Many threads of conversation over the past weeks came together today in a set of coincidences. Conversations with Bedell (Flatiron), Pope (NYU), Luger (Flatiron), and Farr (Flatiron) ranging around stochastic processes and inferring stellar surface features from doppler imaging all overlap at stellar asteroseismic p modes: In principle, with high-resolution, high-signal-to-noise stellar spectral time series (and we have these, in hand!) we should be able not only to see p modes but also see their footprint on the stellar surface. That is, directly read ell and em off the spectral data. In addition, we ought to be able to see the associated temperature variations. This is all possible because the stars are slowly rotating, and each mode projects onto the rotating surface differently. Even cooler than all this: Because the modes are coherent for days in the stars we care about, we can build very precise matched filters to combine the data coherently from many exposures. There are many things to do here.

binary stars and lots more

Today was a very very special Stars meeting, at least from my perspective! I won't do it justice. Carles Badenes (Pitt) led us off with a discussion of how much needs to be done to get a complete picture of binary stars and their evolution. It's a lot! And a lot of the ideas here are very causal. For example: If you find that the binary fraction varies with metallicity, what does it really vary with? Since, after all, stellar age varies with metallicity, as do all the specific abundance ratios. And also star-formation environment! It will take lots of data and theory combined to answer these questions.

Andreas Flörs (ESO) spoke about the problem of fitting models to the nebular phase of late-time supernovae, where you want to see the different elements in emission and figure out what's being produced and decaying. The problem is: There are many un-modeled ions and the fits to the data are technically bad! How to fix this. We discussed Gaussian-process fixes, both stationary and non-stationary. And also model elaboration. And the connection between these two!

Helmer Koppelman (Kapteyn) showed some amazing structure in the overlap of ESA Gaia data and various spectroscopic surveys (including LAMOST and APOGEE and others). He was showing visualizations in the z-max vs azimuthal-action plane. We discussed any ways it could be selection effects. It could be; it is always dangerous to plot the data in derived (rather than more closely observational) properties.

Tyson Littenberg (NASA Marshall) told us about white-dwarf–white-dwarf (see what I did with dashes there?) binaries in ESA LISA. He has performed an information-theoretic analysis for a realistic Milky Way simulation. He showed that many binaries will be very well localized; many thousands will be clearly detected; and some will get full 6-d kinematics because the chirp mass will be visible. Of course there are simplifying assumptions about the binary environments and accelerations, but there is no doubt that it will be incredible. Late in the day we discussed how you might model all the sea of sources that aren't individually detectable. But that said, everything to many tens of kpc in the MW will be visible, so incompleteness isn't a problem until you get seriously extragalactic. Amazing!