  In the age of big data and interpretable machine learning, approaches need to
work at scale and at the same time allow for a clear mathematical understanding
of the method's inner workings. While there exist inherently interpretable
semi-parametric regression techniques for large-scale applications to account
for non-linearity in the data, their model complexity is still often
restricted. One of the main limitations are missing interactions in these
models, which are not included for the sake of better interpretability, but
also due to untenable computational costs. To address this shortcoming, we
derive a scalable high-order tensor product spline model using a factorization
approach. Our method allows to include all (higher-order) interactions of
non-linear feature effects while having computational costs proportional to a
model without interactions. We prove both theoretically and empirically that
our methods scales notably better than existing approaches, derive meaningful
penalization schemes and also discuss further theoretical aspects. We finally
investigate predictive and estimation performance both with synthetic and real
data.
