This paper proposes a Heaviside composite optimization approach and presents a progressive method for solving multi-treatment learning problems with non-convex constraints. A Heaviside composite function is a composite of a Heaviside function (i.e., the indicator function of either the open \(( \, 0,\infty )\) or closed \([ \, 0,\infty \, )\) interval) with a possibly nondifferentiable function. Modeling-wise, we show how Heaviside composite optimization provides a rigorous mathematical formulation for learning multi-treatment rules subject to Gini constraints. A Heaviside composite function has an equivalent discrete formulation and the resulting optimization problem can in principle be solved by integer programming (IP) methods. Nevertheless, for constrained treatment learning problems with large datasets, a straightforward application of off-the-shelf IP solvers is usually ineffective in achieving global optimality. To alleviate such a computational burden, our major contribution is the proposal of the progressive method by leveraging the effectiveness of state-of-the-art IP solvers for problems of modest sizes. We provide the theoretical advantage of the progressive method with the connection to continuous optimization and show that the computed solution is locally optimal for a broad class of Heaviside composite optimization problems. The superior numerical performance of the proposed method is demonstrated by extensive computational experimentation. A brief discussion of how score-based and tree-based multi-classification problems can also be formulated as Heaviside composite optimization problems and thus treated by the same progressive method is presented in an appendix.