Convex bivariate B-spline by construction (PSD Hessian)?

For a univariate spline f(u), I can enforce convexity “by construction” by parametrizing the second derivative spline with positive coefficients:

alpha = exp(theta);               % alpha  >= 0
fdd   = spmak(knots, alpha);      % f''(u) >= 0
% fnint(fnint, fdd) --> f(u)

This parameterization enables unconstrained optimization (e.g. lsqnonlin) in theta-space.

Now I want a similar idea for a bivariate spline f(u,v): can convexity similarly be enforced “by construction”, i.e. the Hessian is PSD on the spline domain, without inequality constraints?

In 2D the mixed term makes it hard (need something like f_uu*f_vv - f_uv^2 >= 0). Are there spline parametrizations / bases that guarantee PSD Hessian?

Any ides are appreciated, thank you!

Answer 1

[John D'Errico]

Even without thinking about it, I'm pretty sure the answer is no. Things get nasty when you move beyond one dimension, and 2d is far enough to cause serious problems. The basic example I'll point out is that of monotonicity.

Consider the simplest shape you can possibly describe, a piecewise linear function in both x and y, defined on the unit square [0,1]X[0,1].

fxy = [1 2;2 5];
syms x y
fun(x,y) = (1-x).*(1-y)*fxy(1,1) + (1-x).*y*fxy(1,2) + x.*(1-y)*fxy(2,1) + x.*y*fxy(2,2);

This form of fun is often called a bilinear interpolant, or sometimes a tensor product linear interpolant. It is easily shown to be monotonic in both x and y, so if you hold x or y fixed at any constant, then the function is linear in the other variable.

fsurf(fun,[0,1,0,1])

Looking at the figure, you can see the shape is linear along any edge. It is also purely linear if you fix either x or y. For example, at y==0.5, we see nicely linear behavior in x.

simplify(fun(x,0.5))

ans = 

But as you can see if we look at the function along some paths through the unit square that are not parallel to an axis, the result need not be so well-behaved. Here, for example, we see that along the anti-diagonal, with y = 1-x, we lose that illusion of monotonicity.

fplot(fun(x,1-x),[0,1])

Again, the point is, what you can get away with in 1-d often fails when you try to do the same in 2-d.

[SA-W]

That example makes the geometric difficulty very clear. Always a good idea to recall "simple" cases.

If you had to fit a smooth bivariate spline surface (bicubic) and impose “Hessian PSD” during the nonlinear fit as cleanly as possible in MATLAB, how would you do it in practice (parameterization / costraints)?