broad interest in incorporating interactions in linear regression. Extensions of linear regression to accommodate pairwise interactions are commonly referred to as quadratic regression. In moderate to high-dimensional settings, it becomes very challenging to implement quadratic regression since the number of parameters to be estimated is $2 p + {p \choose 2}$. Hence, classical methods such as least squares cannot be used and even common penalization and Bayesian methods can encounter computational hurdles. Reliable inferences on main effects and interactions is even more challenging when certain predictors are moderately to highly correlated.

A lot of effort has been focused on estimating pairwise interactions in moderate high-dimensional and ultra high-dimensional problems. We refer to the former when the number of covariates is between $20$ and $100$ and to the latter when $p>100$. When $p = 100$, the number of parameters to be estimated is greater than $5000$. When $p \in [20,100]$, one-stage regularization methods like [@Bien2013] and [@Haris2018] can be successful. Some of these methods require a so-called heredity assumption [@Chipman1996] to reduce dimensionality. Strong heredity means that the interaction between two variables is included in the model only if both main effects are. For weak heredity it suffices to have one 