We consider the least squares regression problem, penalized with a\ncombination of the $\\ell_{0}$ and squared $\\ell_{2}$ penalty functions (a.k.a.\n$\\ell_0 \\ell_2$ regularization). Recent work shows that the resulting\nestimators are of key importance in many high-dimensional statistical settings.\nHowever, exact computation of these estimators remains a major challenge.\nIndeed, modern exact methods, based on mixed integer programming (MIP), face\ndifficulties when the number of features $p \\sim 10^4$. In this work, we\npresent a new exact MIP framework for $\\ell_0\\ell_2$-regularized regression\nthat can scale to $p \\sim 10^7$, achieving speedups of at least $5000$x,\ncompared to state-of-the-art exact methods. Unlike recent work, which relies on\nmodern commercial MIP solvers, we design a specialized nonlinear\nbranch-and-bound (BnB) framework, by critically exploiting the problem\nstructure. A key distinguishing component in our framework lies in efficiently\nsolving the node relaxations using a specialized first-order method, based on\ncoordinate descent (CD). Our CD-based method effectively leverages information\nacross the BnB nodes, through using warm starts, active sets, and gradient\nscreening. In addition, we design a novel method for obtaining dual bounds from\nprimal CD solutions, which certifiably works in high dimensions. Experiments on\nsynthetic and real high-dimensional datasets demonstrate that our framework is\nnot only significantly faster than the state of the art, but can also deliver\ncertifiably optimal solutions to statistically challenging instances that\ncannot be handled with existing methods. We open source the implementation\nthrough our toolkit L0BnB.\n