  We study the problem of minimizing the sum of three convex functions: a
differentiable, twice-differentiable and a non-smooth term in a high
dimensional setting. To this effect we propose and analyze a randomized block
cubic Newton (RBCN) method, which in each iteration builds a model of the
objective function formed as the sum of the natural models of its three
components: a linear model with a quadratic regularizer for the differentiable
term, a quadratic model with a cubic regularizer for the twice differentiable
term, and perfect (proximal) model for the nonsmooth term. Our method in each
iteration minimizes the model over a random subset of blocks of the search
variable. RBCN is the first algorithm with these properties, generalizing
several existing methods, matching the best known bounds in all special cases.
We establish ${\cal O}(1/\epsilon)$, ${\cal O}(1/\sqrt{\epsilon})$ and ${\cal
O}(\log (1/\epsilon))$ rates under different assumptions on the component
functions. Lastly, we show numerically that our method outperforms the
state-of-the-art on a variety of machine learning problems, including cubically
regularized least-squares, logistic regression with constraints, and Poisson
regression.
