# Numerical rank deficiency in P-spline based varying coefficient model

I am studying the following varying coefficient model from [1]:

$$d_x(x,y) = s_x(y)x + c_x(y) + o_x(x)$$ $$d_y(x,y) = s_y(y)x + c_y(y) + o_y(x)$$

where $$x$$ and $$y$$ are covariates, $$d_x$$ and $$d_y$$ are data coming from measurements and $$s_{x,y}(y)$$,$$c_{x,y}(y)$$,$$o_{x,y}(x)$$, are smooth functions modeled by P-splines [2]. The main goal is to recover the coefficients of these P-splines in order to use the model to make predictions. These model equations can be used to assemble the cost function $$c_x(\mathbf{p}) = ||\mathbf{d}-A \mathbf{p}||^2+\lambda_s ||D_s \mathbf{p_s}||^2 +\lambda_c ||D_c \mathbf{p_c}||^2+\lambda_o ||D_o \mathbf{p_o}||^2$$ where $$\mathbf{p}=\left[\mathbf{p_s}\, \mathbf{p_c}\, \mathbf{p_o} \right]$$ is the vector collecting all the spline coefficients to be estimated, $$A$$ is a matrix whose blocks have elements computed by evaluating the B-spline bases in the data points abscissas $$\{ x_i \}$$ and ordinates $$\{ y_i \}$$, $$\mathbf{d}$$ is the vector of measurements, $$D_s$$, $$D_c$$ and $$D_o$$ are the second order divided difference operators used to enforce smoothness of $$\mathbf{p}$$ and the $$\lambda$$'s are tuning parameters. For more details please see [2]. Given the linearity of the cost gradient, we can estimate $$\mathbf{p}$$ by solving the linear least squares problem: $$\hat{\mathbf{p}} = (A^TA + P)^{-1} A^T \mathbf{d}$$ $$P = \text{diag}(\lambda_s D_s^T D_s,\lambda_c D_c^T D_c,\lambda_o D_o^T D_o)$$

When I solve this linear system using synthetic data for $$\mathbf{d}$$ I see that the matrix to be inverted as well as the $$A$$ matrix are numerically rank deficient. To cope with this I use a minimum norm least square solver. The authors of [2] suggest in some of their papers to redefine the $$P$$ matrix as $$P = \text{diag}(\lambda_s D_s^T D_s,\lambda_c D_c^T D_c,\lambda_o D_o^T D_o) + \lambda I$$ with $$\lambda=10^{-6}$$ to reach the same goal.

Now my question is: is the ill-conditioning of this linear system something generally expected in these models? Are there more principled ways to address it?

[1] Saquib, Suhail S.et al. "Spline warp model for registering pushbroom multispectral imagery." Long-Range Imaging II. Vol. 10204. SPIE, 2017.

[2] Eilers, Paul HC, and Brian D. Marx. Practical smoothing: The joys of P-splines. Cambridge University Press, 2021.

• The penalties are what should make the system full rank. From your question, I don't see what $P$ is supposed to be besides just $P = \lambda I$. If that's the case, then this problem is identifiable when $\lambda > 0$. Feb 16, 2023 at 4:54
• Sorry if this was not clear, in the most general case $P = \text{diag}(\lambda_s D_s^T D_s,\lambda_c D_c^T D_c,\lambda_o D_o^T D_o) + \lambda I$ with $\lambda >= 0$. Feb 16, 2023 at 16:47

As long as $$P$$ is a diagonal matrix with positive entries on the diagonal, then $$(A^tA + P)$$ should be invertible and the solution should be unique.
• This is correct, however as I point out the problem is not inversion of $(A^{\rm T}A+P)$ but its very high condition number. I see condition numbers of the order of $10^{20}$. Now I suspect that this is caused by an identifiability problem. Feb 17, 2023 at 19:31
• That is what I was thinking too, however increasing the $\lambda$s did not help. I believe that the problem is identifiability of the model. You can add $\delta$ to $c_x(y)$ and subtract it from $o_x(x)$, for example, and the model will still fit the data with a low error. In my specific problem I see actually two unidentified parameters since there are two singular values close to machine precision. I was able to eliminate one using the reparameterization method presented in Wood's book on GAMs, but there is still a parameter that I cannot identify. Feb 21, 2023 at 20:04