\indent models and estimators for covariance matrices are very well studied. for non-gaussian distributions, simply studying covariance gives an incomplete picture. extending the edgeworth series gives the pxpxp skewness tensor, the pxpxpxp kurtosis tensor, and so on. we describe a strategy for building multilinear factor models of cumulant tensors using subspace varieties. this leads to a difficult optimization problem and a fully implicit, gradient-based numerical optimization method using parallel transport on the grassmannian to perform estimation. we also discuss some of the associated statistical challenges and applications.