  In the current era of vast data and transparent machine learning, it is
essential for techniques to operate at a large scale while providing a clear
mathematical comprehension of the internal workings of the method. Although
there already exist interpretable semi-parametric regression methods for
large-scale applications that take into account non-linearity in the data, the
complexity of the models is still often limited. One of the main challenges is
the absence of interactions in these models, which are left out for the sake of
better interpretability but also due to impractical computational costs. To
overcome this limitation, we propose a new approach using a factorization
method to derive a highly scalable higher-order tensor product spline model.
Our method allows for the incorporation of all (higher-order) interactions of
non-linear feature effects while having computational costs proportional to a
model without interactions. We further develop a meaningful penalization scheme
and examine the induced optimization problem. We conclude by evaluating the
predictive and estimation performance of our method.
