  We consider the problem of learning a sparse graph underlying an undirected
Gaussian graphical model, a key problem in statistical machine learning. Given
$n$ samples from a multivariate Gaussian distribution with $p$ variables, the
goal is to estimate the $p \times p$ inverse covariance matrix (aka precision
matrix), assuming it is sparse (i.e., has a few nonzero entries). We propose
GraphL0BnB, a new estimator based on an $\ell_0$-penalized version of the
pseudolikelihood function, while most earlier approaches are based on the
$\ell_1$-relaxation. Our estimator can be formulated as a convex mixed integer
program (MIP) which can be difficult to compute at scale using off-the-shelf
commercial solvers. To solve the MIP, we propose a custom nonlinear
branch-and-bound (BnB) framework that solves node relaxations with tailored
first-order methods. As a by-product of our BnB framework, we propose
large-scale solvers for obtaining good primal solutions that are of independent
interest. We derive novel statistical guarantees (estimation and variable
selection) for our estimator and discuss how our approach improves upon
existing estimators. Our numerical experiments on real/synthetic datasets
suggest that our method can solve, to near-optimality, problem instances with
$p = 10^4$ -- corresponding to a symmetric matrix of size $p \times p$ with
$p^2/2$ binary variables. We demonstrate the usefulness of GraphL0BnB versus
various state-of-the-art approaches on a range of datasets.
