Fitting

gspice.fit() clips outliers from the training spectra in iterative passes, then estimates the mean and covariance in the working units of The model. Five kinds of model share this function: the default, constant center, S/N bins, low rank, and noise-aware.

Training spectra

A covariance of npix pixels has npix² entries, so it needs many spectra. fit skips a spectrum with more than maxbadpix bad pixels (64 by default) or with the weight 0. It also skips a spectrum whose pixels a clipping pass masks completely. $N_s$ counts the spectra that a fit uses.

With $y_s = a_s \tilde f_s$ the scaled spectrum of The model, the model of S/N bin $b$ holds

$$m_b = \frac{1}{N_s} \sum_s y_s, \qquad C_b = \frac{1}{N_s - 1} \sum_s (y_s - m_b)(y_s - m_b)^\mathsf{T} ,$$

over the spectra $s$ of the bin. With weights $w_s$, both sums are weighted, and $C_b$ divides by $\sum_s w_s - \sum_s w_s^2 / \sum_s w_s$ instead of $N_s - 1$. The weights apply to the final fit only, and the clipping passes use unit weights.

A default model needs $N_s - 1 \ge \mathrm{npix}$, and a good pixel in at least one spectrum at every pixel. Close to that minimum, $C_b$ can still fail the condition check of its precision matrix. What you get back lists the fields of the model.

Outlier clipping

An artifact in a training spectrum adds a large deviation to the sums above, and so widens the covariance at its pixels. fit therefore flags the outliers of the training spectra in clipping passes, with thresholds that fall from pass to pass (paper Alg. 4). The figure follows one spectrum with a bad column through the three default passes:

Pixels 0 to 99 of training spectrum 153, with a bad column at pixels 0 to 12

0.51.01.52.00255075100

Pass 1: |z| > 20 flags pixels 0 to 1

0.11101000255075100

Pass 2: |z| > 8 flags pixels 2 to 11

0.11101000255075100

Pass 3: |z| > 6 flags pixels 12 to 13

0.11101000255075100

• flux, then in-sample |z|
• threshold of the pass
• masked by earlier passes
• flagged by the pass

The panels show fit with its default passes on 4000 simulated spectra of 200 pixels, 30% of them with an artifact. Each pass fits a covariance without the masked pixels and computes |z| of every pixel, with nguard 20. The bad column widens the covariance of pass 1, so only its strongest pixels pass 20. Each refit is narrower, so the next pass finds more of the column. A flag also covers the pixel on each side.

Pass $n$ uses the threshold nsigma[n - 1], from nsigma=(20, 8, 6) by default, and the mask $B_n$, with $B_1$ the bad pixels. It runs four steps:

1. Interpolate the pixels of $B_n$, and scale each spectrum by $a_s$ over its pixels outside $B_n$.
2. Fit the covariance, with unit weights.
3. Compute the pixelwise z-score of every pixel of these spectra, interpolated pixels included, with nguard.
4. Flag $F_n$: each pixel with $|z|$ above the threshold, and its two neighbors.

With cumulative=True, the default, $B_{n+1} = B_n \cup F_n$. The final covariance is fitted on the last of these masks, with the weights. The returned mask is that final mask, and True at every pixel of a spectrum that fit did not use. cumulative=False uses the 2022 rule, $B_{n+1} = \mathrm{bad} \cup F_n$, with the final mask $\mathrm{bad} \cup F_{N-1} \cup F_N$ for $N$ passes. nsigma=() fits without clipping.

The clipping ceiling

A pass can flag only the outliers that its z-scores can reach. These z-scores are in-sample: each spectrum helped to fit the covariance that predicts it, so an outlier widens the covariance that judges it. For a default model, the in-sample $|z|$ of an outlier of any amplitude approaches the clipping ceiling

$$\sqrt{\frac{(N_s - 1)(N_s - N_k - 1)}{N_s}} , \qquad N_k = \max(\mathrm{npix} - 2\,\mathrm{nguard} - 1, 0) ,$$

in a spectrum of typical leverage. With a constant center, the ceiling is $\sqrt{N_s - N_k}$. No ceiling is derived for a low-rank model.

At $N_s = 2500$, npix 2400 and nguard=20, $N_k = 2359$, and the ceiling is $\sqrt{2499 \times 140 / 2500} = 11.8$. A pass at 20 then flags no pixel. fit logs a warning on the gspice logger for each pass whose threshold exceeds the ceiling. A pass above the clipping ceiling gives the remedies.

Constant center

A default model subtracts the mean of its training spectra. center=c subtracts a constant $c$ in flux units instead, such as a continuum level: $x_s = a_s (\tilde f_s - c)$. The model then estimates no mean. With unit weights, the covariance is $\sum_s x_s x_s^\mathsf{T} / N_s$, and the model needs $N_s \ge \mathrm{npix}$.

S/N bins

In scaled units, the noise of every spectrum has a variance near 1. The variation between the spectra scales with $a_s^2$, so it grows with the S/N:

110100100010⁴2050100200400

signal-to-noise ratio of the spectrum (log scale)

• variance of the noise-free flux
• variance of the noise
• edges of snr_bins=3

The steps show the variance between spectra at pixel 37, the core of a line, in scaled units squared, for 20000 simulated spectra of 100 pixels. Each step holds a tenth of the spectra by S/N. The noise variance stays near 0.711, while the variance of the noise-free flux grows from 20.03 to 2972.72. One covariance of all spectra mixes these. The vertical lines are the edges at S/N 37.8 and 59.4 that fit(..., snr_bins=3) puts between three bins of equal count.

Spectra of different S/N therefore have different covariances in scaled units, and one covariance of all of them is a mixture. snr_bins fits one model for each S/N bin.

The S/N of a spectrum is the median of $f \sqrt{\mathrm{ivar}}$ over its good pixels (gspice.snr()). snr_bins=n puts the edges between the bins at quantiles of the S/N, so that each bin holds about the same number of spectra. A sequence gives the interior edges instead.

Each bin clips and fits on its own spectra. Bin $b$ of the model equals the model that fit returns for the spectra of bin $b$ alone. At prediction, each spectrum goes to the bin of its S/N, with the edges in Model.edges. Each bin needs the spectra of a model of its kind, for example $N_s - 1 \ge \mathrm{npix}$ for a default bin. A small bin also has a low clipping ceiling, so it flags fewer outliers.

Low-rank model

A covariance of npix pixels cannot be estimated from fewer than npix + 1 spectra. A low-rank model asks for less. With rank=K, the covariance is a gspice.LowRank:

$$C = V \Lambda V^\mathsf{T} + D ,$$

with $(V, \Lambda)$ the $K$ leading eigenvectors and eigenvalues of the sample covariance $S$, and the diagonal matrix $D$ with $D_{ii} = S_{ii} - (V \Lambda V^\mathsf{T})_{ii}$. So $C$ and $S$ have the same variance at every pixel, and the correlations between pixels come from the $K$ components alone. The clipping passes use the same rank.

A low-rank model needs $N_s - 1 \ge K + 1$, a positive $K$-th eigenvalue and every $D_{ii} > 0$, but no good pixel at every pixel. It suits a training set with fewer spectra than pixels. Other routes gives its precision matrix.

Noise-aware model

A default model is a scaled model: its one covariance holds both the variation between the spectra and their noise. That fails where the noise differs from spectrum to spectrum in level or in shape. A noise-aware model separates the two.

noise=True estimates $\Sigma$, the covariance of the noise-free flux, in flux units, from the ivar of each pixel. It clips with a scaled model first, as above. At prediction, spectrum $s$ has the covariance

$$C_s = \Sigma_+ + \mathrm{diag}(\nu_s) ,$$

with $\nu_{si} = 1/\mathrm{ivar}_{si}$ at the good pixels and $1/a_s^2$ at the bad pixels. $\Sigma_+$ is $\Sigma$ with its negative eigenvalues set to 0, which biases it upward. The working vector is $x_s = \tilde f_s - \mu$, in flux units, with $\mu$ the mean of the training spectra, weighted by $a_s^2$ times the weights. rank=K also works here.

$\Sigma$ is unbiased element by element under four conditions:

1. The spectra are independent, with a common mean and a common signal covariance.
2. The noise has zero mean and the variance 1 / ivar.
3. The noise is uncorrelated between pixels and with the signal.
4. The masks, weights and ivar do not depend on the flux.

The outlier flags of the final mask break condition 4, and so does an ivar computed from the flux.

A noise-aware model has one bin and no center. It needs two spectra good at every pixel, and its clipping passes need the spectra of a scaled model. Each spectrum has its own covariance, so each needs its own factorization (Cost).

What each model needs

Model	Arguments	Spectra needed	Prediction	Finite-sample factor
Default	none	$N_s - 1 \ge \mathrm{npix}$	one $P$ per bin	derived
Constant center	center=c	$N_s \ge \mathrm{npix}$	one $P$ per bin	1
S/N bins	snr_bins=n	per bin, as for its kind	one $P$ per bin	derived for default bins
Low rank	rank=K	$N_s - 1 \ge K + 1$	one $P$ per bin	1
Noise-aware	noise=True	two good spectra at each pixel, and those of a scaled model for clipping	one factorization per spectrum	1

$P$ is the precision matrix of Numerics. A weighted model has no finite-sample factor either (Conditions). Choose a model compares the models on simulated spectra.