We consider the colorization problem of grayscale images when some pixels, called scribbles, with initial colors are given. In this paper, we propose a new multi-layer graph model and an energy formulation that can incorporate higher-order cues for reliable colorization of natural images. In contrast to most existing energy functions with unary and pairwise constraints, we address the problem of imposing a high-order constraint whereby pixels constituting each region tend to have similar colors to the representative color of the region they belong to. The representative colors of the regions that are generated by unsupervised image segmentation algorithms, act as higher-order cues, and they are automatically obtained by a nonparametric learning technique that estimates them from the resulting pixel colors in a recursive fashion. We formulate this problem in terms of two quadratic energy functions of pixel and region colors, that are supplementary to each other, in our proposed multi-la...