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.

  • $\begingroup$ 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$. $\endgroup$
    – Cliff AB
    Feb 16, 2023 at 4:54
  • $\begingroup$ 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$. $\endgroup$ Feb 16, 2023 at 16:47

1 Answer 1


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.

  • $\begingroup$ 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. $\endgroup$ Feb 17, 2023 at 19:31
  • $\begingroup$ @ArrigoBenedetti is it possible that the issue is that there are penalties which are too small then? The larger the penalties, the smaller the conditioning number. $\endgroup$
    – Cliff AB
    Feb 18, 2023 at 2:15
  • $\begingroup$ 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. $\endgroup$ Feb 21, 2023 at 20:04
  • $\begingroup$ Have you removed intercepts from 2 out of 3 of the splines? I believe if you have intercepts in all the spline expansions, you will get non-identifiability. $\endgroup$
    – Cliff AB
    Feb 22, 2023 at 2:29
  • $\begingroup$ P-splines do not have an explicit intercept term, however they can model a generic function with an offset. Usually they are constrained with a method like that presented in Sec. 5.4.1 of Wood, Generalized additive models: an introduction with R. Chapman and Hall/CRC, 2006. $\endgroup$ Feb 22, 2023 at 19:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.