  In many modern regression applications, the response consists of multiple
categorical random variables whose probability mass is a function of a common
set of predictors. In this article, we propose a new method for modeling such a
probability mass function in settings where the number of response variables,
the number of categories per response, and the dimension of the predictor are
large. Our method relies on a functional probability tensor decomposition: a
decomposition of a tensor-valued function such that its range is a restricted
set of low-rank probability tensors. This decomposition is motivated by the
connection between the conditional independence of responses, or lack thereof,
and their probability tensor rank. We show that the model implied by such a
low-rank functional probability tensor decomposition can be interpreted in
terms of a mixture of regressions and can thus be fit using maximum likelihood.
We derive an efficient and scalable penalized expectation maximization
algorithm to fit this model and examine its statistical properties. We
demonstrate the encouraging performance of our method through both simulation
studies and an application to modeling the functional classes of genes.
