R/mgcv: Why do te() and ti() tensor products produce different surfaces?

The mgcv package for R has two functions for fitting tensor product interactions: te() and ti(). I understand the basic division of labour between the two (fitting a non-linear interaction vs. decomposing this interaction into main effects and an interaction). What I don't understand is why te(x1, x2) and ti(x1) + ti(x2) + ti(x1, x2) may produce (slightly) different results.

MWE (adapted from ?ti):

require(mgcv)
test1 <- function(x,z,sx=0.3,sz=0.4) {
x <- x*20
(pi**sx*sz)*(1.2*exp(-(x-0.2)^2/sx^2-(z-0.3)^2/sz^2)+
0.8*exp(-(x-0.7)^2/sx^2-(z-0.8)^2/sz^2))
}
n <- 500

x <- runif(n)/20;z <- runif(n);
xs <- seq(0,1,length=30)/20;zs <- seq(0,1,length=30)
pr <- data.frame(x=rep(xs,30),z=rep(zs,rep(30,30)))
truth <- matrix(test1(pr$x,pr$z),30,30)
f <- test1(x,z)
y <- f + rnorm(n)*0.2

par(mfrow = c(2,2))

# Model with te()
b2 <- gam(y~te(x,z))
vis.gam(b2, plot.type = "contour", color = "terrain", main = "tensor product")

# Model with ti(a) + ti(b) + ti(a,b)
b3 <- gam(y~ ti(x) + ti(z) + ti(x,z))
vis.gam(b3, plot.type = "contour", color = "terrain", main = "tensor anova")

# Scatterplot of prediction b2/b3
plot(predict(b2), predict(b3))


The differences aren't very large in this example, but I'm just wondering why there should be differences at all.

Session info:

 > devtools::session_info("mgcv")
Session info
-----------------------------------------------------------------------------------
setting  value
version  R version 3.3.1 (2016-06-21)
system   x86_64, linux-gnu
ui       RStudio (0.99.491)
language en_US
collate  en_US.UTF-8
tz       <NA>
date     2016-09-13

Packages      ---------------------------------------------------------------------------------------
package * version date       source
lattice   0.20-33 2015-07-14 CRAN (R 3.2.1)
Matrix    1.2-6   2016-05-02 CRAN (R 3.3.0)
mgcv    * 1.8-12  2016-03-03 CRAN (R 3.2.3)
nlme    * 3.1-128 2016-05-10 CRAN (R 3.3.1)

• Seriously people!? Whilst the implementation in clearly an mgcv-specific thing (I'm unaware of any other off-the-shelf software for GAMs that allows this ANOVA-like decomposition of bivariate smooths), the problem and answer are clearly statistical; the models being fitting aren't the same under the hood due to the extra penalty matrices that arise when decomposing the marginal terms from the "interaction" component. This is not specific to mgcv. – Gavin Simpson Sep 14 '16 at 15:58
• @GavinSimpson I have raised a question on Meta about the on-topicality, or otherwise, of this question – Silverfish Sep 17 '16 at 13:31

These are superficially the same model but in practice when fitting there are some subtle differences. One important difference is that the model with ti() terms is estimating more smoothness parameters compared with the te() model:
> b2$sp te(x,z)1 te(x,z)2 3.479997 5.884272 > b3$sp

and this is because there are more penalty matrices associated with the two models; in the ti() model we have one per "term" compared with just two in the te() model, one per marginal smooth.
I see models with ti() as being used to decide whether I want $\hat{y} = \beta_0 + s(x, y)$ or $\hat{y} = \beta_0 + s(x) + s(y)$. I can't compare those models if I use te() terms so I use ti(). Once I've determined if I need $s(x,y)$ I can refit the model with te() if I need it or with separate s() for each marginal effect if I don't need $s(x,y)$.
Note that you can get the models somewhat closer to each other by fitting using method = "ML" (or "REML", but you shouldn't be comparing "fixed" effects with "REML" unless all the terms are fully penalized, which by default they aren't, but would be say with select = TRUE).