The optimal feature selection problem in software product line is typically addressed by the approaches based on Indicator-based Evolutionary Algorithm (IBEA). In this study, we first expose the mathematical nature of this problem — multi-objective binary integer linear programming. Then, we implement/propose three mathematical programming approaches to solve this problem at different scales. For small-scale problems (roughly, less than 100 features), we implement two established approaches to find all exact solutions. For medium-to-large problems (roughly, more than 100 features), we propose one efficient approach that can generate a representation of the entire Pareto front in linear time complexity. The empirical results show that our proposed method can find significantly more non-dominated solutions in similar or less execution time, in comparison with IBEA and its recent enhancement (i.e., IBED that combines IBEA and Differential Evolution).