We propose an ultrasparse graphical model for representing multi-indexed data based on a Kronecker sum representation of the process inverse covariance matrix. This statistical model decomposes the inverse covariance into a linear Kronecker sum representation with sparse Kronecker factors. Under the assumption that the observations are matrix-normal the l1 sparsity regularized log-likelihood function is convex and admits significantly faster statistical rates of convergence than other sparse matrix normal algorithms such as graphical lasso or Kronecker graphical lasso. We establish pointwise convergence of an iterative gradient ascent algorithm to a global maximum and obtain high dimensional statistical convergence rates. We will illustrate TeraLasso on real multi-indexed spatio-temporal data.