  In many learning settings, it is beneficial to augment the main features with
pairwise interactions. Such interaction models can be often enhanced by
performing variable selection under the so-called strong hierarchy constraint:
an interaction is non-zero only if its associated main features are non-zero.
Existing convex optimization based algorithms face difficulties in handling
problems where the number of main features $p \sim 10^3$ (with total number of
features $\sim p^2$). In this paper, we study a convex relaxation which
enforces strong hierarchy and develop a highly scalable algorithm based on
proximal gradient descent. We introduce novel screening rules that allow for
solving the complicated proximal problem in parallel. In addition, we introduce
a specialized active-set strategy with gradient screening for avoiding costly
gradient computations. The framework can handle problems having dense design
matrices, with $p = 50,000$ ($\sim 10^9$ interactions)---instances that are
much larger than current state of the art. Experiments on real and synthetic
data suggest that our toolkit hierScale outperforms the state of the art in
terms of prediction and variable selection and can achieve over a 4900x
speed-up.
