Papers I have read about p splines in multiple dimensions use some variation of the following idea. The spline is a surface $u(x_1, x_2)$ that that minimizes the objective function $\sum w(u - u')^2 + \lambda_{x_1}\sum(D_{x_1} u)^2 + \lambda_{x_2}\sum(D_{x_2} u)^2$. Where $D_{x_1}$ and $D_{x_2}$ are the difference matrices from Whittaker-Henderson, b splines, etc.. There is a separate smoothing penalty in each dimension. In the 2D case, $\lambda_{x_1}\sum(D_{x_1} u)^2$ is the row penalty and $\lambda_{x_2}\sum(D_{x_2} u)^2$ is the column penalty.
What is the benefit of using separate smoothing penalties for each dimension? Why not define $u$ as the surface that minimizes the objective function $\sum w(u - u')^2 + \sum(\nabla^2 u)^2$, where $\nabla^2 u$ is the discrete Laplacian of $u$? Smoothing weights can be applied by adjusting the definition of the units of each dimension. For example, the discretization will be different if the column dimension measures inches vs feet.
It seems that the p spline objective function is a discretization of $$\sum w(u-u')^2 + \lambda_{x_1}\iint_U\left(\frac{d^m u}{d{x_1}^m}\right)^2 \,d{x_1}\,d{x_2} + \lambda_{x_2}\iint_U\left(\frac{d^n u}{d{x_2}^n}\right)^2 \,d{x_1}\,d{x_2}$$
and the Laplacian version is a discretization of $$\sum w(u-u')^2 + \iint_U\left(\lambda_{x_1}\frac{d^2 u}{d{x_1}^2} + \lambda_{x_2} \frac{d^2 u}{d{x_2}^2}\right)^2 \,d{x_1}\,d{x_2}$$
I see many reasons to prefer the Laplacian version. However, smarter people than me are using the traditional objective function, and I want to know why they prefer it. I think the Laplacian objective function that I described above is making a polyharmonic smoothing spline or maybe a thin plate spline.
