extfieldsofstudy:
Computer Science
Mathematics

text:
Featurized Density Ratio Estimation

Density ratio estimation serves as an important technique in the unsupervised machine learning toolbox. However, such ratios are difficult to estimate for complex, high-dimensional data, particularly when the densities of interest are sufficiently different. In our work, we propose to leverage an invertible generative model to map the two distributions into a common feature space prior to estimation. This featurization brings the densities closer together in latent space, sidestepping pathological scenarios where the learned density ratios in input space can be arbitrarily inaccurate. At the same time, the invertibility of our feature map guarantees that the ratios computed in feature space are equivalent to those in input space. Empirically, we demonstrate the efficacy of our approach in a variety of downstream tasks that require access to accurate density ratios such as mutual information estimation, targeted sampling in deep generative models, and classification with data augmentation.

INTRODUCTION
A central problem in unsupervised machine learning is that of density ratio estimation: given two sets of samples drawn from their respective data distributions, we desire an estimate of the ratio of their probability densities [Nguyen et al., 2007, Sugiyama et al., 2012b. Computing this ratio gives us the ability to compare and contrast two distributions, and is of critical importance in settings such as out-of-distribution detection [Smola et al., 2009, Menon andOng, 2016], mutual information estimation [Belghazi et al., 2018, Song andErmon, 2019], importance weighting under covariate shift [Huang et al., 2006, Gretton et al., 2009, You et al., 2019, and hypothesis testing [Gretton et al., 2012]. Related * Denotes equal contribution. Figure 1: Flowchart for the featurized density ratio estimation framework. Direct density ratio estimation using a black-box algorithm DRE on samples leads to poor ratio estimatesr(x) when p and q are sufficiently different. By training a normalizing flow f θ on samples from both densities and encoding them to a shared feature space prior to estimation, we obtain more accurate ratios (r • f θ )(x).
areas of research which require access to accurate density ratios, such as generative modeling [Gutmann and Hyvärinen, 2010, Goodfellow et al., 2014, Nowozin et al., 2016 and unsupervised representation learning [Thomas et al., 2021], have enjoyed tremendous success with the development of more sophisticated techniques for density ratio estimation.
Despite its successes, density ratio estimation is an extremely hard problem when the two distributions of interest are considerably different [Cortes et al., 2010, Yamada et al., 2013, Rhodes et al., 2020. The fundamental challenge in reliably estimating density ratios in this scenario is precisely the access to only a finite number of samples. As the distance between the densities increases, we become less likely to observe samples that lie in low-density regions between the two distributions. Therefore, without an impractically large training set, our learning algorithm is highly likely to converge to a poor estimator of the true underlying ratio.
To address this challenge, we propose a general-purpose framework for improved density ratio estimation that brings the two distributions closer together in latent space. The key component of our approach is an invertible generative model (normalizing flow), which is trained on a mixture of datasets drawn from the two distributions and used to map samples into a shared feature space prior to ratio estimation [Rezende and Mohamed, 2015]. Encoding the data via the normalizing flow transforms the observed samples from the two densities to lie within a unit Gaussian ball. We observe that this contraction helps mitigate pathological scenarios where the learned ratio estimates are wildly inaccurate. The invertibility of our feature map then guarantees that the ratios computed in feature space are equivalent to those in input space. We demonstrate the generality of our framework by pairing it with several existing density ratio estimation techniques, and explore various training procedures in estimation algorithms that require learning a probabilistic classifier. A flowchart of our featurized density ratio estimation algorithm can be found on Figure 1.
Empirically, we evaluate the efficacy of our approach on downstream tasks that require access to accurate density ratios. First, we demonstrate that applying our approach to existing density ratio estimation techniques on synthetic data leads to better performance on downstream domain adaptation and mutual information (MI) estimation. Next, we demonstrate the utility of our framework on a targeted generation task on MNIST [LeCun, 1998]. By leveraging the "featurized" density ratios for importance sampling from a trained generative model, we show that the resulting samples are closer to the target distribution of interest than the synthetic examples generated using input-space density ratios. Finally, we illustrate that our method can be used to improve upon naive data augmentation methods by reweighing synthetic samples, outperforming relevant baselines on multi-class classification on Omniglot [Lake et al., 2015].
The contributions of our work can be summarized as: 1. We introduce a general-purpose algorithm for estimating density ratios in feature space and show its applicability to a suite of existing ratio estimation techniques.
2. By leveraging the invertibility of our feature map, we prove that our featurized density ratio estimator inherits key properties such as unbiasedness and consistency from the original ratio estimation algorithm.
3. On downstream tasks that require access to accurate density ratios, we show that our approach outperforms relevant baselines that compute ratios in input space.

INVERTIBLE TRANSFORMATIONS VIA NORMALIZING FLOWS
Deep invertible generative models, or normalizing flows, are a family of likelihood-based models that describe the two-way transformation between a complex, continuous probability density and a simple one by the change of variables formula [Rezende andMohamed, 2015, Papamakarios et al., 2019]. The flow is parameterized by a deep neural network f θ : X → Z with carefully designed architectures such that the overall transformation is composed of a series of bijective mappings with tractable inverses and Jacobian determinants [Dinh et al., 2016, Kingma et al., 2016, Papamakarios et al., 2017, Kingma and Dhariwal, 2018, Grathwohl et al., 2018, Ho et al., 2019. As a result, the density of the random variable X = f −1 θ (Z) can be evaluated exactly: : Z → X denotes the inverse of the mapping f θ , p denotes the probability density of the random variable X, and t denotes the probability density of Z. The base (prior) distribution t(z) is typically chosen to be an isotropic Gaussian N (0, I) -the simplicity of evaluating this prior density, coupled with the tractability of f −1 θ and its Jacobian, allows us to train the normalizing flow via maximum likelihood.
The key property of normalizing flows that we exploit in our method is their invertibility: the dimensionality of X and Z are the same by design, and data points can be losslessly mapped between the two spaces. As we will demonstrate in Section 3.2, this will be critical for translating the density ratios obtained in latent space back to those in input space.

DENSITY RATIO ESTIMATION TECHNIQUES
Notation and Problem Setup. We denote the input variable as x ∈ X ⊆ R d , and let z ∈ Z ⊆ R d be a latent variable of the same dimensionality as the input. We use capital letters to denote random variables, e.g. X p ∼ p and X q ∼ q. Then, Z p = f θ (X p ) and Z q = f θ (X q ) denote the random variables obtained by transforming X p and X q with f θ . We note that since Z p and Z q are transformed by the same normalizing flow f θ , we can form the mixture density The learning setting we consider is as follows. Given two sets of observed samples , we wish to estimate the ratio of their underlying probability densities r(x) = p(x)/q(x). We focus on direct ratio estimation techniques, where we learn the density ratio estimatorr rather than constructing explicit density estimates ofp(x) andq(x) and computing their ratio [Sugiyama et al., 2012b]. The estimatorr is obtained via an estimation algorithm DRE which takes as input two datasets and returns a function DRE(D p , D q ) =r : X → R. Then, evaluatingr at a particular point x gives us an estimate of the true density ratior(x) ≈ r(x). In the following exposition, we provide background information on the suite of existing density ratio estimation algorithms. Bottom row: Visualizations of the learned encodings for various training strategies for (f-h), with ground truth samples from p(x) and q(x) in (e). We note that using a pretrained flow as an invertible encoder as in (b) leads to the most accurate density ratio estimates.
For KMM, density ratio estimates are obtained by projecting all inputs into a reproducing kernel Hilbert space (RKHS) H induced by a characteristic kernel k : X × X → R.
Although several choices for the kernel are possible, Huang et al. use the Gaussian kernel k(x, x ) = exp(||x − x || 2 ) to arrive at the following objective: where both expectations are approximated via Monte Carlo. Intuitively, KMM attempts to match the mean embedding of the two distributions (where the embedding is produced by the canonical feature map defined by k) in H.
For KLIEP, the goal is to estimate density ratios r(x) = p(x)/q(x) such that the Kullback-Leibler (KL) divergence between p(x) andp(x) =r(x)q(x) is minimized: The solution to this constrained optimization problem can also be obtained in H by parameterizingr and some kernel k, similar in spirit to KMM.
Probabilistic Classification. Another technique to obtain density ratio estimates is via probabilistic classification, in which a binary classifier c φ : X → [0, 1] is trained to discriminate between samples from the two densities p(x) and q(x) [Friedman et al., 2001, Gutmann and Hyvärinen, 2010, Sugiyama et al., 2012b. Concretely, suppose we construct a dataset such that all samples D p are given the pseudolabel y = 1, and those from D q are labeled as y = 0. Assuming that the two datasets are equal in size (n p = n q , though this can be relaxed with a scaling factor), we can use Bayes' rule to arrive at the following expression for the density ratio: where c * φ (x) = P (y = 1|x) denotes the Bayes optimal classifier for this particular task.

MOTIVATING EXAMPLE AND INTUITION
Despite the suite of existing techniques for density ratio estimation, they are of limited use in settings where p(x) and q(x) are mismatched in support [Cortes et al., 2010, Yamada et al., 2013, You et al., 2019, Rhodes et al., 2020 We highlight an illustrative failure case in Figure 2 on a 2-dimensional toy dataset, where p(x) ∼ N ([0, 0] T , I) and q(x) ∼ N ([3, 3] T , I). As shown in Figure 2(e), the two random variables have regions of minimal overlap -when training a binary classifier c φ to distinguish the two sets of samples, the log-ratio estimates logr φ (x) learned by the classifier are noticeably inaccurate (Figure 2(a)).
To develop a solution, we consider a simple example to build intuition about the problem (with more details in Appendix B). Suppose we want to estimate the density ratios between two 1-dimensional Gaussians, p ∼ N (m, 1) and q ∼ N (−m, 1), with a finite number of samples and D q = {x q i } n i=1 of size n from each. The analytic solution for r(x) = p(x)/q(x) = exp − (x−m) 2 −(x+m) 2 2 = exp(2mx) for x ∈ R, which grows exponentially with m. Without access to the parametric forms of p and q, we train a logistic regression model c φ to discriminate between D p and D q , where the maximum likelihood objective is: where σ(z) = 1/(1 + exp(−z)). Although the logistic regression model is well-specified in this setting, and can achieve Bayes optimal risk in the limit of infinite data, we illustrate what can go wrong in the finite sample regime.
Suppose that m > 0 is large -there exists a large separation between p and q. Then, most samples D p ∼ p will take on positive values, and most samples D q ∼ q will be negative. In this situation, the model will be incentivized to push w 1 → ∞ to maximize the objective. This will lead to wildly inaccurate density ratio estimates, as we know that the true values of w 1 = 2m and w 0 = 0 are far from infinity (in fact, r(x) = exp(w 1 · x)). Thus we must see samples between p and q during training: concretely, samples from p such that x p ≤ 0 and samples from q such that x q ≥ 0. But with n samples from p, the probability that D p contains all positive samples is n i=1 P (X p > 0) ≥ (1 − n exp(−m 2 /2)), which means that the number of samples required to avoid pathological solutions is exponential in m 2 . This implies that density ratio estimation via probabilistic classification in input space is near impossible in such scenarios without extremely large amounts of training data.

METHODOLOGY
The motivating example in the previous section suggests that it is critical to bring p and q "closer together" to make the density ratio estimation problem tractable. Our solution is to do so in latent space by leveraging an invertible transformation. Concretely, we consider training an invertible deep generative model f θ on a mixture of p(x) and q(x), such that f θ (X p ) and f θ (X q ) are mapped to a common feature space Z. The result of utilizing f θ as invertible feature map can be visualized in Figure 2(f): the flow compresses all data points to lie in different regions of a unit Gaussian ball. By mapping regions of low density in X into regions of higher density in Z, and training our probabilistic classifier c φ on f θ (x) ∈ Z rather than x ∈ X directly, this contraction leads to learning more accurate density ratios as shown in Figure 2(b). We refer the reader to Appendix F for additional experimental details and results.
Our first observation is that as a direct consequence of the invertibility of f θ , the density ratios obtained in feature space are equivalent to those obtained in input space. We formalize this statement in Lemma 1 below.
Lemma 1. Let X p ∼ p be a random variable with density p, and X q ∼ q be a random variable with density q. Let f θ be any invertible mapping. Let p , q be the densities of Z p = f θ (X p ) and Z q = f θ (X q ) respectively. Then for any x: Proof. We provide the proof in Appendix C.1.
This simple observation is quite powerful, as it lends us a general-purpose algorithm that may improve many existing ratio estimation techniques as a black-box wrapper. We provide the pseudocode for our training procedure in Algorithm 1. Given the two sets of samples, the ratio estimation method, and an invertible generative model family, we first train the normalizing flow on a mixture of the two datasets (Lines 2-3). We then use the trained flow to encode the samples into a common feature space and plug them into the base density ratio estimator algorithm DRE(·) to obtain our featurized density ratio estimatorr • f * θ , which is implicitly composed with the trained normalizing flow (Line 6). This algorithm allows us to lightly modify existing approaches such as KMM and KLIEP as detailed in Appendix A, and we explore their featurized variants in our experiments.

TRAINING PROCEDURE
In practice, there are a variety of ways to implement the training procedure as outlined in Algorithm 1. The most general is separate training, which leverages a pre-trained flow f θ as an invertible encoder to map the inputs into a common feature space, prior to ratio estimation. This approach is capable of handling all parametric and nonparametric techniques which operate directly in input space.
In the probabilistic classification setting, where the density ratio estimation algorithm DRE(·) requires learning a binary classifier c φ to distinguish between D p and D q , we can adapt the normalizing flow f θ to account for the known structure of c φ . We call this procedure joint training. Both the normalizing flow f θ and the discriminative classifier c φ are trained jointly via the following objective: where L sup denotes the standard binary cross entropy (logistic) loss, L flow denotes the maximum likelihood objective for the flow f θ , and α ∈ [0, 1] is a hyperparameter which balances the importance of the two terms in the loss function. This approach is quite common in learning deep hybrid models [Kuleshov andErmon, 2017, Nalisnick et al., 2019].
Finally, we explore discriminative training, where we modify the classifier c φ 's architecture to incorporate that of the flow f θ to build an "invertible" classifier c φ,θ : X → [0, 1] that is trained solely via the logistic loss L sup (θ, φ). This is inspired by the strong performance of invertible networks such as i-

CHARACTERIZATION OF THE LEARNED FEATURE SPACE
At a first glance, Lemma 1 appears to suggest that any feature space induced by an invertible map should work well for density ratio estimation, as long as f θ (X p ) and f θ (X q ) are closer together than X p and X q . To gain further insight into the desirable characteristics of the learned feature space, we visualize the encodings of the various training strategies in Figure 2. For both a pretrained (Figure 2(f)) and jointly trained ( Figure 2(g)) normalizing flow, the data points are mapped to lie closer together in different regions of the unit Gaussian ball. However, for the discriminatively trained classifier equipped with an invertible "encoder" (Figure 2(h)), the encoded examples more closely resemble the shape of the original inputs ( Figure 2(e)). This observation, combined with the low quality density ratio estimates in Figure 2(d) relative to the other training methods (Figure 2(b-c)), suggests that maximum likelihood training of the normalizing flow f θ in addition to shrinking the gap between the densities p and q is crucial for obtaining accurate density ratios in feature space. We hypothesize that mapping the observations into a unit Gaussian ball is an important property of our method, and we save an in-depth theoretical analysis of this phenomenon for future work.

THEORETICAL ANALYSIS
In this section, we provide theoretical justifications for several properties of the featurized density ratio estimator. As a consequence of Lemma 1, we find that our estimator inherits many of the desirable properties of the original estimator.

PROPERTIES OF THE ESTIMATOR
Unbiasedness. Unbiasedness is one of the most fundamental desiderata of a statistical estimator, as it guarantees that the estimated parameter is equivalent to the parameter's true value in expectation. In Corollary 1, we prove that unbiasedness of the featurized ratio estimator follows directly if the original estimator is also unbiased.
and any p, q, and let f θ denote any invertible mapping. Then, is also an unbiased estimator of p(x) q(x) for any p, q.
Proof. We provide the proof in Appendix C.2.
Consistency. Consistency is another key property in a statistical estimator, as it guarantees that in the limit of infinite data used in the estimation procedure, the probability that the estimator becomes arbitrarily close to the true parameter converges to one. We prove in Corollary 2 that consistency of the featurized density ratio estimator also follows if the original density ratio estimator is consistent. This is desirable, as estimators such as the KLIEP and KMM (with universal kernels) are both consistent [Huang et al., 2006, Gretton et al., 2009, Sugiyama et al., 2012b.
Corollary 2. Let D p be n i.i.d samples from density p, and D q be n i.i.d samples from density q. Letr(x) obtained fromr = DRE(D p , D q ) be a consistent estimator of p(x) q(x) for all x ∈ X and for any p, q. Let f θ be any invertible mapping. Then, q(x) for any p, q.
Proof. We provide the proof in Appendix C.3.

EXPERIMENTAL RESULTS
In this section, we are interested in empirically investigating the following questions: 1. Are the density ratios learned in feature space indeed more accurate than those learned in input space?
2. Do estimates in feature space yield better performance on downstream tasks that rely on density ratios?
For conciseness, we report the average over several runs for all experiments and report complete results in Appendix F.

Datasets.
We evaluate the efficacy of featurized density ratio estimation on both synthetic and real-world datasets. Models. We train different classifiers depending on the difficulty of the classification task, but largely keep the same architecture (either an MLP or CNN) across different tasks.
For the normalizing flow, we utilize the Masked Autoregressive Flow (MAF) for all datasets [Papamakarios et al., 2017]. We train the MAF separately on the mixture of the two datasets prior to density ratio estimation for all experiments with the exception of the MI estimation experiment in Section 5.2, where we explore various training strategies mentioned in Section 3.3. For additional details regarding architecture design and relevant hyperparameters, we refer the reader to Appendix E.

DOMAIN ADAPTATION
We first pair our method with two existing techniques, KMM and KLIEP, to assess whether estimating ratios in feature space improves performance on domain adaptation tasks with: 1) 2-D Gaussian mixtures and 2) the UCI Breast Cancer dataset. On the synthetic dataset, our method achieves a lower test error than both baseline logistic regression (without importance weighting) and reweighted logistic regression using density ratios estimated by KMM and KLIEP in input space. See Appendix F.2 for full results.
The UCI Breast Cancer dataset consists of n = 699 examples from 2 classes: benign (y = 1) and malignant (y = −1), where each sample is a vector of 9 features. We replicate the experimental setup of [Huang et al., 2006] to construct a source dataset with a heavily downsampled number of benign labels, while leaving the target dataset as is. After learning the importance weights via density ratio estimation on a mixture of the source and (unlabeled) target datasets, we train a support vector machine (SVM) with a Gaussian kernel of bandwidth σ = 0.1 and varying penalty hyperparameter values C = {0.1, 1, 10, 100} with importance weighting on the source domain. The binary classifier is then tested on the target domain. As shown in Figure 3, when applied to KMM, for nearly all values of C, our method (z-dre) achieves the lowest test error on the target dataset compared to both a vanilla SVM (baseline) and a reweighted SVM with density ratio estimates computed in input space (x-dre). Additionally, we note that our method achieves the absolute lowest test error across varying values of C. We report the average values of our KMM experiments over 30 runs in Figure 3.
All methods performed poorly overall for our KLIEP experiments. This result aligns with many past works with KLIEP importance-weighted classification; empirically, KLIEP only outperforms baseline unweighted classifiers on synthetic data, while on more complex datasets (e.g. UCI), KLIEP shows no significant improvements [Sugiyama et al., 2007, Tsuboi et al., 2008, Yamada and Sugiyama, 2009, Loog, 2012. In order to confirm the consistency of this behavior, we performed an additional experiment with a slightly different dataset-biasing process in which data points that were further from the mean were selected less often, similarly to Huang et al. [2006]; we report more details on the biased subsampling process in Appendix D.2. We used two datasets: 1) the UCI Blood Transfusion dataset and 2) the UCI Wine Quality dataset and found that both reweighted classifiers performed similarly to the baseline. Notably, our z-dre method does not degrade the performance of KLIEP. Table 1 shows our results.

MUTUAL INFORMATION ESTIMATION
Next, we test our approach on a mutual information (MI) estimation task between two correlated 20-dimensional Gaussian random variables, where the ground truth MI is tractable. MI estimation between two random variables X p and X q is a direct application of density ratio estimation, as the problem can be reduced to estimating average density ratios between their joint density and the product of their marginals. If we let v denote the joint density of X p and X q , we can see that: . We adapt the experimental setup of [Belghazi et al., 2018, Poole et al., 2019, Song and Ermon, 2019 to use a correlation coefficient of ρ = 0.9. We further explore the effect of the various training strategies as outlined in Section 3.3. While we use a MAF as the normalizing flow for all configurations, we evaluate our approach against: (a) the baseline classifier (baseline); (b) the two-stage approach (separate), where the flow is trained first on a mixture of D p and D q before training the classifier on the encoded data points; (c) jointly training the flow and the classifier (joint); and (d) a purely discriminative approach where the classifier architecture has a flow component (disc-only). For joint training, we sweep over α = {0.1, 0.5, 0.9}. As shown in Figure 4, the probabilistic classifier trained in feature space (after encoding the data using the normalizing flow) via our method outperforms relevant baselines. Interestingly, we find that for the joint training, higher values of α (which places a greater emphasis on the classification loss L sup rather than L flow as in Eq. 1) leads to more accurate MI estimates. For additional details on the data generation process and experimental setup, we refer the reader to Appendix E. While separate training outperforms all baselines, joint training achieves competitive performance with larger α.

TARGETED GENERATION WITH MNIST
For this experiment, we evaluate the effectiveness of our learned density ratio estimates on a targeted generation task using the MNIST dataset. Our goal is to generate samples according to a target distribution q(x) in a data-efficient manner, given samples from both p(x) and q(x). We test two scenarios: (a) diff-digits: a subset of MNIST in which p(x) is comprised of the digits labeled {1,2}, and q(x) which is comprised of the digits labeled {0,7}; (b) diff-background: a setting in which p(x) contains the original MNIST digits (black background, white digits); and q(x) contains the same examples but with flipped colors (white background, black digits). The second scenario is trickier than the first, since there exists an obvious gap between the two distributions. We also explore the effect of the target dataset size q(x) in learning accurate density ratios. Following the setup of [Choi et al., 2020], we sweep over various sizes of q(x) relative to p(x), which we call perc={0.1, 0.25, 0.5, 1.0} (where perc = 0.5 indicates that D q is 50% the size of D p ). After training a MAF on both D p and D q and obtaining density ratio estimates (importance weights), we sample from the trained MAF via sampling-importance-resampling (SIR) [Liu andChen, 1998, Doucet et al., 2000] at generation time.
As shown in Table 2, we achieve greater success in the targeted generation task when performing SIR with importance weights learned in feature space. Averaged across perc={0.1, 0.25, 0.5, 1.0} with 1000 generated samples each, our method generates 19.1% more samples from q(x) relative to the pretrained flow and 6.7% more samples than the baseline with importance weights learned in input space on the diff-digits task. Similarly for the diff-background task, our framework generates 18.8% more samples from q(x) relative to the pretrained flow and 16.4% more samples than the baseline. For additional experimental details, as well as the generated samples, we refer the reader to Appendix D and F. We first train a MAF on a mixture of the training examples and generated samples, encode all the data using the flow, and obtain importance weights via the encodings. The importance weights are obtained by training a binary classifier on the featurized inputs. We experiment with different baselines: (a) training the classifier without any data augmentation (Data-only); (b) training the classifier on purely synthetic samples (Synthetic-only); (c) training the classifier with data-augmentation without any importance weighting (Mixture-only); (d) the data-augmented classifier with importance weights obtained from input space (Mixture + IW(x)); and (e) the data-augmented classifier with importance weights obtained from feature space (Mixture + IW(z)). As shown in Table 3, the importance weights learned in the feature space show a significant boost in overall downstream classification accuracy as compared to relevant baselines: our method improves 3.7% over the Data-only baseline, and 2.2% over the highest performing baseline. We refer the reader to Appendix F for additional experimental details and results.

RELATED WORK
Density Ratio Estimation in Feature Space. Although density ratio estimation in machine learning has an extremely rich history [Friedman et al., 2001, Huang et al., 2006, Nguyen et al., 2007, Gutmann and Hyvärinen, 2010, Sugiyama et al., 2012b, there is considerably less work exploring the method's counterpart in feature space. [Rhodes et al., 2020], while tackling the same problem of density ratio estimation between two different data distributions, adopts a different approach than our framework. In particular, they propose a divide-and-conquer solution by constructing intermediate distributions between the two densities p(x) and q(x), and requires the training of a multi-task logistic regression model rather than a single binary classifier. Their interpolation technique, which is also conducted in the latent space of a normalizing flow in one of their experiments, is complementary to our framework -investigating the combination of these two approaches would be interesting future work. Additionally, density ratio estimation (in the form of learning importance weights) has been popular in a variety of domain adaptation approaches such as [Bickel et al., 2007, Long et al., 2015, You et al., 2019 which leverage a feature extractor to project the inputs into a lower-dimensional manifold prior to estimation. Although our approach shares a similar idea, the invertibility of our feature map guarantees that the density ratios between input space and feature space are equivalent -this is not necessarily true if the inputs are lossily compressed. Although we empirically verify that this is also the case for our method, it is more general precisely because the best partitioning is learned by the model.

CONCLUSION
In this paper, we proposed a general-purpose framework for improved density ratio estimation in settings where the two underlying data distributions of interest are sufficiently different. The key component of our approach is a normalizing flow that is trained on a mixture of the data sources, which is then used to encode the data into a shared feature space prior to estimating density ratios. By leveraging the invertibility of the flow, we showed that the ratios of the densities in feature space are not only identical to those in input space, but are also easier to learn. Additionally, our method is applicable to a suite of existing density ratio estimation techniques. Empirically, we demonstrated the utility of our framework on various combinations of density ratio estimation techniques and downstream tasks that rely on accurate density ratios for good performance, such as domain adaptation, mutual information estimation, and targeted generation in deep generative models. We provide a reference implementation in PyTorch [Paszke et al., 2017], and the codebase for this work is open-sourced at https://github.com/ermongroup/f-dre.
One limitation of our method is the need to train a normalizing flow on a mixture of the two datasets if a pre-trained model is not available; this may be difficult if the generative model must be extremely high-capacity. For future work, it would be interesting to explore whether the necessity of strict invertibility of the flow can be relaxed, and to gain a deeper theoretical understanding of the role of maximum likelihood training in our framework.

APPENDIX A FEATURIZED KMM AND KLIEP
Similar in spirit to the probabilistic classification approach in Section 2.2, we note that it is quite straightforward to extend this technique to non-parametric density ratio estimation methods. Suppose that (r • f θ ) is obtained fromr = DRE (f θ (D p ), f θ (D q )). Then, we find that the solution to KMM is equivalent after we first map the inputs into the feature space via f θ : X → Z: For KLIEP, we may also solve for the density ratio estimates in feature space: as a straightforward consequence of Lemma 1.

B DERIVATIONS FOR MOTIVATING EXAMPLE
We derive the calculations from the simple example presented in Section 3.1, and restate the problem setting here for completeness. Suppose we want to estimate the density ratios between two Gaussians, p ∼ N (m, 1) and q ∼ N (−m, 1), with a finite number of samples D p and D q of size n from each. We denote the random variable with the density p as X p , and the random variable with density q as X q . Our intuition was that as m grows larger, the probability that we would observe all positive samples from p (and analogously all negative samples from q) would be extremely high.
Without loss of generality, we first compute P (X p ≤ 0) by first using the well-known (lower) tail bound for Gaussian random variables: since the minimum is achieved at θ * = x − m, where ψ(θ) = exp(θm + θ 2 /2) is the moment generating function for N (m, 1). This tells us that the probability of observing a single positive sample from p is P (X p > 0) = 1 − P (X p ≤ 0) ≥ 1−exp(−m 2 /2), so taking account the fact that we have n i.i.d. samples gives us: n i=1 P (X p > 0) ≥ (1−exp(−m 2 /2)) n . Next, we compute the probability of seeing a single sample in our training set such that X p ≤ 0. Our reasoning was that such observed examples would help mitigate the pathological behavior of our learning algorithm driving up the magnitude of the logistic regression parameters to infinity. We find that: which is an extremely low probability. In fact, if we set 1 − (1 − exp(−m 2 /2)) n = δ and solve for δ, we find that: (1 − exp(−m 2 /2)) n = 1 − δ n log(exp(−m 2 /2)) = log(1 − δ) n = log(1 − δ) log(1 − exp(−m 2 /2)) Therefore, we observe a non-positive sample from p with probability at most δ for n < log(1−δ) log(1−exp(−m 2 /2)) . For a perhaps more intuitive bound, we can use Bernoulli's inequality, which states that (1 + x) r ≥ 1 + r · x for x ≥ −1, r ∈ R/ (0, 1). Doing so, we see that: which indicates that we require a training set size that is exponential in the order of m 2 to avoid the pathological scenario described in Section 3.1.

C.1 Proof of Lemma 1
For completeness, we restate Lemma 1 prior to providing the proof.
Lemma 1. Let X p ∼ p be a random variable with density p, and X q ∼ q be a random variable with density q. Let f θ be any invertible mapping. Let p , q be the densities of Z p = f θ (X p ) and Z q = f θ (X q ) respectively. Then for any x: Proof. By the change of variables formula:

C.2 Proof of Unbiasedness for the Featurized Density Ratio Estimator (Corollary 1)
For completeness, we restate Corollary 1 prior to providing the proof.
Corollary 1. Let D p be n i.i.d samples from density p, and D q be n i.i.d samples from density q. Letr(x) obtained from r = DRE (D p , D q ) be an unbiased estimator of r(x) = p(x) q(x) and any p, q, and let f θ denote any invertible mapping. Then, is also an unbiased estimator of p(x) q(x) for any p, q.
Proof. Using the definition of unbiasedness, we have: Let p , q be the densities of f θ (X p ) and f θ (X q ), respectively. Consider the estimatorr = DRE(f θ (D p ), f θ (D q )) which is an unbiased estimator of p (f θ (x)) q (f θ (x)) by assumption. Then: By the definition of p , q , this is equivalent to: where the last equality follows from Lemma 1.

C.3 Proof of Consistency (Corollary 2)
For completeness, we restate Corollary 2 before the proof statement.
Corollary 2. Let D p be n i.i.d samples from density p, and D q be n i.i.d samples from density q. Letr(x) obtained from r = DRE(D p , D q ) be a consistent estimator of p(x) q(x) for all x ∈ X and for any p, q. Additionally, let f θ be any invertible mapping. Then, (r • f θ )(x) obtained fromr = DRE (f θ (D p ), f θ (D q )) is also a consistent estimator of p(x) q(x) for any p, q.
Proof. By the definition of consistency, we have that ∀x ∈ X and ∀ > 0: Let p , q be the densities of f θ (X p ) and f θ (X q ) respectively. Because the estimator is assumed to be consistent for any p, q: and by definition of p , q this is equivalent to: Because the condition holds ∀x ∈ X , we have: where the last equality is due to Lemma 1.

D.1 Miscellaneous Background Information
Data Preprocessing. Prior to training the MAF, we: (a) use uniform dequantization; (b) rescale the pixels to lie within [0,1], and apply the logit transform following [Papamakarios et al., 2017]. For classification, we simply rescale the pixels to lie within [0,1].
Importance Weighting in Practice. As noted in [Grover et al., 2019], we apply two techniques when using the learned density ratio estimates as importance weights in our experiments.
1. Self-normalization: As a way to reduce variance, we ensure that the importance weights in a batch of n examples sum to 1, as in the expression below:r j=1r (x j ) We find that this technique works quite well when estimating density ratios in input space.
2. Flattening: we raise our obtained density ratio estimates to the power of a scaling parameter γ ≥ 0: Empirically, we observe that this approach works best on the ratios obtained in feature space.

D.2 KMM and KLIEP
Code. For our experiments using KMM and KLIEP, we based our code on the following implementations: • https://github.com/sksom/Classification-using-KMM-Kernel-Mean-Matching-• https://github.com/srome/pykliep Datasets. We used two datasets for both the KMM and KLIEP experiments: a generated 2D mixture of Gaussians dataset, and the Breast Cancer Wisconsin Data Set from the UCI Archive [Dua and Graff [2017]]. For each dataset, we construct our source and target splits as follows: • 2D Mixture of Gaussians: For our source dataset, we sampled 10 points from N (0, I) and 990 points from N ([3, 3] T , I), and for our target dataset, we sampled 990 points from N (0, I) and 10 points from N ([3, 3] T , I).
• Breast Cancer: Each sample consists of 9 input features (each of which with values ranging from 0 − 9) and one binary label. For each of n = 30 trials, we first set aside 3/4 of the dataset for our target dataset and then, with the remaining 1/4 of the data, constructed a biased source dataset by subsampling the training data according to P (s = 1|y = 1) = 0.1 and P (s = 1 | y = −1) = 0.9, where s indicates whether or not we include the sample. After subsampling, we normalized each feature value to be mean 0 and variance 1 (the same as in [Huang et al., 2006]).
• Blood Transfusion: This dataset consists of 748 samples (each corresponding to one person) with 5 input features and one binary label that represents whether or not the person is a blood donor. For each of n = 30 trials, as with the Breast Cancer dataset, we set aside 3/4 of the dataset for the target dataset and used the remaining 1/4 of the data to construct a biased source dataset by subsampling Huang et al., 2006]).
• Wine Quality: This dataset consists of 4898 samples with 12 input features and a label between 0 − −10 representing the wine quality. The binary classification task was the prediction of whether or not the wine quality was ≥ 5. We followed the same subsampling setup as for the Blood Transfusion dataset.

D.4 Targeted Generation with MNIST
We note that a normalizing flow model that has been trained on any mixture of D p and D q can be adapted for downstream applications of density ratio estimation. Concretely, we consider importance sampling, where we are interested in computing a statistic of the data g(·) with respect to a target distribution p(x): where the flow has been trained on an equal-sized mixture of D p and D q , the distribution learned by the flow is denoted as h(x) = 1 2 p(x) + 1 2 q(x), and the importance weight (learned density ratio estimate) r(x) has been re-balanced to account for the mixing proportions of p(x) and q(x) in the trained flow: r (x) = r(x) 1 2 (r(x)+1) . In the case that the mixing proportions are different (e.g. D p and D q are of different sizes), the re-balanced importance weight r (x) can be adjusted accordingly. We use this reweighting procedure in the MNIST targeted sampling experiments in Section 5.3.

D.5 Multi-class Classification with Omniglot
For the multi-class classification, we used the CNN-based architecture in [Vinyals et al., 2016] as shown in Table 7. For data augmentation, we randomly sampled 50 examples for each of the 1200 classes from the trained DAGAN -thus for all other models aside from the Data-only baseline, the training set size increased from (1200*10) to (1200*60).
For importance weighting, we trained both binary classifiers and input-space and feature-space to distinguish between the real and synthetic examples. We applied early stopping to the density ratio classifiers based on the validation set, which was comprised of 5 real examples and 5 synthetic examples. For the input-space density ratio estimation classifier, we found that the self-normalization technique worked best. For the feature-space density ratio estimation classifier, however, we found that flattening with γ = 0.2 worked well, and used this configuration. Additional details on self-normalization and flattening can be found in Appendix D.1.  Table 4. We use the sequential input ordering with ReLU activations and batch normalization between the blocks. We build on top of a pre-existing PyTorch implementation (https://github.com/kamenbliznashki/normalizing_flows).
Hyperparameters. During training, we use a batch size of 100 and the PyTorch default values of the Adam optimizer with learning rate = 0.0001 and weight decay of 1e-6 for all datasets. We use early stopping on the best log-likelihood on a held-out validation set.

E.2 MLP Classifier
We utilize the following MLP classifier architecture as shown in Table 5 for several of our experiments: (a) the synthetic 2-D Gaussians setup in Section 3.2; (b) the mutual information estimation experiment; and (c) the attribute classifier for the targeted MNIST generation task.

Name Component Input Layer
Linear in_dim → h_dim, ReLU Hidden Layer #1 Linear h_dim → h_dim, ReLU Hidden Layer #2 Linear h_dim → h_dim, ReLU Output Layer Linear h_dim → out_dim   We note that for the attribute classifier for MNIST, we explored two scenarios: • diff-digits, where all digits of classes {1,2} were given the label y = 0, and digits of classes {0,7} were labeled as y = 1 • diff-background, where all digits from the original dataset were labeled as y = 0 and those with flipped colors (white background, black digits) were labeled as y = 1.
In order to distinguish the separate classes for targeted generation, an MLP-based classifier was trained for each of the diff-digits and diff-background tasks as outlined in Tables 5 and 6.

E.3 Density Ratio Classifier
Depending on the complexity of the dataset, we used either an MLP classifier (Table 5) or CNN-based classifier (Table 7) for the density ratio estimator. For all synthetic experiments including those conducted on the MNIST dataset, we used an MLP for both input-space and feature-space density ratio estimation. For the Omniglot experiments, we used a slightly modified version of the CNN-based classifier where we swap the final output layer to be a Linear layer of dimension 64 → 1.
Hyperparameters. During training, we use a batch size of 64 and the Adam optimizer with learning rate = 0.001. The classifiers learn relatively quickly for both scenarios and we only needed to train for 10 epochs.

E.4 Downstream Classifier for Omniglot
For the multi-class classification task with Omniglot, we leveraged a commonly-used CNN architecture following [Vinyals et al., 2016], as shown in the following  Hyperparameters. During training, we sweep over batch sizes of {32,64,128} and the Adam optimizer with learning rate = 0.001. We also swept over the flattening coefficient for density ratio estimation and found that γ = 0.2 worked best. We trained the classifier for 100 epochs, and used early stopping on the validation set of Omniglot to determine the best model for downstream evaluation.

F.1 Toy Gaussian Mixture Experiment
We provide additional experimental results on the motivating 2-D Gaussian mixture example introduced in Section 3.2, where we sweep through additional values of α = {0.01, 0.1, 0.5, 0.7} on top of the one explored in the main text (α = 0.9). For reference, Figure 5 displays the (a) ground truth data and log density ratios (b-c) that we hope to learn from samples. Results are shown in Figure 6. A visual inspection of the 4 joint training procedures demonstrates that for this experiment, the jointly trained flow with the smallest contribution of the classification loss (α = 0.01 in (a)) outperforms all other methods (b-d). The learned feature space most closely resembles that of the separately trained flow in Figure 2(f), while the boundary separating the two densities p and q for the other models are skewed more to the left.

F.2 2-D Mixture of Gaussians for Featurized KLIEP/KMM
In this experiment, we construct a synthetic domain adaptation task using 2-D Gaussian mixtures. Our goal is to assess whether our featurized density ratio estimation framework improves the performance of KMM and KLIEP, which operate in input space. We construct our source dataset as Then, we train a logistic regression classifier to distinguish between the two classes using 3 methods: 1) an unweighted logistic regression baseline, 2) reweighted logistic regression with importance weights computed in input space, and 3) reweighted logistic regression with importance weights computed in feature space. The importance weights are learned on a mixture of the source and target datasets.
Results are shown in Table 8.

F.3 Domain Adaptation with the UCI Breast Cancer Dataset
We provide full experimental results of our domain adaptation experiment with the UCI Breast Cancer dataset in Table 9 and Figure 7. Results were computed over 30 runs. We note that our method improves upon KMM for most values of C and achieves the best absolute test error out of all combinations of C with different methods. We also note that KLIEP performs poorly on this task, regardless of the method we use. 0.650 ± 0.0109 0.55 ± 0.177 0.590 ± 0.126 0.586 ± 0.119 Table 9: Comparison between test errors of an unweighted SVM and reweighted x-space and z-space SVMs on classification of the UCI Breast Cancer dataset with the biased class label sampling scheme. Standard errors were computed over 30 runs.

F.4 Omniglot samples from DAGAN
A sample of n = 100 examples synthesized by the trained DAGAN, used for the data augmentation experiments in Section 5, are shown in Figure 8.

F.5 Mutual Information Estimation
In Figure 9, we replicate Figure 4 with additional results from joint training procedures using different values of α = {0.1, 0.5, 0.9} in Equation 1. Specifically, we note that α = 0.9 outperforms all other jointly-trained models, indicating that a greater emphasis on the classification loss term helps for this experiment.

F.6 Samples from MNIST Targeted Generation Task
For each DRE in z-space, DRE in x-space, and unweighted settings and for perc={0.1, 0.25, 0.5, 1.0}, Figures 10, 11, and 12 show n = 100 MAF-generated samples from the diff-background experiments and Figures 13, 14, and 15, show n = 100 MAF-generated samples from the diff-digits experiments.


Faster Multi-Goal Simulation-Based Testing Using DoLesS (Domination with Least Square Approximation)

For cyber-physical systems, finding a set of test cases with the least cost by exploring multiple goals is a complex task. For example, Arrieta et al. reported that state-of-the-art optimizers struggle to find minimal test suites for this task. To better manage this task, we propose DoLesS (Domination with Least Squares Approximation) which uses a domination predicate to sort the space of possible goals to a small number of representative examples. Multi-objective domination then divides these examples into a"best"set and the remaining"rest"set. After that, DoLesS applies an inverted least squares approximation approach to learn a minimal set of tests that can distinguish best from rest in the reduced example space. DoLesS has been tested on four cyber-physical models: a tank flow model; a model of electric car windows; a safety feature of an AC engine; and a continuous PID controller combined with a discrete state machine. Comparing to the recent state-of-the-art paper attempted the same task, DoLesS performs as well or even better as state-of-the-art, while running 80-360 times faster on average (seconds instead of hours). Hence, we recommend DoLesSas a fast method to find minimal test suites for multi-goal cyber-physical systems. For replication purposes, all our code is on-line:https://github.com/hellonull123/Test_Selection_2021.

INTRODUCTION
Simulation models play an important role in many domains. Engineers build such models to simulate complex systems [29]. In the case of cyber-physical systems, these models are sometimes shipped along with the actual device, which means that analysts can now access high-fidelity simulations of their systems. Hence, much of the work on cyber-physical testing focuses on taking full advantage of high-fidelity simulators, prior to live testing [3]. For example, analysts can use the simulators for test suite minimization; i.e. they can explore many tests in the simulator in order to remove tests that do not need to be explored in the real world.
However, using these models for test case minimization can be a very difficult process [3]. These simulation models are built to simulate complex systems such as electronic and physical models. Hence executing these simulation models can be very time consuming [5]. This problem gets even worse for multi-goal problems (e.g. minimizing runtime and maximizing the number of bugs found) since it is necessary to run the models multiple times for different subsets of the goals [5]. For example, Arrieta et al. reported that testing a high fidelity simulation model can take hours to days [3,19,37]. Further, they warned that state-of-the-art multi-goal optimizers (e.g. NSGA-III and MOEA/D) struggle to find minimal test suites for this task.
Recently, Chen [10] and Agrawal et al. [1] reported successes with a variant of optimization called "DUO" (data mining using/usedby optimizers). In this approach, a data mining method firstly divides the problem space, and then an optimizer executes in each small division. Inspired by that DUO approach, in this work, we apply a sorting method on the objective space to divide the multiobjective test suite minimization problem into several smaller partitions. Our DoLesS algorithm (Domination with Least Squares Approximation) applies a domination predicate to sort the example space to a small number of representative data points. Multiobjective domination divides these data points into a "best" set and the remaining "rest" set. After all that, DoLesS applies an inverted least squares approach to learn a minimal set of tests that can distinguish the best from the rest in the reduced example space.
To evaluate our proposed test case selection approach for simulation models, we compare DoLesS with the most recent state-ofthe-art approach produced by Arrieta et al. [3]. In that companion, we ask the following research questions.
RQ1: Can we verify that test case selection for multi-goal cyber-physical systems is a hard problem? Arrieta et al. [3] reported that standard multi-goal optimizers such as NSGA-III and MOEA/D failed in this task. This is an important observation since, if otherwise, there will be no clear motivation for this paper. Accordingly, as a first step, we replicate their results in our experiment.
RQ2: Can DoLesS find better ways to select test cases which result in test suites with higher effectiveness measure (objective) scores? Section 3.1 of this paper reviews five effectiveness measurement metrics which are used to select test cases for cyberphysical systems in the previous study [3]. Our results show that on those five metrics, general performances of DoLesS is better than previous state-of-the-art method.
RQ3: Do selected test cases by DoLesS beat the prior stateof-the-art? Apart from the five effectiveness measurement metrics used by Arrieta et al. [3], two other evaluation scores of interest are (a) reduction in the number of test case and (b) faults detection performance with the reduced test cases. As shown in our result section §5, DoLesS usually performs as well, if not better, than the prior state-of-the-art.
RQ4: Is DoLesS far more efficient than prior state-of-theart in terms of running time? For all the reasons stated above, arXiv:2112.01598v1 [cs.SE] 2 Dec 2021 we need methods that offer faster feedback from models of cyberphysical systems. In this regard, it is significant to note that DoLesS runs 80-360 times faster than the prior state-of-the-art.
Based on the above, we say our novel contributions are: (1) We propose a novel test generation method (DoLesS).
(2) We verify that DoLesS solves a hard problem (test case selection for multi-goal cyber-physical systems). This is a problem that defeats state-of-the-art optimizers (NSGA-III and MOEA/D). (3) We clearly document the value of doing DoLesS. When testing on four cyber-physical models, DoLesS finds test suites as good, as even better, than those found by Arrieta et al. 's approach [3]. Further, DoLesS does so while running 80-360 times faster (seconds instead of hours, mean time). Hence, we recommend DoLesS as a fast method to find minimal test cases for multigoal cyber-physical systems.
The rest of this paper is structured as follow. Section 2 introduces the background and related work in test case selection for simulation-based testing. Section 3 introduces the problem of studying effectiveness measurement metrics in cyber-physical systems and illustrates how they are calculated by mathematical formula. Moreover, multi-objective optimizers and our proposed approach are introduced in this section as well. Section 4 introduces the case studies, performance evaluation metrics, and statistical analysis method used in this study. Section 5 shows our experimental results. Section 6 explores threats to validity and Section 7 makes the summary of our study and states the possible future work.
Based on the above, we can conclude that DoLesS is faster, yet more effective, than prior results since: • DoLesS can handle multiple goals (in our experiment, 5 goals) simultaneously. Hence it does not need to loop the algorithm times (where is the number of subsets for the goals) like the prior state-of-the-art method.
• DoLesS's sorting procedure uses continuous domination to very quickly divide candidates into a very small "best" set (that we can focus on) and a much larger "rest" (that we can mostly ignore). Like much research before us, we argue that continuous domination is more informative than binary domination [38,42,52]. • Chen et al. [10] argues that some SE optimization problems can be solved better by over-sampling than via evolutionary methods. For example, the evolutionary NSGA-II method mutates 100 individuals for 250 generations (these parameters were selected to ensure comparability to the prior study). On the other hand, our DoLesS over-sampling method explores 10,000 individuals for one generation. This result suggests that cyber-physical system testing might be another class of problem that better to be solved via the over-sampling methods which stated by Chen et al [10].

BACKGROUND
A repeated result is that test suites can be minimized (i.e. we can run fewer tests) while still being as effective (or better) than running the larger test suite [2,15,44,45,49]. Note that "effective" can mean different things in different domains, depending on the goals of the testing. For example, at FSE'14, Elbaum et al. [16] reported that Google could find similar number of bugs, but after far fewer tests execution. This was an important result since, at that time, the initial Google test suites were taking weeks to execute. Such long test suite runtimes is detrimental to many agile software practices.
Research has found many test case selection techniques such as DejaVu based, firewall based, dependency based, and specification based techniques [17]. We note that different test suite minimization methods need different kinds of data. For example, in 1995, Binkley et al. proposed a semantic-based method which takes the use of differences and similarities of two consecutive versions to select test cases [8]. Rothermel [5]. Also in 2017, Lachmann et al. [23] did an empirical study on several black-box metrics and made comparisons on their performance in selecting test cases in system testing.
Due to the data requirements, many of the above methods are unsuitable for cyber-physical systems, for two reasons.  (1904). Energy is released (as steam) when the weights spin faster and pull away from the shaft. Thus, increased speed leads to energy release, which slows the system. From Wikipedia.
Firstly, cyber-physical systems are embodied in their environment. Hence, it is not enough to explore static features of (e.g.) the code base. Rather, it is required to test how that code base reacts to its surrounding environments. Hence, using just static information such as (e.g.) code coverage metrics is not recommended for testing cyber-physical systems. Secondly, at least for the systems studied here, cyber-physical systems make extensive use of process control theory. In that theory, the feedback controller is used to compare the value or status of process variables with the desired setpoint. This controller then applies the difference as a control signal to bring the process variable output of the plant to the same value as the set-point (for example, see the steam governor of Figure 1). Hence, for test suite minimization of process control applications, the requirement is data collected from the feedback loops inside the cyber-physical systems. Accordingly, here we use input and output signals in the simulation models instead of execution history or coverage information.
In one of the IST'19 journal paper, Arrieta et al. [3] explored issues associated with test suite minimization by using the data extracted from feedback loops. They noted that feedback loops have anti-patterns; i.e. undesirable features that appear in a time series trace of the output of the system. Figure ?? shows three such features correspondingly from left to right in the red dash rectangle: instability, discontinuity, and growth to infinity. Later in this paper we will mathematically define these anti-patterns.
In all, Arrieta et al. [3] explored seven goals for cyber-physical model testing: maximizing the three anti-patterns that shown in Figure ??, maximizing three other measures of effectiveness, as well as minimizing total execution time. Arrieta et al. used mutation testing to check the validity of their minimized test suite. Mutation based testing is a fault-based testing technique which implements "mutation adequacy score" to assess test suite adequacy by creating mutants [22], and then pruning test cases which cannot distinguish the original model from the mutant. In our study, we use the mutants generated in the experiments from Arrieta et al. [3]. Those mutants were generated with Hanh et al. 's technique [7] and some of the mutants are removed if (a) they are not detected by any test case, (b) they are killed by all test cases, and (c) they are equivalent mutants [34]. Like Arrieta et al., we say a test suite is minimal when it retires as many mutants as a larger suite.
Mutation testing is the inner loop of Arrieta et al. 's process and, in their experiments, they found mutation testing to be an effective technique. The problem area in their work was the outer loop that optimized for seven goals. They found that standard optimizers such as NSGA-II [14] can be ineffective for more than three goals (a result that is echoed by prior work [38]). More recent optimizers like NSGA-III [13] and MOEA/D [50] also failed for this multi-goal task. Later in this paper, we replicate their experiment and strengthen that finding (see RQ1).
To address this optimization failure, they resorted to "pairwise" approach based on NSGA-II. That is, they ran NSGA-II with all 21 subsets of "choose two or three from seven" goals then returned the test suite associated with the run which has the best scores (where the "best" here is measured just on a subset of goals). While definitely an extension to the state-of-the-art, Arrieta et al.'s [3] study had two drawbacks. Firstly, the test cases selected in this way was only the best which measured on a subset of the optimization goals. Secondly, the "pairwise" approach increased optimization time by an order of magnitude, which is a major issue for large simulators, especially when we are running these algorithms 20 times (to check the generalizability of this stochastic process).
Hence in this work, we seek to improve the mutation based test suite minimization method from Arrieta et al. 's study [3]. Like them, we will optimize for the anti-patterns and effectiveness measures seen in process control systems. But unlike that prior work, we will offer methods that simultaneously succeed across many goals (without needing anything like the pairwise heuristic used in Arrieta et al.). Further, we show that all this can be achieved without additional runtime cost.

Testing Simulation models
Cyber-physical system developers often use simulation tool (e.g. Simulink) to build cyber-physical models [12]. For an example Simulink models, see Figure 3. This is a model with two hierarchical levels [3]. A complex model will have far more blocks and operators.
In Simulink models, the inputs and outputs are all signals (here signal means a time series function). This means at each simulation time step Δ , there will be a value in each input and in each output regarding to that time step. For example, if we simulate a model for 5 seconds in real time and the time step Δ is 0.05, then there will be 5/0.05 + 1 = 101 simulation steps, which means each input or output should be a vector of length 101.
Our goal for this study is to select representative test cases from the initial test suite to minimize the test execution time, but not influence the testing performance. Here we can define the test selection problem as follow: Given an initial set of test cases = { 1 , 2 , 3 , · · · , }, we want to find a subset of that set of test cases = { 1 , 2 , · · · , } which can test the model as initial test suite does If we search in the space which contains all the subsets of the initial test suite, then the search space will be very large. For example, with only 100 test cases, there will be 2 100 − 1 possible subsets. Thus, cost-effectively selecting test cases is a significant problem.

EXPERIMENTAL METHODS 3.1 Simulation Effectiveness Metrics
In this study, we implement five out of seven effectiveness measurement metrics which Arrieta et al. [3] used in their study. The first three metrics are also widely used in previous studies [27,28,43], and the forth metric is proposed by Arrieta et al. [3].
Aside: We exclude two of the metrics explored by Arrieta et al. (Input & Output-based test similarity metrics) since these two metrics will always result similar normalized values (0.95-0.99) with different test case selections. Such similar normalized values can affect the performance of multi-objective optimization algorithms. We will introduce the remaining five metrics in the rest of this section.

Test Execution Time.
Total test execution time is the first metric we implement in our study. Wang et al. [43] stated that the number of selected test cases can be treated as the measurement for selecting representative test cases from the initial test suite. However, Arrieta et al. pointed out the problem that each test case has different execution time [3]. In our study, we use the similar Figure 3: Simple example of a Simulink model -Cruise Controller of a car [3].
calculation that Arrieta et al. [3] did in their study to deal with the test execution time. The test execution time is calculated as follow: Define a set of test cases = { 1 , · · · , } that are selected from the initial test suite = { , · · · , }. Let denotes the execution time of the test case in , and denotes the execution time of the test case in . The total test execution time of a set of selected test cases is [3] = ∑ In our study, we want to minimize this metric because the goal of test case selection is to decrease the test execution time.

Discontinuity in Output
Signal. Discontinuity is the second metric we implement in our study. As Matinnejad et al. [27] stated, the discontinuity of the output signal is a short duration pulse in the output signal, which means the output signal increases or decreases to a value in a very short time, and recovers back to normal. If executing a test case causes discontinuity in the output signal, then that test case detects the faulty behavior in the model. Assume we have output signals { 1 , 2 , · · · , }. For discontinuity score of each output signal ( ) where 1 ≤ ≤ , Matinnejad et al. calculated it with [27,28] where = | ( · Δ ) − (( − ) · Δ )|/Δ is the left change rate of step and = | (( + ) · Δ ) − ( · Δ )|/Δ is the right change rate of step .
The discontinuity rate of a set of selected test cases is calculated as follow: With same definition of test case above, let denotes the discontinuity score of the test case in , and denotes the discontinuity score of the test case in . The total discontinuity score of a set of selected test cases is [3] = ∑ is the normalized discontinuity score of the test case among output signals In our study, we want to maximize this metric because the goal is to detect more discontinuity in the output signal.

Instability in Output
Signal. Instability is the third metric we implement in our study. As Matinnejad et al. [27] stated, the instability of the output signal is a duration of quick and frequent oscillations in the output signal, which means the output signal increase and decrease repeatedly in a duration of time. If executing a test case causes instability in the output signal, then that test case detects the undesirable impact on physical process [27]. Assume we have output signals where is the total number of simulation steps and Δ is the time stamp in the simulation model for each step.
The instability score of a set of selected test cases is calculated as follow: With same definition of test cases above, let denotes the instability score of the test case in , and denotes the instability score of the test case in . The total instability score of a set of selected test cases is [3] = ∑ is the normalized instability score of the test case among output signals In our study, we want to maximize this metric because the goal is to detect more instability in the output signal.

Growth to Infinity in Output
Signal. This is the forth metric we implemented in our study. As Matinnejad et al. [28] pointed out, the growth to infinity of the output signal is the phenomenon that the output signal increases or decreases to infinity value. If executing a test case causes growth to infinity in the output signal, then that test case detects the faulty behavior in the model. Assume we have output signals { 1 , 2 , · · · , }. For growth to infinity score of each output signal ( ) where 1 ≤ ≤ , Matinnejad et al. calculated it with [28] where is the total number of simulation steps and Δ is the time stamp in the simulation model for each step. The growth to infinity score of a set of selected test cases is calculated as follow: With same definition above, let denotes the infinity score of the test case in , and denotes the infinity score of the test case in . The total growth to infinity score of a set of selected test cases is [3] = ∑ is the normalized infinity score of the test case among output signals { 1 , 2 , · · · , }.
In our study, we want to maximize this metric because the goal is to detect more growth to infinity situation in the output signal.

Output Minimum and Maximum Difference in Output Signal.
This is the last metric we implemented in our study. Arrieta et al. proposed this metric in their work because the difference between maximum output signal and minimum output signal can indicate the level of how a model is being tested. If executing a test case results in large minimum and maximum difference in the output signal, then that test case can detect more parts in the simulation model. Assume we have output signals { 1 , 2 , · · · , }. For output minimum and maximum difference score of each output signal where is the total number of simulation steps and Δ is the time stamp in the simulation model for each step. The difference of output minimum and maximum of a set of selected test cases is calculated as follow: With same definition of test case above, let denotes the output minimum and maximum difference of the test case in , and denotes the output minimum and maximum difference of the test case in . The total output minimum and maximum difference score of a set of selected test cases is [3] = ∑ is the normalized output minimum and maximum difference score of the test case among output signals { 1 , 2 , · · · , }.
In our study, we want to maximize this metric because the goal is to coverage more parts that can be tested.

Binary vs Continuous Domination.
In the following, all the algorithms use binary domination except for DoLesS that uses continuous domination.
Binary domination decides one individual is better than another if it is better on at least one goal and worse on none. Numerous studies [38,42,52] warn that binary domination is hard to distinguish candidates once the number of goals grows to three or more.
For many-goal problems, Zilter's continuous domination predicate [52] is useful [38,42,52]. Continuous domination judges the domination status of pair of individuals by running a "what-if" query which checks the situation when we jump from one individual to another, and back again. Specifically: • For the forward jump, we compute 1 = − * ( − )/ .
where and are the values on the same index from two individuals, is the number of goals (in our case = 5), and is the weight {-1,1} if we are minimization or maximizing the goal correspondingly. According to Zitler [52], one example is preferred to another if we lost the least jumping to it; i.e. 1 < 2 .
Specifically, in this work, we use this predicate to select better goal sets that (a) minimize test execution time, (b) maximize discontinuity score, (c) maximize instability score, (d) maximize growth to infinity score, and (e) maximize output minimum & maximum difference.

NSGA-II.
NSGA-II is a common evolutionary genetic algorithm [14]. Firstly, it generates an initial set of population as the starter of the entire algorithm. Secondly, these candidates will evolve to offsprings in a series of generations by implementing the crossover and mutation operators with their individual probability. In our reproduction experiment, we use single point crossover with 0.8 probability and bit-flip mutation with 1/ probability ( is the number of test cases). Thirdly, parents for next generation will be selected by selection operator, which utilizes a non-dominated sorting algorithm to select top non-dominated solutions [32]. In the situation where a front needs to be divided because it exceeds the total number of population, NSGA-II uses the crowding distance to split candidates in that group.

NSGA-III.
NSGA-III is an improved NSGA-II algorithm [13]. In NSGA-III, all procedures such as initial population generation, crossover, and mutation are similar to NSGA-II, except selection procedure. In NSGA-III, the selection procedure is applied based on a set of reference points. The reference points are uniformly distributed on the normalized hyper-plane with some division number [13]. After that, each objective point is normalized adaptively and associated with a reference point by calculating its distance to the corresponding reference line. Niche-Preservation operation is then applied to select candidates which will be used in the next generation [13].
3.2.4 MOEA/D. MOEA/D is the first multi-objective optimization algorithm which utilizes decomposition technique [50]. More specifically, MOEA/D explicitly decomposes the problem into multiple sub-problems with less objectives in each subgroup and solves these sub-problems simultaneously [50]. To so, prior to inference, all examples get random weights assigned to their goals.  itself, its local neighborhood will move in the same direction as 0 does. Figure 4 shows the entire framework of our approach. Unlike above evolutionary algorithms, our proposed approach DoLesS (Domination with Least Squares Approximation):

DoLesS.
• Uses continuous domination (defined below, see the first block in Figure 4) to reduce the size of initial large random sets of goals and find a "best" group of representative samples. • For each data entry in the "best" group, DoLesS then uses a least square approximation technique (the third block in Figure 4) to inversely predict the test selection outcomes which can fit the representative sets of goals best. This least square approximation technique is discussed below.
After sorting on the domination score, DoLesS divides data into: • The √ "best" items. In our case study, we randomly generate 10000 initial candidates, hence the first 100 candidates with highest domination score are grouped into the "best" group.
• And the remaining "rest" items.
Here we select the number of final population as 100 with two reasons: (a) 10000 random initial population is large enough to cover a wide range of possible outcomes and (b) to make comparison fair, we select same number of final candidates as previous work [3].
In the data processing stage of Figure 4, we take data from each model (which Arrieta et al. also used in their study [3]), and then processes it into the form of least square approximation structure by combining with the representative goals generated from continuous domination. Table 1(i) shows a simple example of effectiveness measurement data collected from the models. We can find that each test case will have a single score for all effectiveness measurement data ( means the score of test case in effectiveness measure ). The corresponding matrix equation system for the above example is shown in Table 1(ii). This equation shows the linear relationship of test selection outcomes and the final effectiveness measure scores (e.g. the final score of effectiveness measure 1 can be obtained by 1 = 11 · 1 + 12 · 2 + 13 · 3 + 14 · 4 where is the outcome of test selection). In this example, our goal is to find the best outcomes of 1 to 4 which can result 1 to 5 . To summarize the above example, in our approach, we collect effectiveness measurement data for test cases (like Table 1(i)) and want to find the best set of outcomes for 1 to which can get the closest scores to representative goals which are selected by continuous domination. Linear Least Squares Approximation is a method which predicts the best value of a set of unknown variables that fits the relationship between expected and observed sets of data. In general, solving a system of linear equations ( = ) will result no solution or infinite solutions. This always happens when (a) the number of constraints (equations) greater than the number of variables (overdetermined) or (b) the number of variables greater than the number of constraints (underdetermined). The way to find the best approximate solution is called the linear least square approximation. As mentioned in Section §3.2.5, in our study, we have 5 goals and number of test cases (where 100 ≤ ≤ 150). Thus, our equation system contains 5 equations (5 goals) and variables (where must greater than 5). In this case, finding possible selections of test cases becomes a underdetermined least square approximation.
In the formulation = , where is an outcome vector of test cases, we want to predict the value (0/1) for each entry of . The final outcome of is a vector of float numbers (ranged from 0 to 1 by controller) which indicates lower effect to the final score with coefficient → higher effect to the final score with coefficient from 0 → 1. Since test selection outcomes can only have 0 (discard that test) and 1 (select that test), we use the threshold of 0.5 to indicate higher probability or lower probability. A value < 0.5 means higher chance to be 0 and a value > 0.5 means higher chance to be 1. For each representative candidate found by continuous domination, DoLesS finds the test selection which can get the closest score to that candidate. Although there exists delta between original ideal scores and truth scores generated by predicted results because of the approximation procedure, our results show that least square approximation can find adequate test cases which perform as well or better as the previous state-of-the-art. Our implementation of above step uses scipy.optimize, a python library, and uses the function called lsq_linear, which solves the above problem by using either dense QR decomposition technique or Singular Value Decomposition technique. Finally, in the Evaluation stage of Figure 4, DoLesS selects the Pareto front set 1 from the final population as Arrieta et al. did in their study [3]. All evaluations are made through 20 repeats for each of algorithm.

EXPERIMENTAL SETUP 4.1 Case Studies
We use four Cyber-physical system (CPS) models to evaluate our proposed approach. These four models come from the previous state-of-the-art study [3]. We implement the test cases and mutants that Arrieta et al. generated from these four models 2 . The summary of number of initial test cases and number of mutants are shown in Table 2. In that table: • Two Tanks project is a model that simulate the incoming and outgoing flows of the tanks [30]; • CW project is a model that simulate the electrics and mechanics of four car windows [3]; • AC Engine project is a model that simulate some safety functionalities in the AC engine [4,6]; • and EMB project simulates the software model controller which includes a continuous PID controller and a discrete state machine [27].
At first glance, the case studies in Table 2 may appear to contain very small test cases. But appearances can be deceiving; e.g. the number of input signals is a poor measure of the internal complexity of a cyber-physical system. As shown in our RQ1 results, the systems of Table 2 are so complex that, for the purposes of test suite minimization, they defeated state-of-the-art optimizers (NSGA-III and MOEA/D).

Performance Criteria
To evaluate the selected test cases, we use two evaluation metrics from prior work [3]. These two evaluation metrics are (a) normalized test execution time and (b) mutant detection score. Previous study used these two metrics to calculate the hypervolumne indicator and average weighted sum of mutation score and normalized test execution time [3], while in our study, we directly compare the performance of algorithms in these two metrics. Normalized test execution time (TET-): Our goal for selecting test cases from the initial test suite is to speed up the testing process. Thus, test execution time is a very important indicator which can indicate whether selected test cases can significantly reduce the cost of testing. In this study we want to minimize this value since, as discussed in our introduction, the whole point of this paper is to reduce the time required for testing cyber-physical systems.
Mutant detection score (MS+): If a set of selected test cases can significantly reduce the test execution time, but cannot detect most of the mutants, then such a selection is a bad choice. Our goal for selecting test cases is detecting as much as mutants when minimizing the test execution time. Therefore, mutant detection score becomes another important evaluation metric. In this study we want to maximize this value since higher value means the test suite is better since it can detect more mutants.
That is, a good test case selection approach can both (a) minimize the test execution time and (b) maximize the mutant detection score.

Statistical Analysis
In our study, we record the value of above two evaluation metrics in 20 repeats. To compare the total performance of different algorithms, we implement A Scott-Knott analysis [31]. The Scott-Knott analysis can sort the candidates by their values, and assign candidates to different ranks if the values of candidate at position is significantly  [3]. The metric with "-" means less is better while "+" means more is better. The light gray cell in each project means that approach wins others significantly (as computed by the statistical method in §4.3).  Table 4: RQ2 results: Scores of five effectiveness measurement metrics calculated by sets of selected test cases. All entries report the median score of 20 repeats. In the title row, the metric with "-" means we want to minimize that metric and the metric with "+" means we want to maximize that metric. The light gray cell in each project means that approach wins over another approach significantly (as computed by the statistical method of §4.3) in that metric. Last column counts the number of wins for each approach. different (by more than a small effect size) to the values of candidate at position − 1 [24]. More precisely, Scott-Knott sorts the candidates by their median scores (and in our study, the candidates are the test case selection approaches). Scott-Knott method will split the sorted candidates into two sub-lists which maximize the expected value of differences in the observed performances before and after division [41]. After that, Scott-Knott will declare the one of the split as the best split. The best split should maximize the difference (Δ) in the expected mean value before and after the split [40,46]: where | |, | 1 |, and | 2 | are size of list , 1 , and 2 . , 1 , and 2 are mean value of list , 1 , and 2 . After the best split, Scott-Knott then implements some statistical hypothesis tests to check the division. If two items 1 and 2 after division differ significantly by applying hypothesis test , then such division is defined as a "useful" division. Scott-Knott will run recursively on each half of the best division until no division can be made. In our study, we use cliff's delta non-parametric effect size measure as the hypothesis test. Cliff's delta quantifies the number of difference between two lists of observations beyond p-values interpolation [46]. The division passes the hypothesis test if it is not a "small" effect ( ≥ 0.147). The cliff's delta non-parametric effect size test explores two list and with size | | and | |: Cliff's delta estimates the probability that a value in the list is greater than a value in the list , minus the reverse probability [25] in the above formula. This hypothesis test and its effect size is supported by Hess and Kromery [21].

RESULTS
Returning now to the research questions offered in the introduction, we offer the following results. RQ1: Can We Verify that Test Case Selection for Multigoal Cyber-physical Systems is a Hard Problem? To answer RQ1, we replicate Arrieta et al.'s experiment [3]. Table 3 shows our replication results with 20 repeats (with different random number seeds). Two algorithms differ significantly if they separate in different ranks in the Scott-Knott analysis.
As seen in Table 3, we can found the approach with NSGA-II that Arrieta et al. proposed [3] performs better than other multigoal optimizers (MOEA/D and NSGA-III). Specifically, NSGA-II has higher performance in both test execution time and mutation score in four case studies. However, even though NSGA-II beats the other methods, we cannot sanction its use. NSGA-II has all the problems discussed in §2. Specifically, NSGA-II can only handle pairs of goals. Hence, it has to be re-run multiple times to explore all pairs of five goals. As shown below, we can achieve better results, orders of magnitude faster. Therefore, we can answer RQ1 as follow: Test case selection for multi-goal cyber-physical models is a hard problem that cannot be solved by merely applying, off-the-shelf, the latest optimizer technology.
These RQ1 results motivate the rest of this paper where we develop a fast approach, which can handle multiple goals simultaneously for the cyber-physical systems.
RQ2: Can DoLesS Find Better Ways to Select Test Cases which Result in Test Suites with Higher Effectiveness Measure (Objective) Scores? To answer RQ2, we calculate the scores of effectiveness measures (test execution time, discontinuity, growth to infinity, instability, minimum and maximum difference) after we generate the subsets of selected test cases. Table 4 shows our simulation results. For each project, we use Scott-Knott statistical analysis to compare the performance across 20 repeats. Table 4 reports the median scores for each effectiveness measurement metric. To visualize the final results, we mark the winning approach in each metric by light gray , and count the number of wins in the last column.
As seen in Table 4, DoLesS wins in three out of four projects (in Twotanks, CW, and EMB). Moreover, in Twotanks and EMB, our proposed approach results higher scores on most of the effectiveness measurement metrics while previous approach only wins in one (DoLesS wins in all effectiveness measurement metrics in Twotanks and wins in 4 out of 5 effectiveness measurement metrics in EMB). These outcomes can indicate that in most of the cases, our proposed approach gets significant improvement than previous approach in these five metrics. Moreover, in half projects, the state-of-the-art approach concentrates on optimizing the Time metric and Instability metric (because of the algorithm design). However, by using DoLesS, we can find that most of the "goals" are equally optimized. Hence, DoLesS can handle multiple goals all the time while state-of-the-art method can only handle one or two goals in some cases.
By summarizing above findings, we answer RQ2 as follow: DoLesS can find better test selections than state-of-the-art approach in terms of five effectiveness measurement metrics.
RQ3: Do selected test cases by DoLesS beat the prior stateof-the-art? To answer RQ3, we compare scores of TET-(normalized test execution time) and MS+ (mutant detection score) between the prior state-of-the-art and DoLesS. It is important to have higher performance on these two evaluation metrics since they directly indicate whether a test case selection is good or not. Table 5 shows our simulation results (note that TET -test execution time & MS -Mutation Score). For each method in each project, we repeat experiment 20 times and calculate the value of two evaluation metrics for each repeat. To obtain the final conclusion, we implement Scott-Knott statistical method to check if our approach significantly differ to state-of-the-art approach in each metric. The light gray cells mark the winning method (in the first rank) resulted from the Scott-Knott test. Moreover, we count the number of higher performance in these two evaluation metrics for each algorithm and record the number of wins in the last column.
As seen in Table 5, DoLesS gets better performance in both two evaluation metrics in two out of four projects (ACEngine & EMB). Moreover, in Twotanks, DoLesS and state-of-the-art method achieve the same performance (both two evaluation metrics are tied in the first rank). In the CW project, DoLesS has higher performance in minimizing test execution time when state-of-the-art method gets a little bit higher mutation score than DoLesS. Taking above comparisons together, we can conclude that test cases selected by DoLesS can achieve similar or better performance in minimizing test execution time and detecting more mutants in all projects.
By summarizing above findings, we answer RQ3 as follow:  In all projects, comparing to the state-of-the-art, DoLesS can get similar or better performance on minimizing the execution time of selected test cases while keeping to detect most of the mutants.
RQ4: Is DoLesS far more efficient than prior state-of-theart in terms of running time? To answer RQ4, we count the execution time for both algorithms during the experiment. To make comparison fair enough, we run both two algorithms on the same 64-bit Windows 10 machine with a 4.2 GHz 8-core Intel Core i7 processor and 16 GB of RAM. Moreover, when running experiments, we make sure no huge process is starting or ending in our machine. Table 6 shows the recorded runtime for each project. For each method, we repeat experiments 20 times and record the total runtime. The light gray cells mark the fastest approach.
As seen in Table 6, in both four projects, DoLesS runs significantly faster (80-360 times faster) than the previous method. By analyzing our proposed algorithm and state-of-the-art approach, we find previous approach implemented NSGA-II as their multiobjective optimization, which designed for 2 or 3 objectives [32]. To handle this issue, Arrieta et al. group objectives into 21 different combinations with two or three objectives in each combination, and select one of the best combinations by repeating their approach in those 21 groups [3]. However, in our approach, we just need continuous domination to find "ideal goals" and approximate corresponding test cases inversely.
By summarizing above findings, we answer RQ4 as follow: In both four projects, DoLesS run significantly faster (80-360 times) than the previous method. In other words, DoLesS is far more efficient than state-of-the-art approach.
Even though our current empirical results can only boost a speed of (up to) 300 times faster, we can make a theoretical case that, if the number of goals increases, our technique would be even more comparatively faster (please note that in the following counts, test execution time is the metric that must be in every combination): • 5 goals will result 10 different combinations of 2 or 3 objectives.
• 7 goals will result 21 different combinations of 2 or 3 objectives.
• 9 goals will result 36 different combinations of 2 or 3 objectives.
• The above pattern shows that with more goals being utilized, the number of repeats for state-of-the-art NSGA-II approach increases exponentially. However, our approach can handle multiple goals simultaneously in one time. This can show the efficiency of our approach.

THREATS TO VALIDITY
This section discusses issues raised by Feldt et al. [18]. Construct validity: The construct validity threat mainly exists in the parameter settings of algorithms. For example, in our replication experiment, we use one point crossover with 0.8 crossover probability and bitflip mutation with 1/(number of variables) mutate probability as prior studies did in order to keep consistent. For another example, in least square approximation, we use 0.5 threshold to indicate whether a test case has large probability to be chose or not. Moreover, we use the default setting of Python least square solver in our algorithm. Changing these parameters can result differences in selecting test cases. Therefore, our observation may differ when different parameters are used. We would consider hyper-parameter tuning [39,40] in future work to mitigate this threat.
Conclusion validity: The conclusion validity threat in this study is related to the random variations of our algorithm. To reduce the effect caused by this threat, we repeat all experiments 20 times in the same machine. Moreover, we apply Scott-Knott statistical test to compare if the outcomes of our proposed approach and the previous methods differ significantly.
Internal validity: Internal validity focuses on the correctness of the treatment caused the outcome. In this study, we constraint our simulations to the same data set. Moreover, we evaluate our approach and the previous approach in the same workflow. Another internal validity threat can refer to the mutants generated from the projects. To mitigate this threat, we use the same mutants that Arrieta et al. [3] used in their study, which they have removed duplicated mutants.
External validity: External validity concerns the application of our algorithm in other problems. In this study, we generate our conclusion from four real-world Simulink cyber-physical systems. When applying our method into other case studies, these concerns may raised: (a) DoLesS may not be applicable in projects which cannot obtain the effectiveness measurement data; (b) DoLesS may needs modification for those projects which effectiveness measurement data and test cases are not in the linear relationship. For those projects, we would consider non-linear least square approximation as future work to mitigate this threat; (c) Other algorithms may perform better if other information is utilized.

CONCLUSION & FUTURE WORK
Finding representative test cases from the initial test suite is an important task in simulation based model. Better test case selection methods can not only reduce the test execution time in the future testing, but also maintain the same testing performance as usual.
In other word, a good test case selection approach can (a) minimize the test execution time and (b) maximize the mutation score.
Previous literature by Arrieta et al. [3] has shown the great success on using NSGA-II as the multi-objective optimization method to select representative test cases. However, their design has a deficiency where they need to evaluate 21 combinations first to select the best subset. Moreover, since NSGA-II is a randomized algorithm, repeats are necessary during the experiment. Therefore, their approach will execute NSGA-II 420 times with 20 repeats.
In this study, we address this deficiency by selecting test cases from all effectiveness measurement metrics. To do that, we use a very fast approach -continuous domination to select representative goals. Moreover, we make a better use of the linear relationship between test cases and goals to find the best test selections correspond to the representative goals (by linear least square approximation). Our experimental results show that our proposed approach DoLesS can (a) achieve better scores in terms of five effectiveness measurement metrics, (b) achieve better scores in terms of two evaluation metrics, and (c) solve test selection problem in a very fast speed (80-360 times faster than the state-of-the-art method) for cyber-physical models.
We conjecture that our method would be a better candidate for scaling to large systems than the method proposed by Arrieta et al. [3]. To see that, consider the following scenario: To successfully perform test case selection on selected cyber-physical case studies, Arrieta et al.'s approach required 3-5 hours execution time. Now imagine in some higher complexity simulation models (e.g. drone simulation models) with dozens more of test cases in the initial test suite, and these models have more signal processing criteria in I/O signals, both evaluation time and objectives are increased. In such scenario, the execution time of running NSGA-II in all subsets of objectives will exponentially increase as we mentioned at the end of RQ4. Moreover, such models (e.g. drone simulation models) requires far more fast feedback than usual cyber-physical models. Due to the above reasons, the ideal test case selection approach for complex simulation models needs to handle multiple goals (more than 4 goals) in the same time and perform selection in a very short execution time for the fast feedback. For that task, we recommend DoLesS.
As to further work, apart from extending this exploration of feedback loop anti-patterns, we conjecture that our methods could be useful for other multi-objective reasoning tasks. Standard practice in this area is to mutate large populations across a Pareto frontier. This has certainly been a fruitful research agenda [26,33,48,51]. But perhaps the testing community could reason about more goals, faster, if used our domination methods and least squares methods to "reason backwards" from goal space to decision space. Hence, future works can be conducted with • Finding more simulation projects which can strength our approach. • Developing more effectiveness measurement metrics which can better indicate representative test cases. • Adjusting our approach based on the testing scenarios of different projects. Moreover, in some cases, a new test case selection approach is needed for those projects.