Optimizer’s dilemma: optimization strongly influences model selection in transcriptomic prediction

Motivation Most models can be fit to data using various optimization approaches. While model choice is frequently reported in machine-learning-based research, optimizers are not often noted. We applied two different implementations of LASSO logistic regression implemented in Python’s scikit-learn package, using two different optimization approaches (coordinate descent and stochastic gradient descent), to predict driver mutation presence or absence from gene expression across 84 pan-cancer driver genes. Across varying levels of regularization, we compared performance and model sparsity between optimizers. Results After model selection and tuning, we found that coordinate descent (implemented in the liblinear library) and SGD tended to perform comparably. liblinear models required more extensive tuning of regularization strength, performing best for high model sparsities (more nonzero coefficients), but did not require selection of a learning rate parameter. SGD models required tuning of the learning rate to perform well, but generally performed more robustly across different model sparsities as regularization strength decreased. Given these tradeoffs, we believe that the choice of optimizers should be clearly reported as a part of the model selection and validation process, to allow readers and reviewers to better understand the context in which results have been generated. Availability and implementation The code used to carry out the analyses in this study is available at https://github.com/greenelab/pancancer-evaluation/tree/master/01_stratified_classification. Performance/regularization strength curves for all genes in the Vogelstein et al. 2013 dataset are available at https://doi.org/10.6084/m9.figshare.22728644.


Introduction
Gene expression pro les are widely used to classify samples or patients into relevant groups or categories, both preclinically [1,2] and clinically [3,4]. To extract informative gene features and to perform classi cation, a diverse array of algorithms exist, and di erent algorithms perform well across varying datasets and tasks [1]. Even within a given model class, multiple optimization methods can often be applied to nd well-performing model parameters or to optimize a model's loss function. One commonly used example is logistic regression. The widely used scikit-learn Python package for machine learning [5] provides two modules for tting logistic regression classi ers: LogisticRegression , which uses the liblinear coordinate descent method [6] to nd parameters that optimize the logistic loss function, and SGDClassifier , which uses stochastic gradient descent [7] to optimize the same loss function.
Using scikit-learn, we compared the liblinear (coordinate descent) and SGD optimization techniques for prediction of driver mutation status in tumor samples, across a wide variety of genes implicated in cancer initiation and development [8]. We applied LASSO (L1-regularized) logistic regression, and tuned the strength of the regularization to compare model selection between optimizers. We found that across a variety of models (i.e. varying regularization strengths), the training dynamics of the optimizers were considerably di erent: models t using liblinear tended to perform best at fairly high regularization strengths (100-1000 nonzero features in the model) and over t easily with low regularization strengths. On the other hand, after tuning the learning rate, models t using SGD tended to perform well across both higher and lower regularization strengths, and over tting was less common.
Our results caution against viewing optimizer choice as a "black box" component of machine learning modeling. The observation that LASSO logistic regression models t using SGD tended to perform well for low levels of regularization, across diverse driver genes, runs counter to conventional wisdom in machine learning for high-dimensional data which generally states that explicit regularization and/or feature selection is necessary. Comparing optimizers or model implementations directly is rare in applications of machine learning for genomics, and our work shows that this choice can a ect generalization and interpretation properties of the model signi cantly. Based on our results, we recommend considering the appropriate optimization approach carefully based on the goals of each individual analysis.

Data download and preprocessing
To generate binary mutated/non-mutated gene labels for our machine learning model, we used mutation calls for TCGA Pan-Cancer Atlas samples from MC3 [9] and copy number threshold calls from GISTIC2.0 [10]. MC3 mutation calls were downloaded from the Genomic Data Commons (GDC) of the National Cancer Institute, at https://gdc.cancer.gov/about-data/publications/pancanatlas. Thresholded copy number calls are from an older version of the GDC data and are available here: https:// gshare.com/articles/dataset/TCGA_PanCanAtlas_Copy_Number_Data/6144122. We removed hypermutated samples, de ned as two or more standard deviations above the mean non-silent somatic mutation count, from our dataset to reduce the number of false positives (i.e., non-driver mutations). Any sample with either a non-silent somatic variant or a copy number variation (copy number gain in the target gene for oncogenes and copy number loss in the target gene for tumor suppressor genes) was included in the positive set; all remaining samples were considered negative for mutation in the target gene.
RNA sequencing data for TCGA was downloaded from GDC at the same link provided above for the Pan-Cancer Atlas. We discarded non-protein-coding genes and genes that failed to map, and removed tumors that were measured from multiple sites. After ltering to remove hypermutated samples and taking the intersection of samples with both mutation and gene expression data, 9074 total TCGA samples remained.

Cancer gene set construction
In order to study mutation status classi cation for a diverse set of cancer driver genes, we started with the set of 125 frequently altered genes from Vogelstein et al. [8] (all genes from Table S2A). For each target gene, in order to ensure that the training dataset was reasonably balanced (i.e., that there would be enough mutated samples to train an e ective classi er), we included only cancer types with at least 15 mutated samples and at least 5% mutated samples, which we refer to here as "valid" cancer types. In some cases, this resulted in genes with no valid cancer types, which we dropped from the analysis. Out of the 125 genes originally listed in the Vogelstein et al. cancer gene set, we retained 84 target genes after ltering for valid cancer types.

Classi er setup and optimizer comparison details
We trained logistic regression classi ers to predict whether or not a given sample had a mutational event in a given target gene using gene expression features as explanatory variables. Based on our previous work, gene expression is generally e ective for this problem across many target genes, although other -omics types can be equally e ective in many cases [11]. Our model was trained on gene expression data (X) to predict mutation presence or absence (y) in a target gene. To control for varying mutation burden per sample and to adjust for potential cancer type-speci c expression patterns, we included one-hot encoded cancer type and log 10 (sample mutation count) in the model as covariates. Since gene expression datasets tend to have many dimensions and comparatively few samples, we used a LASSO penalty to perform feature selection [12]. LASSO logistic regression has the advantage of generating sparse models (some or most coe cients are 0), as well as having a single tunable hyperparameter which can be easily interpreted as an indicator of regularization strength, or model complexity.
To compare model selection across optimizers, we rst split the "valid" cancer types into train (75%) and test (25%) sets. We then split the training data into "subtrain" (66% of the training set) data to train the model on, and "holdout" (33% of the training set) data to perform model selection, i.e. to use to select the best-performing regularization parameter, and the best-performing learning rate for SGD in the cases where multiple learning rates were considered. In each case, these splits were strati ed by cancer type, i.e. each split had as close as possible to equal proportions of each cancer type included in the dataset for the given driver gene.

LASSO parameter range selection and comparison between optimizers
The scikit-learn implementations of coordinate descent (in liblinear / LogisticRegression ) and stochastic gradient descent (in SGDClassifier ) use slightly di erent parameterizations of the LASSO regularization strength parameter. liblinear 's logistic regression solver optimizes the following loss function: where denotes the negative log-likelihood of the observed data given a particular choice of feature weights . SGDClassifier optimizes the following loss function: which is equivalent with the exception of the LASSO parameter which is formulated slightly di erently, as . The result of this slight di erence in parameterization is that liblinear values vary inversely with regularization strength (higher values = less regularization, or greater model complexity) and SGDClassifier values vary directly with regularization strength (lower values = less regularization, or greater model complexity).
For the liblinear optimizer, we trained models using values evenly spaced on a logarithmic scale between (10 -3 , 10 7 ); i.e. the output of numpy.logspace (-3, 7, 21) . For the SGD optimizer, we trained models using the inverse range of values between (10 -7 , 10 3 ), or numpy.logspace (-7, 3, 21) . These hyperparameter ranges were intended to give evenly distributed coverage across genes that included "under t" models (predicting only the mean or using very few features, poor performance on all datasets), "over t" models (performing perfectly on training data but comparatively poorly on cross-validation and test data), and a wide variety of models in between that typically included the best ts to the cross-validation and test data.
For ease of visual comparison in our gures, we plot the SGD parameter directly, and the liblinear parameter inversely (i.e. ). This orients the x-axes of the relevant plots in the same direction: lower values represent lower regularization strength or higher model complexity, and higher values represent higher regularization strength or lower model complexity, for both optimizers.

SGD learning rate selection
scikit-learn's SGDClassifier provides four built-in approaches to learning rate scheduling: constant (a single, constant learning rate), optimal (a learning rate with an initial value selected using a heuristic based on the regularization parameter and the data loss, that decreases across epochs), invscaling (a learning rate that decreases exponentially by epoch), and adaptive (a learning rate that starts at a constant value, which is divided by 5 each time the training loss fails to decrease for 5 straight epochs). The optimal learning rate schedule is used by default.
When we compared these four approaches, we used a constant learning rate of 0.0005, and an initial learning rate of 0.1 for the adaptive and invscaling schedules. We also tested a fth approach that we called " constant_search ", in which we tested a range of constant learning rates in a grid search on a validation dataset, then evaluated the model on the test data using the best-performing constant learning rate by validation AUPR. For the grid search, we used the following range of constant learning rates: {0.000005, 0.00001, 0.00005, 0.0001, 0.0005, 0.001, 0.01}. Unless otherwise speci ed, results for SGD in the main paper gures used the constant_search approach, which performed the best in our comparison between schedulers.

liblinear and SGD LASSO models perform comparably, but liblinear is sensitive to regularization strength
For each of the 125 driver genes from the Vogelstein et al. 2013 paper, we trained models to predict mutation status (presence or absence) from RNA-seq data, derived from the TCGA Pan-Cancer Atlas.
For each optimizer, we trained LASSO logistic regression models across a variety of regularization parameters (see Methods for parameter range details), achieving a variety of di erent levels of model sparsity (Supplementary Figure S1). We repeated model tting/evaluation across 4 cross-validation splits x 2 replicates (random seeds) for a total of 8 di erent models per parameter. Cross-validation splits were strati ed by cancer type.
Previous work has shown that pan-cancer classi ers of Ras mutation status are accurate and biologically informative [13]. We rst evaluated models for KRAS mutation prediction. As model complexity increases (more nonzero coe cients) for the liblinear optimizer, we observed that performance increases then decreases, corresponding to over tting for high model complexities/numbers of nonzero coe cients ( Figure 1A). On the other hand, for the SGD optimizer, we observed consistent performance as model complexity increases, with models having no nonzero coe cients performing comparably to the best ( Figure 1B). In this case, top performance for SGD (a regularization parameter of 10 -1 ) is slightly better than top performance for liblinear (a regularization parameter of 1 / 3.16 x 10 2 ): we observed a mean test AUPR of 0.722 for SGD vs. mean AUPR of 0.692 for liblinear .
To determine how relative performance trends with liblinear tend to compare across the genes in the Vogelstein dataset at large, we looked at the di erence in performance between optimizers for the best-performing models for each gene ( Figure 1C). The distribution is centered around 0 and more or less symmetrical, suggesting that across the gene set, liblinear and SGD tend to perform comparably to one another. We saw that for 52/84 genes, performance for the best-performing model was better using SGD than liblinear , and for the other 32 genes performance was better using liblinear . In order to quantify whether the over tting tendencies (or lack thereof) also hold across the gene set, we plotted the di erence in performance between the best-performing model and the largest (least regularized) model; classi ers with a large di erence in performance exhibit strong over tting, and classi ers with a small di erence in performance do not over t ( Figure 1D). For SGD, the least regularized models tend to perform comparably to the best-performing models, whereas for liblinear the distribution is wider suggesting that over tting is more common. Performance vs. regularization parameter for KRAS mutation status prediction, using the SGD optimizer. "Holdout" dataset is used for SGD learning rate selection, "test" data is completely held out from model selection and used for evaluation. C. Distribution of performance di erence between best-performing model for liblinear and SGD optimizers, across all 84 genes in Vogelstein driver gene set. Positive numbers on the x-axis indicate better performance using liblinear , and negative numbers indicate better performance using SGD. D. Distribution of performance di erence between best-performing model and largest (least regularized) model, for liblinear and SGD, across all 84 genes. Smaller numbers on the y-axis indicate less over tting, and larger numbers indicate more over tting.

SGD is sensitive to learning rate selection
The SGD results shown in Figure 1 select the best-performing learning rate using a grid search on the holdout dataset, independently for each regularization parameter. We also compared against other learning rate scheduling approaches implemented in scikit-learn (see Methods for implementation details and grid search speci cations). For KRAS mutation prediction, we observed that the choice of initial learning rate and scheduling approach a ects performance signi cantly, and other approaches to selecting the learning rate performed poorly relative to liblinear (black dotted lines in Figure 2) and to the grid search. We did not observe an improvement in performance over liblinear or the grid search for learning rate schedulers that decrease across epochs (Figure 2A, C, and D), nor did we see comparable performance when we selected a single constant learning rate for all levels of regularization without the preceding grid search ( Figure 2B). Notably, scikit-learn's default "optimal" learning rate schedule performed relatively poorly for this problem, suggesting that tuning the learning rate and selecting a well-performing scheduler is a critical component of applying SGD successfully for this problem ( Figure 2D). We observed similar trends across all genes in the Vogelstein gene set, with other learning rate scheduling approaches performing poorly in aggregate relative to both liblinear and SGD with the learning rate grid search (Supplementary Figure S2). Performance vs. regularization parameter, using SGD optimizer with constant learning rate scheduler and a learning rate of 0.0005. C. Performance vs. regularization parameter, using SGD optimizer with inverse scaling learning rate scheduler. D. Performance vs. regularization parameter, using SGD optimizer with "optimal" learning rate scheduler.

liblinear and SGD result in di erent models, with varying loss dynamics
We sought to determine whether there was a di erence in the sparsity of the models resulting from the di erent optimization schemes. In general across all genes, the best-performing SGD models mostly tend to have many nonzero coe cients, but with a distinct positive tail, sometimes having few nonzero coe cients. By contrast, the liblinear models are generally sparser with fewer than 2500 nonzero coe cients, out of ~16100 total input features, and a much narrower tail ( Figure 3A). The sum of the coe cient magnitudes, however, tends to be smaller on average across all levels of regularization for SGD than for liblinear ( Figure 3B). This e ect is less pronounced for the other learning rate schedules shown in Figure 2, with the other options resulting in larger coe cient magnitudes (Supplementary Figure S3). These results suggest that the models t by liblinear and SGD navigate the tradeo between bias and variance in slightly di erent ways: liblinear tends to produce sparser models (more zero coe cients) as regularization increases, but if the learning rate is properly tuned, SGD coe cients tend to have smaller overall magnitudes as regularization increases.
We also compared the components of the loss function across di erent levels of regularization between optimizers. The LASSO logistic regression loss function can be broken down into a datadependent component (the log-loss) and a parameter magnitude dependent component (the L1 penalty), which are added to get the total loss that is minimized by each optimizer; see Methods for additional details. As regularization strength decreases for liblinear , the data loss collapses to near 0, and the L1 penalty dominates the overall loss ( Figure 3C). For SGD, on the other hand, the data loss decreases slightly as regularization strength decreases but remains relatively high ( Figure 3D). Other SGD learning rate schedules have similar loss curves to the liblinear results, although this does not result in improved classi cation performance (Supplementary Figure S4).

Discussion
Our work shows that optimizer choice presents tradeo s in model selection for cancer transcriptomics. We observed that LASSO logistic regression models for mutation status prediction t using stochastic gradient descent were highly sensitive to learning rate tuning, but they tended to perform robustly across diverse levels of regularization and sparsity. Coordinate descent implemented in liblinear did not require learning rate tuning, but generally performed best for a narrow range of fairly sparse models, over tting as regularization strength decreased. Tuning of regularization strength for liblinear , and learning rate (and regularization strength to a lesser degree) for SGD, are critical steps which must be considered as part of analysis pipelines. The sensitivity we observed to these details highlights the importance of reporting exactly what optimizer was used, and how the relevant hyperparameters were selected, in studies that use machine learning models for transcriptomic data.
To our knowledge, the phenomenon we observed with SGD has not been documented in other applications of machine learning to genomic or transcriptomic data. In recent years, however, the broader machine learning research community has identi ed and characterized implicit regularization for SGD in many settings, including overparametrized or feature-rich problems as is often the case in transcriptomics [14,15,16]. The resistance we observed of SGD-optimized models to decreased performance on held-out data as model complexity increases is often termed "benign over tting": over t models, in the sense that they t the training data perfectly and perform worse on the test data, can still outperform models that do not t the training data as well or that have stronger explicit regularization. Benign over tting has been attributed to optimization using SGD [16,17], and similar patterns have been observed for both linear models and deep neural networks [18,19].
Existing gene expression prediction benchmarks and pipelines typically use a single model implementation, and thus a single optimizer. We recommend thinking critically about optimizer choice, but this can be challenging for researchers that are inexperienced with machine learning or unfamiliar with how certain models are t under the hood. For example, R's glmnet package uses a cyclical coordinate descent algorithm to t logistic regression models [20], which would presumably behave similarly to liblinear , but this is somewhat opaque in the glmnet documentation itself. Increased transparency and documentation in popular machine learning packages with respect to optimization, especially for models that are di cult to t or sensitive to hyperparameter settings, would bene t new and unfamiliar users.
Related to what we see in our SGD-optimized models, there exist other problems in gene expression analysis where using all available features is comparable to, or better than, using a subset. For example, using the full gene set improves correlations between preclinical cancer models and their tissue of origin, as compared to selecting genes based on variability or tissue-speci city [21]. On the other hand, when predicting cell line viability from gene expression pro les, selecting features by Pearson correlation improves performance over using all features, similar to our liblinear classi ers [22]. In future work, it could be useful to explore if the coe cients found by liblinear and SGD emphasize the same pathways or functional gene sets, or if there are patterns to which mutation status classi ers (or other cancer transcriptomics classi ers) perform better with more/fewer nonzero coe cients.

Data and code availability
The data analyzed during this study were previously published as part of the TCGA Pan-Cancer Atlas project [23], and are available from the NIH NCI Genomic Data Commons (GDC). The scripts used to download and preprocess the datasets for this study are available at https://github.com/greenelab/pancancer-evaluation/tree/master/00_process_data, and the code used to carry out the analyses in this study is available at https://github.com/greenelab/pancancerevaluation/tree/master/01_strati ed_classi cation, both under the open-source BSD 3-clause license. Equivalent versions of Figure 1A and 1B