In this article, we consider a bilaterally constrained optimization model arising from the semisupervised multiple-class image segmentation problem. We prove that the solution of the corresponding unconstrained problem satisfies a discrete maximum principle. This implies that the bilateral constraints are satisfied automatically and that the solution is unique. Although the structure of the coefficient matrices arising from the optimality conditions of the segmentation problem is different for different input images, we show that they are M-matrices in general. Therefore, we study several numerical methods for solving such linear systems and demonstrate that domain decomposition with block relaxation methods are quite effective and outperform other tested methods. We also carry out a numerical study of condition numbers on the effect of boundary conditions on the optimization problems, which provides some insights into the specification of boundary conditions as an input knowledge in the learning context. © 2010 Wiley Periodicals, Inc. Int J Imaging Syst Technol, 20, 191–201, 2010