Analytic marginalization of absorption line continua

TL;DR

This paper introduces an analytic method for marginalizing over continuum parameters in absorption line spectroscopy, enabling faster, more accurate, and automated analysis of large datasets compared to traditional numerical methods.

Contribution

It develops a closed-form analytic marginalization technique for linear continuum models, improving speed and accuracy in automated spectral analysis.

Findings

Analytic marginalization speeds up spectral analysis.

Improves accuracy of absorption line parameter recovery.

Facilitates automated continuum placement in large datasets.

Abstract

Absorption line spectroscopy is a powerful way of measuring properties of stars and the interstellar medium. Absorption spectra are often analyzed manually, an approach that limits reproducibility and which cannot practically be applied to modern datasets consisting of thousands or even millions of spectra. Simultaneous probabilistic modeling of absorption features and continuum shape is a promising approach for automating this analysis. Existing implementations of this approach use numerical methods such as Markov chain Monte Carlo (MCMC) to marginalize over the continuum parameters. Numerical marginalization over large numbers of continuum parameters is too slow to be convenient for exploratory analysis, can increase the dimensionality of an inference problem beyond the capacity of simple MCMC samplers, and is in general impractical for the analysis of large datasets. When continua are parameterized as linear functions such as polynomials or splines, it is possible to reduce continuum parameter marginalization to an integral over a multivariate normal distribution, which has a known closed form. In addition to speeding up probabilistic modeling, analytic marginalization makes it trivial to marginalize over continuum parameterizations and to combine continuum description marginalization with optimization for absorption line parameters. These new possibilities allow automatic, probabilistically justified continuum placement in analyses of large spectroscopic datasets. We compare the accuracy to within which absorption line parameters can be recovered using different continuum placement methods and find that marginalization is in many cases an improvement over other methods. We implement analytic marginalization over linear continuum parameters in the open-source package amlc.

Figures

Equations

$$\mathbf{y}(\theta )$$

$$\mathbf{y}$$

$$\equiv \mathbf{L}\left({\mathbf{\mu }}_b(\theta )+{\mathbf{D}}_{\theta }{\mathbf{\mu }}_m(\theta )+\mathbf{B}\mathbf{c}\right)+\mathbf{\epsilon }.$$

$$\mathbf{B}$$

$$p_n(\mathbf{c})=\mathcal{N}\left(0,\mathbf{\Lambda }\right)\text{ (normal) and }{p}_u(\mathbf{c})=\prod _{i=1}^{P+Q}{Z}_i^{-1}\text{ (uniform)},$$

$${\hat{\mathbf{c}}}_{n/u}$$

$$\mathbf{r}$$

$${\mathbf{C}}_n$$

$${\mathbf{C}}_u$$

$$p_n(\mathbf{y}|\theta ,\mathbf{L},\mathbf{B},\mathbf{K},\mathbf{\Lambda })$$

$$=(2\pi {)}^{-\frac{M}{2}}\det (\mathbf{K}{)}^{-\frac{1}{2}}\det (\mathbf{\Lambda }{)}^{-\frac{1}{2}}\det ({\mathbf{C}}_n{)}^{-\frac{1}{2}}\exp \left[-\frac{1}{2}{\mathbf{r}}^T{\mathbf{K}}^{-1}\left(\mathbf{r}-{\hat{\mathbf{r}}}_n\right)\right],$$

$${\hat{\mathbf{r}}}_{n/u}$$

$$p_u(\mathbf{y}|\theta ,\mathbf{L},\mathbf{B},\mathbf{K})$$

$$=\left(\prod _{i=1}^{P+Q}{Z}_i^{-1}\right)(2\pi {)}^{-\frac{M-(P+Q)}{2}}\det (\mathbf{K}{)}^{-\frac{1}{2}}\det ({\mathbf{C}}_u{)}^{-\frac{1}{2}}\exp \left[-\frac{1}{2}{\mathbf{r}}^T{\mathbf{K}}^{-1}\left(\mathbf{r}-{\hat{\mathbf{r}}}_u\right)\right].$$

$$\begin{array}{rl}\nabla \log (p)& (\mathbf{d}(\theta ))=\left({\mathbf{L}}^T{\mathbf{K}}^{-1}\left(\mathbf{r}-{\hat{\mathbf{r}}}_{n/u}\right)\right)\odot \left({\mathbf{B}}^{\prime }\hat{\mathbf{c}}+{\mathbf{\mu }}_m\right)\\ & -\frac{1}{2}\left(\left({\mathbf{C}}_{n/u}^{-1}{\mathbf{B}}^{\prime T}\right)\odot \left({\mathbf{B}}^T{\mathbf{L}}^T{\mathbf{K}}^{-1}\mathbf{L}\right)+\left({\mathbf{C}}_{n/u}^{-1}{\mathbf{B}}^T{\mathbf{L}}^T{\mathbf{K}}^{-1}\mathbf{L}\right)\odot {\mathbf{B}}^{\prime T}\right)1,\end{array}$$

$${\mathbf{B}}^{\prime }$$

$$=\sum _{i=1}^N\left[\begin{array}{cc}{\mathbf{J}}^{i,i}{\mathbf{A}}_{\mathbf{m}}& 0\times {\mathbf{A}}_{\mathbf{b}}\end{array}\right]$$

$$=\left[\begin{array}{cc}{\mathbf{A}}_{\mathbf{m}}& 0\end{array}\right],$$

$$\nabla \log (p)({\mathbf{\mu }}_m(\theta ))$$

$$\nabla \log (p)({\mathbf{\mu }}_b(\theta ))$$

Full text

methods: statistical, techniques: spectroscopic

††software: emcee (Foreman-Mackey et al., 2013), matplotlib (Hunter, 2007), numpy (van der Walt et al., 2011), scipy (Jones et al., 2001)

1 Introduction

Absorption lines contain information on the composition and properties of interstellar matter (ISM) and stellar atmospheres. To extract this information, it is necessary to decompose the spectrum into absorption features and the intrinsic flux, typically referred to as the continuum, produced by the illuminating background source towards which the absorption is seen. The most common way of doing this separation has been manually finding regions in a spectrum that do not contain absorption features, fitting a function to these regions, and using this function to interpolate over the absorption features. Given the longevity and popularity of this approach, it is clear that it can produce acceptable results. It does, however, have two important weaknesses. The first is that every spectrum must be examined and interacted with by a human. This cannot efficiently be done for datasets containing thousands or even millions of spectra. The second is that it is unlikely that the absorption parameter estimator this procedure implicitly defines uses data efficiently. There is variance between analyses done by different humans as well as between analyses done by the same human at different times. If there is a subset of analysts whose estimates are the most accurate and precise, then the estimates of the rest are using the available data inefficiently.

An alternative approach is to infer absorption line and continuum parameters simultaneously. To improve the accuracy of the inferred absorption line parameters, it can be useful to marginalize over, rather than fit for, the continuum parameters. This has been done in packages meant for the analysis of absorption lines from both the ISM (BayesVP, Liang et al. 2018) and stellar atmospheres (Starfish, Czekala et al. 2015, and sick, Casey 2016). In these packages, continuum parameter marginalization is done numerically, using Markov chain Monte Carlo (MCMC). As the authors of two of these packages point out, including large numbers of continuum parameters in MCMC sampling leads to long convergence and autocorrelation times. To keep the number of continuum parameters low, these packages either do not support (BayesVP) or advise against (sick) including continuum parameters when simultaneously analyzing multiple spectral segments.

Even when there are few continuum parameters, using MCMC to analyze a modern spectroscopic dataset consisting of thousands or even millions of spectra will be computationally demanding. Probablistic analyses of comparable numbers of photometric observations (e.g. Green et al., 2015; Gordon et al., 2016) require months or years of computation time. A single spectrum contains orders of magnitude more data points than a single set of photometric observations, which comes with an at least linearly proportional increase in required computation time.

In the packages discussed above and in much of the absorption line analysis literature, the continuum is assumed to be a low order polynomial or spline. While these are non-linear functions of wavelength, they are linear functions of the polynomial or spline coefficients. This linearity means that if some additional assumptions hold, it is possible to marginalize over these coefficients analytically.

Analytic marginalization has several advantages over numerical marginalization. First, it can speed up MCMC-based inference for absorption line parameters by reducing the dimensionality of the problem. Second, the ability to evaluate the continuum parameter-marginalized likelihood function and its gradient makes it possible to do optimization, rather than sampling, for the absorption parameters while still keeping the robustness provided by marginalizing over continuum parameters. Finally, it makes marginalization over different possible continuum parameterizations computationally trivial—simply sum together continuum-marginalized likelihoods that assume different parameterizations. Marginalization over parameterizations allows an even greater degree of automation and systematization of absorption line inference.

The assumptions required for this particular form of analytic marginalization are: that the continuum can be expressed as a linear function (not necessarily a polynomial or spline); that the priors on the parameters of this linear function are either improper uniform or (multivariate) normal; and that residuals from the model are normally distributed. If these assumptions hold, then, given a model for the absorption, the posterior probability distribution function of the continuum parameters is itself a multivariate normal distribution. The result of marginalizing over the continuum parameters is simply an update to the covariance matrix of the model residuals. Conceptually, when a set of absorption line parameters is specified, the continuum parameters can be treated as additive, rather than multiplicative, linear nuisance parameters. An explanation of marginalization over additive linear nuisance parameters in an astronomical context is given in Luger et al. (2017). This approach to marginalizing over multiplicative linear nuisance parameters has already been used for several astronomical applications, for example for analyzing sparsely sampled radial velocity measurements (Price-Whelan et al., 2017).

Models for absorption line spectra have features, such as the presence of a line spread function (LSF), which should be accounted for to more efficiently compute marginalized likelihoods and likelihood gradients. In this work, we derive expressions for these quantities that account for these features. This derivation is given in Section 2. We have created a package, amlc111amlc is available at https://github.com/ktchrn/amlc., for evaluating these expressions. The package is described in Appendix A. The performance of continuum parameter and parameterization marginalization is explored in Section 3. We discuss strengths and weaknesses of this method in Section 4 and conclude in Section 5.

2 Assumptions and formalism

We assume the following model for a spectrum ${\bf{y}}$ given parameters $\theta$, ${\bf{m}}$, and ${\bf{b}}$:

[TABLE]

The background source emits a continuum, which is expressed as the sum of a possibly non-linear term, $\boldsymbol{\mu}_{m}(\theta)$, and a linear combination of basis elements ${\bf{a}}_{m,i}$ with coefficients $m_{i}$. Intervening matter absorbs part of this continuum with transmittance function ${\bf{d}}(\theta)$. The absorption happens independently at each wavelength. This is indicated by the elementwise product $\odot$ between the transmittance and continuum. Foregrounds, such as sky lines or instrumental artifacts are, like the continuum, expressed as the sum of a possibly non-linear term, $\boldsymbol{\mu}_{b}(\theta)$, and a linear combination of basis elements ${\bf{a}}_{b,i}$ with coefficients $b_{i}$. The resulting spectrum is convolved with an LSF ${\bf{L}}$ and observed. $\boldsymbol{\varepsilon}$ are the residuals between the observed ${\bf{y}}$ and the LSF-convolved spectrum and are assumed to be normally distributed with mean zero and covariance matrix ${\bf{K}}$. The length of the observed spectrum is $M$, the length of the pre-LSF model spectrum is $N$, the number of continuum basis elements is $P$, and the number of foreground basis elements is $Q$.

Collecting the multiplicative (continuum) and additive (foreground) basis elements ${\bf{a}}_{m,i}$ and ${\bf{a}}_{b,i}$ into matrices ${\bf{A_{m}}}$ and ${\bf{A_{b}}}$ and converting the transmittance vector ${\bf{d}}(\theta)$ into the diagonal matrix ${\bf{D}}_{\theta}\equiv\text{diag}\left({\bf{d}}(\theta)\right)$,

[TABLE]

In the second expression, ${\bf{B}}$ and ${\bf{c}}$ are defined as:

[TABLE]

We consider two possible priors for the nuisance parameter vector ${\bf{c}}$, a multivariate normal distribution with mean zero and covariance matrix $\boldsymbol{\Lambda}$ and an improper uniform distribution:

[TABLE]

where $Z_{i}$ can be any positive real number.

2.1 Conditional probability of the nuisance parameters

For both priors, the conditional distribution of ${\bf{c}}$ at fixed $\theta$ is proportional to a multivariate normal distribution. The mean ${\bf{\hat{c}}}$ of this normal distribution is

[TABLE]

where ${\bf{r}}$ is the vector of residuals

[TABLE]

and ${\bf{C}}_{n/u}$ is

[TABLE]

if the prior on ${\bf{c}}$ is normal, and

[TABLE]

if the prior on ${\bf{c}}$ is uniform. The covariance matrix of the conditional distribution of ${\bf{c}}$ is ${\bf{C}}_{n/u}^{-1}$.

The conditional distribution of ${\bf{c}}$ can be used for visualization and predictive checks. The mean of the conditional distribution is also its mode, so ${\bf{L}}{\bf{B}}{\bf{\hat{c}}}$ is the best-fit model for ${\bf{y}}$ at a given value of $\theta$. Samples drawn from the conditional distribution of ${\bf{c}}$ can be used to visualize the effect and extent of nuisance parameter variation.

2.2 Marginal likelihood

Assuming the normal prior $p_{n}({\bf{c}})$, marginalizing over ${\bf{c}}$ following e.g. Luger et al. (2017) or Rasmussen & Williams (2006) gives

[TABLE]

where

[TABLE]

If we instead assume the improper prior $p_{u}({\bf{c}})$,

[TABLE]

The marginal likelihood $p_{u}$ will be proper if ${\bf{C}}_{u}$ is positive definite, which will be the case when ${\bf{L}}{\bf{B}}$ is full rank and $M\geq P+Q$. The marginal likelihood $p_{n}$ is always proper because ${\bf{C}}_{n}$ is always positive definite. ${\bf{C}}_{n}$ is always positive definite because $\boldsymbol{\Lambda}^{-1}$ is always positive definite and ${\bf{B}}^{T}{\bf{L}}^{T}{\bf{K}}^{-1}{\bf{L}}{\bf{B}}$ is always at least positive semi-definite.

2.3 Gradients

We give expressions for the gradients of $\log(p_{n})$ and $\log(p_{u})$ with respect to ${\bf{d}}(\theta)$, $\boldsymbol{\mu}_{b}(\theta)$, and $\boldsymbol{\mu}_{m}(\theta)$. The gradient of $\log(p)$ with respect to the parameters $\theta$ can be obtained by evaluating each of these gradients, computing the Jacobians of ${\bf{d}}(\theta)$, $\boldsymbol{\mu}_{b}(\theta)$, and $\boldsymbol{\mu}_{m}(\theta)$ with respect to $\theta$, and applying the chain rule.

The gradient of $\log(p)$ with respect to ${\bf d}(\theta)$ is

[TABLE]

where ${\bf{1}}$ is a column vector of ones of length $P+Q$. ${\bf{B^{\prime}}}$ is the sum of derivatives of ${\bf{B}}$ with respect to each element of ${\bf{d}}(\theta)$:

[TABLE]

where ${\bf{J}}^{i,i}$ is a square matrix whose $(i,i)$-th entry is 1 and whose other entries are all 0. The first row of Equation 15 is the gradient of the argument of the exponentials in Equations 10 and 13. The second row is the gradient of $\log\left(\det\left({\bf{C}}_{n/u}\right)\right)$.

The gradient of $\log(p)$ with respect to $\boldsymbol{\mu}_{m}(\theta)$ is

[TABLE]

and the gradient of $\log(p)$ with respect to $\boldsymbol{\mu}_{b}(\theta)$ is

[TABLE]

3 Practical test cases

Does marginalization over continuum parameters, analytic or numerical, have benefits beyond having a clear probablistic justification? We consider two metrics: how marginalizing over, instead of fitting for, continuum parameters and parameterizations affects the error to within which absorption line parameters can be measured; and how marginalizing analytically, instead of numerically, affects the speed of MCMC-based inference.

For the error metric, we consider two possible cases: one in which the continuum parameterization is known and one in which it is necessary to choose a continuum parameterization from a set of possibilites. The suggested method is marginalizing over continuum parameters, in the first case, and marginalizing over continuum parameters and parameterizations, in the second. When the parameterization is known, all methods we consider are equally effective at high signal-to-noise ratios (SNRs) but marginalization is consistently more robust at low SNRs. When the parameterization is not known, marginalization over parameterizations has the lowest over-all error rate among the methods we consider. Furthermore, its error rate is close to that of parameter marginalization given the correct parameterization. These two cases are examined in Sections 3.1 and 3.2.

To assess the relative speed of analytic and numerical marginalization, we examine how two metrics change as the complexity of an inference problem increases: the number of iterations required for MCMC to converge and the number of independent samples generated per unit time. The basic problem is analyzing a single line on a continuum whose parameterization is known. To build up more complicated problems, we add more spectral segments each of which has its own continuum and contains another absorption line. All of these absortion lines share widths and central velocities but have independent column densities. Problems with this structure arise when analyzing multiple lines from a single species or from multiple species that can be assumed to share a common component structure. The convergence speedup from using analytic marginalization is dramatic, reaching a full order of magnitude difference in the number of required iterations with as few as three spectral segments. Analytic marginalization yields more independent samples per unit time when there are multiple spectral segments with high-order continua. For example, when there are six spectral segments with 3rd order continua, analytic marginalization is three times faster than numerical marginalization. When there are few spectral segments, analytic marginalization is slightly slower or of comparable speed to numerical marginalization. These metrics are examined in Section 3.3.

3.1 Marginalization over parameters

Here, we consider the problem of measuring the column density of a single well-resolved, unsaturated absorption line superimposed on a continuum whose parameterization is known but whose parameters are not known. To do this, we generate spectra containing an absorption line with fixed absorption parameters but with varying continuum parameterizations, continuum parameters, and SNRs. The continuum parameterizations we consider are polynomials of order 0, 1, and 2 and the continuum-level SNRs we use are between 5 and 25. We generate 1000 spectra at each combination of parameterization and SNR and measure the column density of the absorption line in each spectrum. From the recovered column densities, we compute the root mean square error (RMSE) of the base ten logarithm of the column density ($\log_{10}{\rm N}$). We use the logarithm of the column density because its RMSE scale-free. In particular, physical constants such as oscillator strengths cancel in the $\log_{10}{\rm N}$ RMSE calculation.

We measure the column density using four methods: (1) supply the correct continuum parameters and only fit for the absorption line parameters; (2) simultaneously fit for continuum and absorption line parameters; (3) analytically marginalize over continuum parameters while fitting for absorption line parameters; and (4) use the absorption line parameters recovered using method (1) to define a line-free spectral region, fit continuum parameters just to this region, and with those continuum parameters fit for the absorption line parameters. The first method provides a lower limit on the RMSE of $\log_{10}{\rm N}$ as a function of SNR. The second and third methods are two possible ways of automatically modeling the continuum. The fourth method approximates the actions of a human manually analyzing a spectrum. We assume the human can correctly estimate the continuum by eye, correctly estimate the best-fit absorption line profile by eye given this continuum, and use this profile to determine which part of the spectrum is not affected by the line.

The RMSEs obtained using these methods are shown in Figure 1. Above an SNR of 10-15, all methods where the continuum is not a priori known yield equal results. At and below that SNR range, marginalization has a lower RMSE than both simultaneous fitting and the human-like analysis. This advantage becomes greater as the continuum parameterization becomes more complex.

3.2 Marginalization over parameterizations

Next, we consider a problem where there is still a single well-resolved and unsaturated absorption line but where it is only known that the continuum belongs to a family of possible continuum parameterizations. As in the previous section, we consider three possible continuum parameterizations: 0th, 1st, and 2nd order polynomials. The approach, simulating spectra, measuring $\log_{10}N$ for each simulated spectrum, and computing the RMSE of $\log_{10}N$, is also the same. However, we consider a different range in continuum-level SNR: 10 to 100.

We measure the column density using four methods: (1) supply the correct parameterization and marginalize over its parameters; (2) assume the most complicated of the three parameterizations and marginalize over its parameters; (3) use an iterated likelihood ratio test to choose a parameterization and fit for its parameters; (4) marginalize over parameterizations as well as parameters. Method (1) establishes a reference minimum RMSE for this test case. Method (2) is a conservative assumption that can be made when the family of possible parameterizations is nested—a polynomial of order $n$ with leading coefficient 0 is a polynomial of order $n-1$. Methods (3) and (4) are different ways of automatically accounting for the different possible parameterizations, in one case (3) by selecting a parameterization and in the other (4) by averaging over the possible parameterizations. The likelihood function that we maximize when using method (4) is the weighted sum of the continuum-marginalized likelihoods of the three continuum models. The weights in this sum are the prior probabilities of each of the models; we assume all three are equally likely.

The RMSEs of the four methods are shown in Figure 2. We compare the methods using three criteria: their robustness to decreasing the SNR and increasing the number of parameters relative to the number of observations; the value of their RMSE at fixed SNR; and the SNR they require to obtain the same RMSE. We consider a method to be robust if its RMSE scales consistently with SNR and true continuum order rather than rapidly increasing at some critical value. The reference, conservative, and parameterization-marginalization methods are robust for SNRs between 10 and 100 while the parameterization-selection method is not. Its RMSE blows up below an SNR of 10 for spectra with 1st order continua and at all SNRs considered for 2nd order continua. While parameterization selection performs as well as the other methods on high SNR spectra, its lack of robustness means that it is not as generally applicable as the other methods.

In terms of RMSE, the parameterization-marginalized estimator is nearly as good as the reference estimator. The ratio RMSEm/RMSEr of the parameterization-marginalized RMSE, RMSEm, to the reference RMSE, RMSEr, is 1, 1.04, and 1.1 for spectra with 0th, 1st, and 2nd order continua. For spectra with 0th or 1st order continua, the conservative estimator is significantly worse than the reference estimator. RMSEc/RMSEr is 1.5, 1.5, and 1, respectively, where RMSEc is the RMSE of the conservative estimator. These ratios are approximately constant across the entire considered SNR range. When analyzing an already acquired set of observations in which a variety of continuum parameterizations are present, using the parameterization-marginalization estimator rather than the conservative estimator will, on average, yield higher accuracy and precision.

It is also useful to compare the SNRs the conservative and parameterization-marginalization estimators require to achieve the same RMSE. For spectra with 0th, 1st, and 2nd continua, the ratio of the required SNRs SNRc/SNRm is 1.6, 1.5, and 0.9. These ratios are, again, consistent across the entire considered SNR range. Assuming that SNR is proportional to the square root of observing time, as is the case for Poisson noise-limited data, these SNR ratios can be converted to required observing time ratios. Reaching the RMSE of the parameterization-marginalized estimator with the conservative estimator takes 1.26, 1.22, and 0.95 times as much observing time. When designing an observing strategy to meet a column density RMSE requirement, using the parameterization-marginalized estimator rather than the conservative estimator can save observing time given a fixed sample or increase the size of a sample given a fixed amount of observing time.

3.3 MCMC efficiency

In ISM absorption spectra, it is common to have multiple lines in a spectrum with shared parameters. These lines can be from the same species, e.g. the Lyman series, or from different species, e.g. from Mg I, Zn II, and Cr II which have overlapping lines in the near ultraviolet. When these lines are in different parts of a spectrum, each part needs its own continuum parameters. This is a case in which analytic marginalization can potentially be more efficient than MCMC marginalization.

We compare how quickly MCMC done using each of the two methods converges and how efficient MCMC done using each method is post-convergence. Which comparison is more informative for choosing a method to use will depend on the purpose of the MCMC run. If the goal of an MCMC run is to estimate some value at low-to-moderate precision, the rate of convergence will be the more important factor. If the goal is instead to estimate some value at high precision, the burn-in period will usually be a small fraction of the total chain and post-convergency efficiency will be more important.

We consider a case where there are $N$ absorption lines with shared central velocities and widths and independent column densities. Each absorption line is in a different spectral region. The continuum in each spectral region is a polynomial of order $M$. The marginalized likelihood has $2+N$ absorption line parameters. The unmarginalized likelihood has $2+N$ absorption line parameters and $N\times M$ continuum parameters. We use the emcee implementation of the Goodman and Weare affine-invariant MCMC ensemble sampler to generate draws from the posterior corresponding to each of these likelihoods. We use the minimum number of “walkers,” which is twice the number of parameters.

We use the Rubin-Gelman statistic $\hat{R}$ (Gelman & Rubin, 1992) to assess convergence. The Rubin-Gelman statistic compares the variance between and within different MCMC instances. If the instances have all converged, these two variances should be approximately equal. We run ten MCMC instances for 12800 (per-walker) steps and compute the Rubin-Gelman statistic from the second half of sub-chains of length $2^{p}\times 100$ for $p=0,1,\ldots,7$. $\hat{R}$ is computed separately for each parameter. Following common usage, we consider convergence to be reached when the $\hat{R}$ of all parameters is less than $1.1$. We run this test for 1, 2, and 3 regions and absorption lines assuming a continuum of order 1. The value of the $\hat{R}$ as a function of (total) number of steps is shown in Figure 3. When there is a single region and line, the MCMC marginalization chain takes twice as many steps as the analytic marginalization chain to converge; when there are two regions, it takes eight times as many steps; when there are three, the MCMC marginalization chain has not converged by the maximum chain length of 12800 while the analytic marginalization chain converges within 1600 steps.

We use the number of independent samples per unit time to assess efficiency. We run MCMC with the marginalized likelihood for 2000 burn-in steps and 8000 converged steps and record the average time per sample, $t_{s}$. Because MCMC with the unmarginalized likelihood takes many steps to converge, we use draws from the converged part of the marginalized likelihood chain as a starting point; these draws only have values for the absorption line parameters. At each set of absorption line parameters, we sample a set of continuum parameters from the conditional distribution discussed in Section 2.1. From this starting point, we run MCMC with the unmarginalized likelihood for 4000 burn-in steps and 36000 converged steps and record the average time per sample. We then compute the average integrated autocorrelation times $\tau_{f}$ of the walkers in both chains. The number of independent samples per unit time is $n_{i}=\left(\tau_{f}\,t_{s}\right)^{-1}$.

We compute $n_{i}$ for a number of regions $N=1,2,\ldots,6$, continua of polynomial order $M=$ 1, 2, and 3, and either a trivial LSF or a banded LSF. The ratio $n_{ind}^{(\text{m})}/n_{ind}^{(\text{u})}$ for each of these cases is shown in Figure 4. When this ratio is greater than 1, running MCMC with the marginalized likelihood for a fixed amount of time will produce more independent samples than running MCMC with the unmarginalized likelihood for the same amount of time. The greater the number of regions and the order of the continuum, the greater the efficiency advantage of the marginalized likelihood over the unmarginalized likelihood. This advantage will not depend on the number of datapoints in each spectral region so long as the LSF is trivial or banded, since in these cases the evaluation time of both likelihoods grows linearly with dataset length (see Section A.3).

4 Discussion

4.1 Assumptions and consequences

The explicit assumptions of the analytic marginalization method are that the continuum is a linear function, that the prior on the coefficients is the improper uniform or multivariate normal distribution, that residuals between the data and model are normally distributed, and that the covariance matrix of the residuals does not depend on the continuum. It is obvious that these assumptions do not hold strictly for any dataset. For example, both possible priors require that there not be constraints on the coefficients even though no background source produces negative flux. A less trivial example is data in the low photon count regime, which are better described by a Poisson distribution than a Gaussian distribution. This is particularly important when the uncertainties on the measurements are themselves highly uncertain and should be explicitly modeled. In that case, the uncertainties will depend on the Poisson intensity function, which explicitly depends on the continuum. Analytic marginalization of the kind described in this work should not be applied to low SNR X-ray or UV spectra.

An implicit assumption of the method is that the absorption model is realistic. For analytic marginalization to be useful, it must be possible for the absorption model to correctly describe the actually present absorption features. For example, if a region of a spectrum contains two clearly distinct absorption lines but the model only allows for a single line, the presence of the un-modeled line will bias the continuum model. In short, improvements in continuum modeling cannot solve problems of absorption model misspecification.

The continuum models envisaged in this work will usually be effective descriptions rather than (often non-linear) physical descriptions. Most continua that vary over longer wavelength scales than the width of absorption lines in question can be approximated in this way. Examples of background sources with slowly varying continua include quasars and (particularly rapidly rotating) hot stars. With flexible linear models such as splines, it is even possible to describe more complicated pseudo-continua such as stellar wind lines. For even more complicated pseudo-continua such as those of cool stars (Zasowski et al., 2015, e.g.), it is necessary to use a non-linear model. Marginalizable linear models can still be useful even in this case as a way of introducing small corrections for pseudo-continuum features that are not perfectly described by the non-linear model.

4.2 Applications of analytic marginalization

The test cases in Section 3 showed that marginalization over continuum parameters and parameterizations is more precise, accurate, and robust than the alternatives. Considered purely as a replacement for numerical marginalization, analytic marginalization is just a potentially more computationally efficient way of implementing an existing inference approach. However, it also allows two qualitatively new approaches: continuum model averaging and absorption parameter optimization with a continuum-marginalized likelihood function.

The test case in Section 3.2 combines both of these approaches—optimizing an absorption parameter likelihood function where the parameterization and parameters of the continuum have been marginalized over. Analytic marginalization makes this possible in two ways: availability of closed form likelihoods and availability of gradients of closed form likelihoods. Optimization with continuum-marginalized likelihoods is useful for analyzing large surveys. Analyses of absorption lines in tens of thousands of spectra (Zhu & Ménard, 2013; Zasowski et al., 2015, e.g.) cannot practically be done with MCMC. With analytic marginalization, it is possible to at least marginalize over continuum parameters. The results of the test cases suggest that this approach could mean a non-trivial improvement in the accuracy and precision of absorption line measurements.

In cases where MCMC is possible, combining continuum parameterization marginalization with a probabilistic specification of absorption component structure would allow absorption line analysis with human intervention only at the level of specifying priors and candidate continuum parameterizations. Component structure specification can be done using trans-dimensional inference, in which the dimensionality of parameter space (in this case the number of sets of absorption line parameters) is itself a parameter of the model. This way of doing absorption line analysis has two potential advantages. First, marginalizing over the velocity structure of the absorption as well as the continuum should automatically include effects such as unresolved saturated structure in parameter uncertainties. Second, because inference approach is almost completely automatic, it allows blinding, which improves reproducibility.

5 Conclusion

Absorption lines are an important source of information about stars and the ISM. As larger spectroscopic datasets become available and as reproducibility becomes more standard in astronomy, it becomes necessary to move beyond ad-hoc analysis methods, particularly ones in which a human directly interacts with data. In multiple recent works, there have been attempts to partially automate continuum placement by including and marginalizing over continuum parameters in probabilistic spectral models. Marginalizing over continuum parameters has, in these works, been hypothesized to also improve the accuracy of the recovered absorption line parameters. Despite these advantages, this approach has so far not become popular, in part due to the computational expense of numerically marginalizing over these additional parameters.

In this work, we have shown that it is possible in many cases to replace this numerical marginalization with analytic marginalization (Section 2). Analytic marginalization speeds up MCMC-based analyses in problems with many continuum parameters (Section 3.3). The continuum parameter-marginalized likelihood can also be used for optimization over absorption line parameters. This approach combines the speed of optimization with the advantages of continuum marginalization. Analytic marginalization over continuum parameters makes it trivial to also marginalize over continuum parameterizations. As with parameter marginalizaton, parameterization marginalization can be combined with optimization over absorption line parameters. parameterization marginalization further reduces the amount of direct human interference in the analysis of individual spectra and will be especially useful in analyses of datasets containing spectra with different continuum shapes.

We have also confirmed that marginalization over continuum parameters and parameterizations indeed improves the accuracy of absorption line parameter measurements. The advantage of parameter marignalization is only significant at low SNRs (Section 3.1). On the other hand, parameterization marginalization is significantly more accurate than alternative methods of deciding on a continuum parameterization at all SNRs (Section 3.2).

We have released an open-source python package, amlc, which can be used to evaluate continuum parameter-marginalized likelihoods and related quantities. Features of this package are described in Appendix A. It is meant to be used as a drop-in replacement for likelihood functions in existing absorption spectrum analysis tools.

The author thanks Andrew Casey, Andrew Fox, Cameron Liang, and Yong Zheng for useful discussions about use cases and Joshua Peek and Linda Tchernyshyov for helpful comments.

Appendix A Implementation and demonstration

In this Appendix, we describe how amlc is implemented (Section A.1), list some of its capabilities (Section A.2), and show how the computation time of different calculations grows with dataset and continuum model size (Section A.3).

A.1 Implementation

We have implemented amlc as a pure-Python package with numpy and scipy as dependencies. amlc does not contain functionality for building LSFs or computing transmittances from absorption parameters and is not intended to be a stand-alone analysis tool. It is meant to be used as a drop-in likelihood function replacement in analysis packages or scripts.

A.2 Package functionality

This package was designed for a use case where the log marginal likelihood and its gradient are evaluated at many different values of the $\theta$-dependent parameters (see Section 2) while the $\theta$-independent parameters are held constant. The core feature of the package is the MarginalizedLikelihood class. A MarginalizedLikelihood instance stores $\theta$-independent parts of the model and pre-computes quantities that are re-used during repeated marginalized likelihood evaluations. In particular, it stores the data covariance matrix ${\bf{K}}$; the ${\bf{c}}$ prior covariance matrix ${\bf{\Lambda}}$ and its explicit inverse, if applicable; and the LSF mapping ${\bf{L}}$ and its transpose.

Both covariance matrices can be diagonal or fully general. The package includes the CovarianceMatrix class, which defines a consistent interface for calculations, and two subclasses, DiagonalCovarianceMatrix and GeneralCovarianceMatrix. DiagonalCovarianceMatrix wraps the simple, one-dimensional determinant and inverse calculations possible with a covariance matrix consisting purely of variances and does the book-keeping required to produce output with the correct shape. GeneralCovarianceMatrix uses the Cholesky decomposition of the supplied covariance matrix to calculate its determinant and to left multipy matrices and vectors by its inverse. Computing the Cholesky decomposition of a general covariance matrix of size $M$ by $M$ takes $\mathcal{O}(M^{3})$ calculations, making it prohibitively computationally expensive for large $M$.

The LSF mapping ${\bf{L}}$ can be any object that implements the matrix multiplication interface, i.e. has a matmul or matmul method. For example, ${\bf{L}}$ can be a dense matrix represented by a numpy array, a sparse matrix represented by a scipy.sparse matrix, or a convolution operator represented by a scipy.sparse.linalg LinearOperator. ${\bf{L}}$ can also be the identity mapping (indicated by None), in which case it is left out of any likelihood calculations.

A.3 Computation time as a function of dataset and basis size

The most time-consuming step in computing all of the quantities derived in Section 2 is forming the matrix ${\bf{C}}_{n/u}$. This step requires matrix-matrix products, while most other steps only involve matrix-vector products. These expensive products are ${\bf{L}}{\bf{B}}$ and ${\bf{K}}^{-1}\left({\bf{L}}{\bf{B}}\right)$. The amount of time required to compute these products depends on the structure ${\bf{L}}$ and ${\bf{K}}$.

${\bf{L}}$ can be the identity matrix, a dense matrix, a sparse matrix, or a linear mapping such as convolution. The fastest case is when ${\bf{L}}$ is the identity matrix, since then ${\bf{L}}{\bf{B}}$ does not need to be computed. The slowest case is when it is a dense matrix, in which case computation time grows as $\mathcal{O}(MN(P+Q))$. When ${\bf{L}}$ is a sparse matrix or linear mapping, the scaling depends on its exact structure. An LSF that varies with wavelength can be represented by a banded matrix, which will be sparse if the spectrum spans many resolution elements. If the bandwidth of ${\bf{L}}$ is independent of the size of the dataset, the computation time of this product grows as $\mathcal{O}(M(P+Q))$.

We consider covariance matrices ${\bf{K}}$ that are either diagonal or general. If ${\bf{K}}$ is diagonal, ${\bf{K}}^{-1}\left({\bf{L}}{\bf{B}}\right)$ requires exactly $M(P+Q)$ multiplications. When ${\bf{K}}$ is a general covariance matrix, we decompose it into its Cholesky factors and left-multiply ${\bf{L}}{\bf{B}}$ by ${\bf{K}}^{-1}$ by solving the linear problem ${\bf{{\bf{L}}{\bf{B}}}}={\bf{K}}{\bf{X}}$. The time needed to factor ${\bf{K}}$ grows as $\mathcal{O}\left(M^{3}\right)$ but only needs to be done once per set of observations. The time needed to solve the linear problem grows as $\mathcal{O}\left(M^{2}(P+Q)\right)$.

To empirically confirm these growth rates, we timed how long it takes to evaluate the log-likelihood and its gradient for a range of dataset sizes $M$ and basis sizes $P+Q$ and three ${\bf{L}}$ and ${\bf{K}}$ structure scenarios. The scenarios are: ${\bf{L}}$ is the identity mapping, ${\bf{K}}$ is diagonal; ${\bf{L}}$ is a dense matrix, ${\bf{K}}$ is general; and ${\bf{L}}$ is a sparse, banded matrix and ${\bf{K}}$ is diagonal. The first two scenarios are the fastest and slowest combination. The third scenario is more typical for a spectrum; the data uncertainty is diagonal and the LSF has finite extent. The evaluation time of the log-likelihood as a function of $M$ and $P+Q$ for these three scenarios is shown in Figures 5, 6, and 7. We do not show the evaluation time of the gradient because it behaves in the same way as the evaluation time of the log-likelihood in all three scenarios; the most expensive step of the two calculations is the same.

The dependence of computation time on $M$ and $P+Q$ generally agrees with the predictions based on the two most time-consuming steps. At low $M$ and in particular at low $P+Q$, the computation time is either overhead-dominated or evenly split between the most time-consuming steps and other steps. When $M\gtrsim 10^{5}$, computation time increases faster than expected purely from the growth rate of the required number of operations (see e.g. the left panel of Figure 5). This excess increase in computation time is most likely due to changes in memory bandwidth, as the size of matrix rows and columns increases past the size of the highest-level CPU cache on the laptop used to run these tests.

To put these dataset sizes into context, a Sloan Digital Sky Survey (SDSS) BOSS or APOGEE spectrum is $\sim 10^{3}$ pixels long, a Hubble Space Telescope Cosmic Origins Spectrograph (HST-COS) spectrum is $\sim 10^{4}$ pixels long, and a spectrum from an echelle spectrograph such as the Ultraviolet and Visual Echelle Spectrograph on the Very Large Telescope or the Magellan Inamori Kyocera Echelle spectrograph is $\sim 10^{5}-10^{6}$ pixels long. The uncertainties associated with these spectra are usually assumed to be diagonal and the LSFs are acceptably described by sparse, banded matrices, so the computation times given in Figure 5 and 7 should apply.