Bivariate isotonic regression on unbalanced grid I would like to do classic isotonic regression on bivariate data (that is to say, x with two columns). The function biviso of package Iso is very fast but requires a balanced grid. I would like to do isotonic regression on unbalanced grid. For example:
x1 <- runif(100)
x2 <- runif(100)
y <- x1*x2 + rnorm(100)

Is there an R package or custom R code that can do this? By "classic" isotonic regression, I attempt to make reference to the following definition, which is in the Wikipedia article:
$$
{\displaystyle \min \sum _{i=1}^{n}w_{i}(x_{i}-a_{i})^{2}}
$$
$$
{{\text{subject to }}x_{i}\leq x_{j}{\text{ for all }}(i,j).}
$$
The weights $w_i$ are not important for my usage so we can define $w_i=1$.
 A: The easiest way I think of is using monotonic splines through scam. That way we will constraint the fit of each term to be a monotonic function giving us a monotonic response. I prefer it to other functions that solve the constrained optimisation problems directly because it gives us access to (most of) the diagnostic goodies we get from mgcv::gam.
set.seed(5)
x1 <- runif(100)
x2 <- runif(100)
y <- x1*x2 + rnorm(100)
my_data = data.frame(x1=x1, x2=x2, y=y)
library(scam) 
my_model_sc <- scam(y ~ s(x1, bs="mpi")+ s(x2, bs="mpi"), data=my_data)

A quick visual inspection confirm that our fit is monotonic along x1 and x2 too.
library(pdp)
pd <- partial(my_model_sc, pred.var = c("x1", "x2"), smooth = FALSE,
              grid.resolution = 100, progress = "text" )
plotPartial(pd)


If we are looking to replicate something closer to the step changes seen in QP solutions, we can look into cgam. 
library(cgam) 
my_model_cg <- cgam(y ~ incr(x1)+ incr(x2), data=my_data)

Here we can specify an increasing but smoothed fit through incr (they are other options to have smooth convex or concave fits, etc). 
I append a bit of code below  showing how both functions allow us to have a monotonic fit across values of x2 for a fixed x1. 
my_data_2=my_data[order(my_data$x2),]
my_data_2$x1=0.5;

par(mfrow=c(1,2))
plot(x=my_data_2$x2, y = predict(my_model_cg, newData = my_data_2)$fit, lwd=2,
     t='l', xlab="x2", ylab="f(x1=0.5, x2)", ylim=c(-0.4,0.5), main="CGAM fit")
grid()
plot(x=my_data_2$x2, y = predict(my_model_sc,newdata =my_data_2), lwd=2,  
     t='l', xlab="x2", ylab="f(x1=0.5, x2)", ylim=c(-0.4,0.5), main="SCAM fit")
grid()


