This work deals with the improvement of multi-target prediction models through a proposed optimization called Selection Of medical Features by Induced Alterations in numeric labels (SOFIA). This method performs a data transformation when: (1) weighting the features, (2) performing small perturbations on numeric labels and (3) selecting the features that are relevant in the trained multi-target prediction models. With the purpose of decreasing the computational cost in the SOFIA method, we consider those multi-objective optimization metaheuristics that support parallelization. In this sense, we propose an extension of the Natural Optimization (NO) approach for Simulated Annealing to support a multi-objective (MO) optimization. This proposed extension, called MONO, and some multiobjective evolutionary algorithms (MOEAs) are considered when performing the SOFIA method to improve prediction models in a multi-stage migraine treatment. This work also considers the adaptation of these metaheuristics to run on GPUs for accelerating the exploration of a larger space of solutions and improving results at the same time. The obtained results show that accuracies close to 88% are obtained with the MONO metaheuristic when employing eight threads and when running on a GPU. In addition, training times have been decreased from more than 8 h to less than 45 min when running the algorithms on a GPU. Besides, classification models trained with the SOFIA method only require 15 medical features or fewer to predict treatment responses. All in all, the methods proposed in this work remarkably improve the accuracy of multi-target prediction models for the OnabotulinumtoxinA (BoNT-A) treatment, while selecting those relevant features that allow us to know in advance the response to every stage of the treatment.