We propose a novel approach to solve exactly the sparse nonnegative least squares problem, under hard l0 sparsity constraints. This approach is based on a dedicated branch-and-bound algorithm. This simple strategy is able to compute the optimal solution even in complicated cases such as noisy or ill-conditioned data, where traditional approaches fail. We also show that our algorithm scales well, despite the combinatorial nature of the problem. We illustrate the advantages of the proposed technique on synthetic data sets, as well as a real-world hyperspectral image.