The model

GSPICE models an ensemble of spectra on a shared wavelength grid as draws from a multivariate Gaussian distribution and predicts masked or corrupted pixels from the remaining pixels by Gaussian conditioning (Finkbeiner et al., 2025).

A spectrum as one draw from a Gaussian

A model must say which spectra are likely and which are not. Take many stars observed on one wavelength grid. At pixel $i$, the flux differs from star to star, so pixel $i$ is one variable, measured once in every spectrum. A spectrum of npix pixels is then one point in npix dimensions.

GSPICE describes where these points lie with a multivariate Gaussian distribution. It has two parameters: a mean spectrum, and a covariance $C$ between every pair of pixels. gspice.fit() estimates both from the training spectra, and Fitting explains how. What GSPICE needs lists what this asks of your spectra.

The spectra need not be Gaussian. A population of stars is far from a Gaussian distribution. The conditional estimates of a Gaussian fitted to it still predict its pixels well (Finkbeiner et al., 2025).

Covariance: lines vary together

One pixel can tell you about another only if the two vary together. Take two lines of one element. A star with a deep line at pixel $i$ also tends to have a deep line at pixel $i'$. Across the spectra, the fluxes at $i$ and $i'$ then rise and fall together.

The covariance $C_{ii'}$ measures this. It is positive where the two fluxes tend to lie above their means together. It is negative where one tends to lie above its mean while the other lies below, and 0 where they are unrelated. $C_{ii}$ is the variance at pixel $i$. The correlation $C_{ii'} / \sqrt{C_{ii} C_{i'i'}}$ puts this on a scale from −1 to 1.

The explainer below maps the correlation of a fitted model, and plots the training spectra at any two of its pixels against each other.

Point at a cell of the map, or tap it, to choose two pixels. With the map focused, the arrow keys move the cell, and Shift with an arrow key moves it 10 pixels.

03060901200306090120

0.20.30.40.50.60.70.20.30.40.50.60.70.8pixel 12pixel 83

0.50.60.70.80.91.003060901201283

Pixels 12 and 83: correlation −0.98.

• −1
• −0.5
• 0
• 0.5
• 1

• training spectrum
• 1σ and 2σ ellipses of the model
• mean spectrum

The map shows the correlation of each pair of pixels in a model fitted to 4000 simulated spectra of 120 pixels, seed 11. The scatter shows 100 of those training spectra at the two pixels.

Conditioning on two pixels, then on many

A prediction answers this question: given the flux that you observed at some pixels, what is the flux at the others? Start with two pixels, $i$ and $i'$, and measure each flux from its mean, so that both means are 0. If you observe $x_{i'}$ at pixel $i'$, the flux at pixel $i$ has a Gaussian distribution again (Rasmussen & Williams, 2006):

$$x_i \mid x_{i'} \sim \mathcal{N}\!\left(\frac{C_{ii'}}{C_{i'i'}}\, x_{i'},\; C_{ii} - \frac{C_{ii'}^2}{C_{i'i'}}\right) .$$

Its mean moves away from 0 in proportion to $x_{i'}$. Its variance is smaller than $C_{ii}$, the more so the stronger the correlation. In the explainer below, you fix the flux of one pixel and see the distribution of the other.

Move the slider, or drag across the left plot, to fix the flux of one pixel. The plot on the right shows the distribution of the other pixel.

0.20.40.60.81.00.50.60.70.80.9pixel 14, fixedpixel 12

Flux of pixel 14:

0.825

Next pairPair 1 of 4: correlation 0.97

Pixel 14 fixed at 0.825. Pixel 12: 0.713 ± 0.031, against 0.541 ± 0.118 without pixel 14.

• 1σ and 2σ ellipses; distribution without the fixed pixel
• fixed flux
• conditional mean and distribution
• ±1σ
• ±2σ

The figure compares two pixels of a model fitted to 4000 simulated spectra of 120 pixels, seed 11. Pixel 12 has the largest variance, and the 4 pairs span correlations from about 1 to about −1. The prediction is that of gspice.infill from the fixed pixel alone.

The same holds for any number of pixels. Split the pixels into two sets. $G$ holds the pixels to predict or to leave out, and $O$ holds the conditioning pixels, whose values the prediction uses. With $x$ measured from the mean again, the conditional distribution of $x_G$ is

$$x_G \mid x_O \sim \mathcal{N}\!\left(C_{GO} C_{OO}^{-1} x_O,\; C_{GG} - C_{GO} C_{OO}^{-1} C_{OG}\right) .$$

$C_{GO}$ is the block of $C$ with the rows $G$ and the columns $O$, and the other blocks follow the same rule. The conditional mean is the prediction. It is the average of all spectra in the model that agree with the conditioning pixels, and a linear function of $x_O$. Where $G$ holds all pixels, the result is the prior $\mathcal{N}(0, C)$.

In the explainer below, the masked block is $G$, and the square beside the spectrum marks the rows and columns of $G$ in the correlation of $C$.

Drag the gray block, or drag its edges to resize it. With the plot focused, the arrow keys move the block, and Shift with an arrow key resizes it.

0.20.40.60.81.00306090120

03060901200306090120

Show flux without noiseNext spectrumSpectrum 1 of 8

Pixels 50 to 69 masked (20). Mean σ 0.019. RMS error 0.003.

• observed flux
• prediction
• ±1σ
• ±2σ
• flux without noise

• −1
• −0.5
• 0
• 0.5
• 1

The prediction is that of gspice.infill with a model fitted to 4000 simulated spectra of 120 pixels. The ±1σ and ±2σ bands include the noise of the observed flux. The square shows the correlation of each pair of pixels. The prediction reads the marked rows and columns of the masked pixels.

A prediction can report only some pixels of $G$. The others are marginalized: their values do not enter, and the conditional covariance of the reported pixels includes the uncertainty of the others. Prediction shows how each function chooses $G$. The code computes the same distribution from the inverse of $C$, as The precision form explains.

The uncertainty of a prediction

A prediction needs an uncertainty. The conditional covariance $C_{GG} - C_{GO} C_{OO}^{-1} C_{OG}$ provides it. It depends on which pixels you observed, but not on their values. Its diagonal holds the conditional variance of each predicted pixel.

The training spectra carry their noise, so $C$ holds the noise as well as the variation between the spectra. Every predicted variance therefore includes the noise of the predicted pixel. The prediction describes the observed flux at the pixel, not its noise-free flux.

Two checks compare the predicted uncertainty with the data:

• The z-score of one pixel is its flux minus the prediction, divided by the square root of the conditional variance.
• The χ² of several pixels uses their full conditional covariance (χ² and its degrees of freedom).

Both have a known distribution only if the model is right. A model estimated from a finite number of spectra predicts new spectra less well than it says (Held-out uncertainty). Diagnostics measures both effects on your spectra.

Working units

Spectra differ in brightness and in signal-to-noise ratio, and some pixels are bad. One covariance of the raw flux would mix these differences with the variation that predicts pixels. The model therefore describes the working vector $x_s$ of spectrum $s$, not its flux. It computes $x_s$ in three steps:

1. Interpolate over the bad pixels: f̃_s (flux)

0.20.40.60.81.00255075100

2. Multiply by a_s = 34.88: a_s f̃_s (scaled units)

2040600255075100

3. Subtract the mean of the training spectra: x_s = a_s f̃_s − m_b (scaled units)

−40−35−30−25−20−15−100255075100

• observed flux
• interpolated flux
• mean of the training spectra
• bad pixels

The panels trace a simulated spectrum of 100 pixels at a signal-to-noise ratio of 31, with pixels 56 to 68 marked bad. The model is fitted to 4000 simulated spectra. Step 1 replaces the bad pixels by a straight line. Step 2 brings the noise to a variance near 1. Step 3 leaves the variation about the mean.

1. Interpolate the flux linearly over the bad pixels, and hold it constant beyond the first and the last good pixel. This gives $\tilde f_s$. A bad pixel is never read, so its flux may be NaN.
2. Multiply by $a_s$, the square root of the mean of ivar over the good pixels of the spectrum, in flux⁻¹. This gives scaled units.
3. Subtract the mean $m_b$ of the training spectra in the S/N bin $b$ of the spectrum. There is one bin by default.

So

$$x_s = a_s \tilde f_s - m_b .$$

In scaled units, the noise at pixel $i$ has the variance $a_s^2 / \mathrm{ivar}_{si}$. That variance is 1 where ivar is equal at every pixel of the spectrum, whatever its S/N. A variance in scaled units divided by $a_s^2$ is a variance in flux units. Continuum normalization says what the scaling means for your spectra.

Two kinds of model use other working vectors:

Model	Working vector	Units
Default, or with S/N bins	$x_s = a_s \tilde f_s - m_b$	scaled
Constant center $c$	$x_s = a_s (\tilde f_s - c)$	scaled
Noise-aware	$x_s = \tilde f_s - \mu$	flux

Fitting describes both. The public functions take flux and return flux, flux² or unit-free values. The working units appear only in gspice.scale(), gspice.unscale() and the fields of gspice.Model.