\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.