  Quadratic regression involves modeling the response as a (generalized) linear
function of not only the features $x^{j_1}$ but also of quadratic terms
$x^{j_1}x^{j_2}$. The inclusion of such higher-order "interaction terms" in
regression often provides an easy way to increase accuracy in
already-high-dimensional problems. However, this explodes the problem dimension
from linear $O(p)$ to quadratic $O(p^2)$, and it is common to look for sparse
interactions (typically via heuristics). In this paper, we provide a new
algorithm - Interaction Hard Thresholding (IntHT) which is the first one to
provably accurately solve this problem in sub-quadratic time and space. It is a
variant of Iterative Hard Thresholding; one that uses the special quadratic
structure to devise a new way to (approx.) extract the top elements of a $p^2$
size gradient in sub-$p^2$ time and space. Our main result is to theoretically
prove that, in spite of the many speedup-related approximations, IntHT linearly
converges to a consistent estimate under standard high-dimensional sparse
recovery assumptions. We also demonstrate its value via synthetic experiments.
Moreover, we numerically show that IntHT can be extended to higher-order
regression problems, and also theoretically analyze an SVRG variant of IntHT.
