Kronecker product of two base-learners with anisotropic penalty  #28

When %O% is called with a specification of dfin both base-learners,
e.g. bbs(x1, df = df1) %O% bbs(x2, df = df2), the global dffor the
Kroneckered base-learner is computed as df = df1 * df2.
And thus the penalty has only one smoothness parameter lambdaresulting in an isotropic penalty,
K = lambda * ( kronecker( diag, K1) + kronecker(K2, diag ) )
with global penalty K, marginal penalty matrices K1, K2and identity matrices diag.

But I would like to specify a Kroneckered effect with anisotropic penalty. Thus, I implemented a new operator %A%, which allows for a different amount of smoothness in the two directions.
For example bbs(x1, df = df1) %A% bbs(x2, df = df2) results in computing, two
different smoothing parameters lambda1, lambda2 for the two marginal effects and a global lambdato adjust for the global df, i.e.
K = lambda * ( kronecker( diag, lambda1*K1) + kronecker(lambda2*K2, diag ) )
where lambda1 is computed for df1 and lambda2 is computed for df2.
lambdais computed such that the global dfhold again df = df1 * df2.

The current implementation of %A% is experimental and cannot deal with weights in the model call.
See anisotropicKronecker.txt for the current EXPERIMENTAL implementation in R.

For a MWE on how to use %A% instead of %O% (does not make much sense in that case):

library(mboost)
library(Matrix)

## estimate mean of image treating image as matrix
x1 <- 1:nrow(volcano)
x2 <- 1:ncol(volcano)
vol <- as.vector(volcano)

## model with isotropic penalty
mod <- mboost(vol ~ bbs(x1, df = 3, knots = 10) %O% bbs(x2, df = 3, knots = 10),
              control = boost_control(nu = 0.25))
mod[250]

## model with anisotropic penalty
mod_a <- mboost(vol ~ bbs(x1, df = 2, knots = 10) %A% bbs(x2, df = 5, knots = 10),
              control = boost_control(nu = 0.25))
mod_a[250]

par(mfrow=c(1,3))
image(volcano, main = "data")
volf <- matrix(fitted(mod), nrow = nrow(volcano))
image(volf, main = "fitted isotropic")
volf_a <- matrix(fitted(mod_a), nrow = nrow(volcano))
image(volf_a, main = "fitted anisotropic")

In issue #19 I covered the special case where lambda1=0 , which is covered here, if df1 = 'rank of design-matrix for x1'.

Do you think it would be generally helpful to have an anisotropic effect like %A% in mboost? Do you think the current implementation is sensible? Any comments and suggestions are welcome.

Best
Sarah