Spline Baselines

Introduction

The two dimensional extension of penalized splines (P-splines) for baseline correction within pybaselines follows the framework of Eilers, Currie, and Durbán from [1].

Let the number of rows be \(M\) and the number of columns \(N\) within the matrix of measured data \(Y\). Note that \(y\) is the flattened array of matrix \(Y\) with length \(M * N\). Let \(Y\) be a function of \(x\) along the rows and \(z\) along the columns, ie. \(Y_{ij} = f(x_i, z_j)\), and \(B_r(x)\) and \(B_c(z)\) represent the spline basis matrices along the rows and columns, respectively, each with a number of knots \(g\) and h. Analogous to the 1D case, the goal is to make the baseline, \(V\) match the measured data as well as it can while also penalizing the difference between spline coefficients, resulting in the following minimization:

\[\sum\limits_{i}^M \sum\limits_{j}^N W_{ij} (Y_{ij} - \sum\limits_{k}^g \sum\limits_{l}^h B_{r,k}(x_i) B_{c,l}(z_j) \alpha_{kl})^2 + \lambda_r \sum\limits_{i}^{g - d_r} (\alpha_{i\bullet} \Delta^{d_r})^2 + \lambda_c \sum\limits_{j}^{h - d_c} (\Delta^{d_c} \alpha_{j\bullet})^2\]

and

\[V = \sum\limits_{i}^g \sum\limits_{j}^h B_{r,i} B_{c,j} \alpha_{ij}\]

where \(Y_{ij}\) is the measured data, \(\alpha\) is the matrix of spline coefficients, \(\lambda_r\) is the penalty along the rows, \(\lambda_c\) is the penalty along the columns, \(W_{ij}\) is the weighting, \(\Delta^{d_r}\) is the finite-difference operator of order \(d_r\) along each row of \(\alpha\), \(\alpha_{i\bullet}\), and \(\Delta^{d_c}\) is the finite-difference operator of order \(d_c\) along each column of \(\alpha\), \(\alpha_{j\bullet}\).

Let \(B = B_c \otimes B_r\) denote the kronecker product of the basis matrices for the columns and rows, which represents the overall two dimensional tensor product spline basis. The resulting linear equation for solving the above minimization is:

\[(B^{\top} W_{diag} B + \lambda_r I_h \otimes D_{d_r}^{\top} D_{d_r} + \lambda_c D_{d_c}^{\top} D_{d_c} \otimes I_g) \alpha = B^{\top} W_{diag} y\]

and the baseline is then:

\[v = B \alpha\]

where \(W_{diag}\) is the diagaonal matrix of the flattened weights, \(v\) is the flattened estimated baseline, and \(D_d\) is the matrix version of \(\Delta^d\), as already explained for the 1D case. Further, \(\otimes\) denotes the Kronecker product, and \(I_g\) and \(I_h\) are the identity matrices of length \(g\) and \(h\), respectively. After solving, the array \(v\) can then be reshaped into the matrix \(V\).

Since experimental data is measured on gridded data (ie. \(Y_{ij} = f(x_i, z_j)\)), the above equation can be optimized following [1] and expressed as a generalized linear array model which allows directly using the matrices of the measured data, \(Y\), and the weights, \(W\), rather than flattening them, which significantly reduces the required memory and computation time.

Let \(F\) be the face-splitting product operator of a matrix with itself such that \(F(B_r) = (B_r \otimes 1_{g}^{\top}) \odot (1_{g}^{\top} \otimes B_r)\) and \(F(B_c) = (B_c \otimes 1_{h}^{\top}) \odot (1_{h}^{\top} \otimes B_c)\), where \(1_g\) and \(1_h\) are vectors of ones of length \(g\) and \(h\), respecitvely, and \(\odot\) signifies elementwise multiplication. Then the linear equation can be rewritten as:

\[(F(B_r)^{\top} W F(B_c) + \lambda_r I_h \otimes D_{d_r}^{\top} D_{d_r} + \lambda_c D_{d_c}^{\top} D_{d_c} \otimes I_g) \alpha = B_{r}^{\top} (W \odot Y) B_c\]

and the baseline is:

\[V = B_r \alpha B_{c}^{\top}\]

Algorithms

Spline Baselines

Introduction

Algorithms

mixture_model (Mixture Model)

irsqr (Iterative Reweighted Spline Quantile Regression)

pspline_asls (Penalized Spline Version of asls)

pspline_iasls (Penalized Spline Version of iasls)

pspline_airpls (Penalized Spline Version of airpls)

pspline_arpls (Penalized Spline Version of arpls)

pspline_iarpls (Penalized Spline Version of iarpls)

pspline_psalsa (Penalized Spline Version of psalsa)