divide by your selection function, or multiply by it?

With Kate Storey-Fisher (San Sebastián), Abby Williams (Caltech) is working on a paper about large-angular-scale power, or anisotropy, in the distribution of quasars. It is a great subject; we need to estimate this power in the context of a very non-trivial all-sky selection function. The tradition in cosmology is to divide the data by this selection function. But of course you shouldn't manipulate your data. Instead, you could multiply your model by the selection function. You can guess which one I prefer! In fact you can do either, as long as you weight the data in the right way in the fit. I promised to write up a few words and equations about this for Williams.

double periodogram

Cole Johnston (Leuven) is in New York this week. We discussed the problem of finding oscillation modes in the photometry of stars in the presence of a large, binary-induced periodicity. What he kind-of wants is a simultaneous fitting of a flexible periodic function plus a periodogram. We did some experiments (very promising!) and discussed the elements that will come together to make this all happen. The final method will look like a double fourier transform, in which one frequency grid gets the periodic part, and the other grid gets the rest of the modes and noise.

linear regression

Valentina Tardugno (NYU) and I are looking at the NASA TESS housekeeping data: What parts of it are relevant to understanding the light curves? The weird thing is: We are asking this by asking: What housekeeping data can be reliably predicted using the light curves? Why this way? Because the light curves are higher in signal-to-noise (in general) than most channels of the housekeeping data. Today we went through all the relevant linear algebra for big linear models (which is where we are starting, of course!).

making linear algebra faster in practice

The key thing to make your code run faster is to avoid building large linear-algebra objects. For example, if you need to get the matrix product A.x, where A is a huge matrix and x is a long vector, and you only ever use the matrix A to do this one multiply by x, there is no reason to actually create A. Just create a function that evaluates A.x for any input x. That should be way faster, because of less memory allocation, and because you don't have to make the parts of the matrix that are all zeros (for example). Matt Daunt (NYU) and I discussed all this at length today, as we profiled code. This comment has some overlap with Section 9 of this paper.

CMB component separation with linear fitting

Today I sat down with Fiona McCarthy (Flatiron) to look at data-driven methods for separating cosmic microwave background data into different components. We implemented a simple polynomial regression to fit foregrounds, using (observed) difference maps as inputs (features) that are designed to contain foregrounds only. We obtained some preliminary results that looked exciting but we’ve only just started. Part of the motivation is that CNNs are hard to train, but linear combinations of image monomials are easy! I realized in all this that there are connections to the group-equivariant stuff I’ve done with Villar’s group, because we use invariants, and also to the causal inference things that Schölkopf’s group does, because we’re trying to impose some causal structure on our functions.

applying the SVD to nontrivial objects

Singular value decomposition (SVD) is a method for finding the equivalent of eigenvalues and eigenvectors for a rectangular matrix. It is what we use when we want to know the rank of a rectangular matrix, or make a low-rank matrix factorization (indeed, it is precisely what is used in principal components analysis or PCA).

The cool thing is: The method is exceedingly general; it can find the rank or a low-rank approximation to any space; it doesn't have to be a vector space exactly. It just has to obey certain algebra rules. So in my work with Soledad Villar (JHU) we use it to find a basis to represent all the linearly independent geometric images (grids of tensors) possible subject to constraints (like symmetries). I wrote words about using the SVD in this context in our nascent paper. Here is some example output of my SVDs in eye-candy form:

how many equivariant linear functions are there?

O. M. G. As my loyal reader knows, Soledad Villar (JHU) and I are trying to build a replacement for convolutional neural networks that can handle geometric objects (scalars, vectors, pseudovectors, tensors of any order, so on) and that can create functions that exactly (or approximately if you like) obey the symmetries of classical physics (rotation, translation, parity, boost, maybe even gauge). Our method produces polynomial functions of images (functions where both the input and the output are images), making use of convolutions, outer products, index contractions, index permutations, cross-products, pooling, and so on.

Meanwhile, Ben Blum-Smith (NYU, JHU) has been doing (scary to me) group theory stuff in which he has been computing the number of unique polynomial functions of images of fixed polynomial degree that are possible, given image inputs and outputs (of some tensor orders), when those polynomial functions obey the symmetries of classical physics. And he has results! He can tell us how many unique linear and quadratic, say, functions there are of vector images (say) that output vector images. It's a formula that depends on the image size and the degree of the polynomial.

Today Villar and I had a breakthrough: We used our geometric generalization of convolutions to produce all possible linear functions of small images, deduplicated the results using a singular value decomposition, and (thereby) counted all linearly independent group-equivariant linear functions there are that go from vector images to vector images. And our results agree with Blum-Smith's formula. So we may actually have a complete basis for all image functions that could ever exist in the context of classical physics?

making Fourier fitting super fast

At the request of Conor Sayres (UW), I have been looking at distortion patterns in the SDSS-V Focal Viewing Camera (FVC), which is the part of the system that looks at whether the focal-plane fiber robots are where they need to be. The distortions are extremely bad; they are large in amplitude and vary on extremely small scales on the focal plane. So I have to fit an extremely flexible model. Here are my comments:

First, you should use mixtures of sines and cosines for problems like this. Not polynomials! Why? Because sines and cosines do not blow up at the edges.

Second, you should punk fast Fourier transform (FFT) codes to speed up your regressions. I wrote code to do this, which wraps the finufft best-in-class non-uniform FFT code in scipy.sparse linear-algebra code. This wrapping makes the FFT operators into linear-algebra operators and permits me to do solve() operations. That move (wrapping FFT in linear algebra) sped up my code by factors of many!

GPRV, day 3

Day 3 of GPRV continued great! There were a few talks and discussions of very young stars that got everyone in the room quite excited, from Di Maio (INAF), Suárez Mascareño (IAC), and Nielsen (Oxford). The activity signals are huge, but the planets are extremely interesting, so how do we approach this? Tons of of observing time? Cleverness? Give up? Of course I think it is so important to understand how planetary systems form and evolve, I would be willing to spend the telescope time.

In the morning, Luger (Flatiron) gave a seminar and then a tutorial about modeling stellar surfaces and predicting spectroscopic quantities. The tutorial was fun; his code Starry does everything an astronomer could want, and beautifully (and, of course, blazingly fast). We had fun playing with it in a group hacking session.

is an eigenvector a vector?

I spent time today talking to Kate Storey-Fisher about features to use in her cosmological regression projects. The point is that we are only considering features that have well-defined, coordinate-free meanings, because we are trying to do regressions that are invariant to coordinate transformations. These features include scalars, vectors, and tensors, which we contract into scalars. But what can you do with a tensor? At order 2, a tensor has a trace (its self-contraction); it can be contracted with two vectors; it has eigenvalues and eigenvectors. The eigenvalues are classical scalars; good! But are the eigenvectors classical vectors? No, they aren't, because they don't have signs. What can you do with them? I have some theories...

feature engineering for dark-matter halos

Kate Storey-Fisher and I spoke about adding eigenvalues (scalars) and eigenvectors (vectors) to our geometric features of dark-matter halos. For regression! But the problem with this is that the eigenvectors have ambiguous sign; is that an issue? Yes! They are descriptions of an order-2 tensor, not order-1 vectors. Hmmm. We also spoke about whether the feature engineering should be deterministic, or chosen by hand.

convert an L1 minimization into a linear program

Okay, linear programs are amazing things! And okay, L1 minimization and related sparse methods are magical. But can these be related? I always thought “no!” After all, the L1 norm involves an absolute value operator, and that isn't linear; linear programs are linear (it's in the name!). But today Soledad Villar (JHU) blew me away with the following (retrospectively simple) observation: If you augment any variable that you are going to L1 (that variable could be a parameter or a residual) with an auxiliary variable that is constrained to be greater than both the original variable and also the variable multiplied by negative one, then you can set up linear programs such that the auxiliary variables become the absolute values of their associated varaibles. And the optimization performed by the linear program becomes an L1 optimization. It's beautiful, we coded it up today in a few different contexts, and it worked!

linear programming FTW

As my loyal reader knows, I am violating all of my principles and engaging in a linear-regression version of a symbolic regression project to test how ideas in dimensional analysis (or units equivariance) might impact machine learning. I have been struggling to get a sparse regression working, because when problems get large, optimizing combined L1 and L2 losses can be sticky and tricky. But Soledad Villar (JHU) saved me today by pointing out that in the over-parameterized regime (when you have more terms in your linear regression than data), you can do a sparse regression with a cleverly designed linear program! Woah, we coded it up and it Just Worked (tm)! We can get the exactly correct total mechanical energy expression in our toy problem with a very small amount of training data, far less than we needed when we were using L2 as our objective. Far less.

a conditional form of domain adaptation?

I have been discussing with Soledad Villar (JHU) projects related to domain adaptation (also with Thabo Samakhoana and Katherine Alsfelder), in which you have two different instruments (say) taking data and you want to find the transformation between the instruments such that the data are the same from the two sources. The modification we want (or need) to make is to create a conditional version of this: In SDSS-V, there are two observatories; they take very similar (but not identical) data. What are the transformations that make the data identical? The problem is: The two observatories also observe different stars on average (because they see different parts of the Galaxy). So we need to find the transformations that make the data identical, conditional on other data (like the ESA Gaia data) that we have for the stars. Great problem, and we came up with some non-elegant solutions. Are there also elegant solutions?

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.

predicting wavelength calibration from housekeeping data

I did some of my favorite thing today, which is fitting flexible models. The context was my attempt to predict the wavelength solution for the SDSS-IV BOSS spectrographs using only the housekeeping data, like the state of the telescope, temperatures, and so on. It doesn't work accurately enough, or at least not with the housekeeping data I've tried so far. It looks like there might be a hysteresis or a clank or something like that. If this is right, it bodes poorly for reducing the number of arcs we need to take in SDSS-V, which is supposed to move fast and not break things.

But all that said, I still have one card left to play, which is to see if we can look at sky lines in science frames and learn enough from sky lines—plus the historical behavior of all arcs ever taken—such that we can lock everything else down without taking an arc for every visit.

Dr Dou Liu

Today I had the pleasure of sitting on the PhD defense of Dou Liu (NYU), who has been working on AGN in the centers of galaxies, using MANGA data from SDSS-IV. The part of Liu's thesis that is most exciting to me (perhaps not surprisingly) is the technical chapter, in which he finds a new method for combining irregularly dithered integral-field-unit spectroscopy exposures into a full data cube, with sky coordinates on two axes, and wavelength on the third. In this final data cube, his method gets much better final resolution, and lower pixel-to-pixel covariances in the noise, relative to the standard pipelines. His trick? He has generalized spectro-perfectionism (invented for spectral extraction) to the multi-dimensional spectral domain. It's beautiful stuff, and has implications for all sorts of imaging and spectroscopy projects going forward. Congratulations Dr. Liu and thank you!

my own special FMM algorithm

I was in a location with no internet and no computing, so I spent an hour or so writing down what I think is the fast multipole method. It involves building an octtree, balanced in volume (not necessarily point content), and computing recursively the multipole amplitudes in all nodes (starting with the points in the leaves). Once that is done, at evaluation time, you do different things at different levels, depending on the radius out to which things are computed exactly. One thing I'm interested in is: Can you simplify that evaluation if you, say, build multiple trees?

PCA and convex optimization

Soledad Villar (JHU) and I talked about bi-linear problems today, in the context of instrument calibration and computer vision. We looked at the kinds of optimizations that are involved in these problems. She showed me that, if you think of PCA as delivering a projection operator (that is, not the specific eigenvectors and amplitudes, but just the projection operator you would construct from those eigenvectors), that projection operator can be derived from a convex optimization in which the objective is purely mean squared error. That was news to me!

k near n?

Teresa Huang (JHU) has a nice paper (with Villar and me) that shows the risk and regularization of linear regression involving PCA. We discussed it more today, in particular whether we can say more about the regime in which the PCA dimensionality reduction (to k dimensions) doesn't do much (because k is close to the number of data points n). We think we can, because the Marchenko-Pastur distribution of eigenvalues is so skew: Cutting off even one small eigenvalue (k=n-1) can be useful!