Tensor product between an ispline and a bspline for fitting data that should be monotonic in one dimension

I’m not very familiar with the process for solving tensor product basis fittings. I’ve done some work with fitting an ispline basis with a non-negative-least-squares solver to fit a monotonic spline to noisy data where the underlying function is guaranteed to be monotonic.

Now I’m interested in taking a tensor product of this ispline basis with a bspline basis. But I’m not certain the best way to do this and still enforce the non-negative constraint. My one thought is to double each of the b-spline basis with its negative. I can take the tensor product of this with the ispline basis, then I should still be able to use a non-negative-least-squares solver. Is this a reasonable approach? Is there a better approach?