NEJSDS The New England Journal of Statistics in Data Science 2693-71662693-7166 New England Statistical Society NEJSDS34 10.51387/23-NEJSDS34 Methodology Article Machine Learning and Data Mining Construction of Supersaturated Designs with Small Coherence for Variable Selection QiYouranyqi28@wisc.edu ChienPeterpeter.chien@wisc.edu∗ Department of Statistics, University of Wisconsin-Madison, 1300 University Ave, Madison, WI 53706, USA. E-mail address: yqi28@wisc.edu Department of Statistics, University of Wisconsin-Madison, 1300 University Ave, Madison, WI 53706, USA. E-mail address: peter.chien@wisc.edu Corresponding author. 2023562023133233331942023 © 2023 New England Statistical Society2023 Open access article under the CC BY license.

The supersaturated design is often used to discover important factors in an experiment with a large number of factors and a small number of runs. We propose a method for constructing supersaturated designs with small coherence. Such designs are useful for variable selection methods such as the Lasso. Examples are provided to illustrate the proposed method.

 Design of Experiments Supersaturated Design Unbalanced Design Coherence Linear Model Variable Selection Introduction

Modern experiments can involve a large number of factors. Assuming the effect sparsity principle [3], an experiment with high-dimensional inputs only has a small number of significant factors. Two-level supersaturated designs are often used to identify important factors in screening experiments. Data from such a design can be modeled by a linear model [33] given by y = X β + ϵ , where X is the n × p design matrix at two levels − 1 and + 1, denoted by − and +, y is the response vector of the n observations, β is the regression vector of p coefficients, the error vector ϵ follows the multivariate normal distribution N n ( 0 , σ 2 I ) and I is the n × n identity matrix. In X, the first column is the intercept with all +’s and the remaining columns represent the p − 1 factors. We call X balanced if all its factors consist of an equal number of − and +. If p = n, X is a saturated design since all the degrees of freedom are used to estimate β. If p > n, X is a supersaturated design since there are not enough degrees of freedom to estimate all components of β. Construction of supersaturated designs that can accommodate a large number of factors relative to the number of runs has been gaining more and more attention [5, 16, 30, 26]. Popular criteria for comparing supersaturated designs include minimax [2], E ( s 2 ) [18, 27] and mean square correlation [8].

It is increasingly common to use modern variable selection methods to analyze data from a supersaturated design and identify important factors [18, 32, 17]. For example, when the number of variables is greater than the number of observations, Osborne, Presnell and Turlach [22] advised to use the Lasso [28] for an initial selection, and then conduct a best subset selection (or similar method) on the variables selected by the Lasso. We focus on discussing and demonstrating our method for the Lasso given its popularity. The proposed designs also work with other penalized regression methods. For variable selection methods such as the Lasso, Zhao and Yu [34] showed that their performance depends more critically on the worst case column correlation than the average column correlation of the design matrix. More details about this point can be found in Appendix B.

Motivated by the importance of controlling the worst case column correlation, we propose a new method to construct supersaturated designs according to the coherence criterion [7], which measures the worst case column correlation. The coherence of an n × p matrix X is μ ( X ) = max 1 ≤ i < j ≤ p | ⟨ x i , x j ⟩ | ‖ x i ‖ ‖ x j ‖ , where the subscript i denotes the ith column, ⟨ · , · ⟩ is the dot product and ‖ · ‖ is the L 2 -norm. Whenever there is no confusion, hereinafter we will drop the symbol X in μ ( X ). Clearly, 0 ≤ μ ≤ 1. The smaller μ is, the closer the matrix is to orthogonal. If the columns of X are centered, μ is reduced to the maximum absolute column correlation of X, which is equivalent to the correlation criterion used in [19]. Note that μ ( X ) = max 1 ≤ i < j ≤ p | x i ⊤ x j | n = s max n , where s max = max 1 ≤ i < j ≤ p | x i ⊤ x j | is the value used in the minimax criterion [2]. Different from s max , μ ( X ) reflects the sample size n in its definition. Because the sample size matters in many data analysis procedures, we will construct designs under μ ( X ) instead of s max (we will also compare designs with different sample sizes in Section 3.1).

Our construction method allows some columns of the design matrix to be unbalanced. Factor balance in a supersaturated design guarantees an accurate estimate of the intercept but unbalance provides other significant benefits in constructing designs. Being unbalanced provides flexibility in controlling the column correlations between every pair of the main effects to achieve a lower average column correlation [29]. Sacrificing the precision of the estimate of the intercept can better estimate the main effects. In addition, allowing unbalance enables us to construct designs with more flexible sizes.

The remainder of the article will unfold as follows. In Section 2, we detail our construction method. In Section 3, we use several examples to demonstrate the advantage of the constructed designs against some benchmark designs for a variable selection problem. We conclude the paper and provide some discussions in Section 4.

 The Construction Method Notation and Definitions

Throughout, “#Balance” denotes the number of balanced columns of a design. Let 1 n denote the n-dimensional unit vector with + 1’s. Let I denote the identity matrix. For a matrix A with entries in { − 1 , + 1 } and an even number of rows, let A ∗ denote a matrix obtained by changing the signs of entries in all odd rows and A ∗ ∗ denote a matrix obtained by changing the signs of entries in all even rows, where ∗ and ∗ ∗ are two matrix operators. Note that ( − A ) ∗ ∗ = A ∗ . An n × n Hadamard matrix H is a matrix with entries in { − 1 , + 1 } and H ⊤ H = n I, where n can only be 1, 2 or a multiple of 4.

 Construction Steps

Our construction method expands a 6 m × p supersaturated design D 0 with μ ≤ 1 / 3 to a 12 m × 4 p design with μ ≤ 1 / 3, where m , p are positive integers with 6 m < p. Partition D 0 as D 0 = U L , where U and L have 4 m and 2 m rows, respectively. The proposed method has three steps.

Step 1: Use two copies of D 0 to obtain a matrix

Step 2: Expand D 1 to obtain a matrix

Step 3: Expand D 2 to obtain a matrix

If D 1 is viewed as a block matrix with four rows, Step 2 applies the operator ∗ to the second row and changes the signs of all entries in the fourth row to obtain D 2 . If D 2 is viewed as a block matrix with four rows and two columns, Step 3 applies the operator ∗ ∗ to the first and fourth rows and applies the operator ∗ to the third row to obtain D 3 . It can be proved that the coherence of D 1 , D 2 and D 3 is no greater than 1 / 3 by Lemmas 1, 2, 3 and Theorem 1 in Appendix A.

Since the design D 3 from the construction is a 12 m × 4 p design with μ ≤ 1 / 3, repeating the above procedure multiple times with D 0 being D 3 produces a class of designs with ( 2 k × 6 m ) rows, ( 4 k × p ) columns and μ ≤ 1 / 3 for every positive integer k.

The only requirement for D 0 is to be a 6 m × p two-level supersaturated design with μ ≤ 1 / 3. For a given pair of m and p, different choices of D 0 yield different forms of D 3 but all of them have guaranteed small coherence. As suggested by Chen and Lin [6] and Liu, Ruan and Dean [20], a supersaturated design with coherence of 1 / 3 can identify most of the important factors.

 Examples

We now provide several examples for the proposed method with different choices of D 0 . In the following examples, we show the first several designs constructed by the proposed method and choose D 0 to have as many columns as possible.

Let D 0 be a 6 × 16 design given by + + + + + + + + + + + + + + + + + + + + − − − − − − − − + + + + + + − − + + − − + + − − + + − − + + − − + + − − − − + + − − + + + − + − + − + − + − + − + − + − + − + − − + − + + − + − − + − + , where − and + denote − 1 and + 1, respectively, and its coherence is 1 / 3. The coherence of the above design attains the lower bound derived in [31]. The sizes, coherence values and numbers of the balanced columns of the first several designs constructed by the proposed method are given in Table 1.

Designs Constructed from a 6 × 16 D 0 .

Size μ #Balance
12 × 64 x x 1/3 24
24 × 256 x 1/3 168
48 × 1024 1/3 840
96 × 4096 1/3 3720
192 × 16384 1/3 15624

Let

where D 01 is the 24 × 23 Plackett and Burman design [23], D 02 is formed by taking the two-order interaction terms of D 01 and D 03 is formed by taking the three-order interaction terms of D 01 . Then D 0 is a 24 × 2048 design with μ = 1 / 3. The first 277 columns of this design form the design in [32]. The sizes, coherence values and numbers of the balanced columns of the first several designs constructed by the proposed method are given in Table 2.

Designs Constructed from a 24 × 2048 D 0 .

Size μ #Balance
48 × 8192 x 1/3 2968
96 × 32768 1/3 18616
192 × 131072 1/3 99064

Let D 0 be the 30 × 59 design from Lin [18], where an additional column 1 30 is added as the first column. This design has μ = 0.2. The sizes, coherence values and numbers of the balanced columns of the first several designs constructed by the proposed method are given in Table 3. This example illustrates that our method can use a design D 0 with μ < 1 / 3.

Designs Constructed from a 30 × 59 D 0 .

Size μ #Balance
60 × 236 x 1/3 101
120 × 944 x x 1/3 550
240 × 3776 x 1/3 2864
480 × 15104 1/3 13156
960 × 60416 1/3 56396

Whenever necessary, one can select some columns of the designs constructed for use. Selecting a subset of columns could remain coherence the same or further decrease it, but will never increase coherence. For example, for a given p, first, one can select as many balanced columns as possible. If there are no sufficient balanced columns, then one can select some unbalanced columns to achieve p columns. In addition, one may also select the columns according to an additional criterion such as E ( s 2 ).

 Generalization

We now generalize the proposed method to obtain a supersaturated design with μ smaller than 1 / 3. Recall that our method partitions D 0 into two parts. The upper part U has 6 m r rows and the lower part L has 6 m ( 1 − r ) rows with a partition ratio r = 2 / 3. We generalize the original construction by using a different partition ratio and stopping at Step 2. Suppose that D 0 is an n × p two-level supersaturated design with μ ≤ t / n for t ∈ { 2 , 4 , … , n / 2 } and an even n. Partition D 0 as in (2.1), where U has 2 t rows and L has n − 2 t rows with r = 2 t / n. Define D 2 as in (2.2). Then D 2 is the supersaturated design constructed by our generalized method. It can be proved that the coherence of D 2 constructed here is no greater than t / n by Theorem 2 in Appendix A.

The generalization indicates the possibility of expanding any two-level supersaturated design with μ ≤ 1 / 2 and an even number of runs while retaining its coherence. For the original construction, coherence must be no greater than 1 / 3 and the initial design must have 6 m runs but the initial design with p columns can be expanded to a design with 4 p columns and μ ≤ 1 / 3. The generalization only requires that the coherence is no greater than 1 / 2 and the initial design has an even number of runs but the initial design with p columns can only be expanded to a design with 2 p columns and the same coherence. In addition, the number of balanced columns can be fewer than the number of runs for the designs constructed by the generalization. Here is an example.

Let D 0 be the 30 × 59 design in Example 3 with μ = 0.2. The sizes, coherence values and numbers of the balanced columns of the first several designs constructed by the generalization are given in Table 4.

Designs Constructed from a 30 × 59 D 0 by the Generalization.

Size μ #Balance
60 × 118 0.2 76
120 × 236 . . 0.2 112
240 × 472 x 0.2 184
. 480 × 944 x . 0.2 328
960 × 1888 0.2 616
Simulation Study

In this section, we compare the proposed designs with four popular classes of supersaturated designs: Lin’s designs [18], Wu’s designs [32], the Bayesian D-optimal supersaturated designs (Jones, Lin, and Nachtsheim 2008) and the U E ( s 2 )-optimal designs [15] for the Lasso problem by simulations. These designs are denoted by “LIN”, “WU”, “BAYES” and “JM” respectively and our proposed designs are denoted by “Proposed”.

Lin’s designs are constructed from the Plackett and Burman design [23]. According to Bulutoglu and Cheng [4], Lin’s designs are E ( s 2 )-optimal among all the balanced two-level supersaturated designs. Wu’s designs can be obtained by appending interaction columns to a Hadamard matrix. Bulutoglu and Cheng [4] showed that in certain cases, Wu’s designs are E ( s 2 )-optimal among all the balanced two-level supersaturated designs. The Bayesian D-optimal supersaturated design is obtained by the coordinate exchange algorithm [9, 14]. Any n rows of a p × p Hadamard matrix form a type of U E ( s 2 )-optimal design [15, 5].

We simulate data from the linear model in (1.1) with different designs and active coefficients, where the error vector ϵ follows the standard multivariate normal distribution N n ( 0 , I ). We use the Lasso to select the active factors and fix the intercept term β 1 as active. We repeat the above data generation and variable selection procedure N = 300 times. We use the “cv.glmnet” function [12] in the R software [24] and calculate β ˆ with “s=lambda.min”.

We conduct simulations using the eight active coefficients settings in Table 5 and denote them as Case 1, … , Case 8, respectively. We use “#Active” to denote the number of active factors. For each group of active coefficients, we consider both the case where they take large values, and the case where they take small values. For example, for β 1 , β 21 and β 23 , we consider both the case where they take values of 1, 5 and 10, and the case where they take values of 0.1, 1 and 1.5. We let the regression coefficients of all inactive factors be zero and do not write them in Table 5. We neither use a single setting of active coefficients for all the comparisons, nor use all settings in Table 5 for each comparison, because of the following reasons. First, the designs under consideration have different numbers of columns and some designs do not have enough columns for a given setting, which makes the setting not applicable. Second, in practice, the more factors we are considering, the more active factors there tends to be, so we use settings with more active coefficients for designs with large numbers of columns, and settings with fewer active coefficients for designs with small numbers of columns.

Active Coefficients Settings.

Case #Active Active Coefficients and Their Values
1 2 β 1 00 = 1 .0 x β 21 0 = 5 .0 x β 23 0 = 10 .
2 2 β 1 00 = 0.1 x β 21 0 = 1 .0 x β 23 0 = 1.5
3 4 β 1 00 = 1 .0 x β 6 00 = 5 .0 x β 8 00 = 10 . x β 9 00 = 10 . x β 26 0 = 15 .
4 4 β 1 00 = 0.1 x β 6 00 = 0.6 x β 8 00 = 1.7 x β 9 00 = 0.9 x β 26 0 = 1 .0
5 4 β 1 00 = 1 .0 x β 19 0 = 5 .0 x β 59 0 = 14 . x β 143 = 10 . x β 214 = 13 .
6 4 β 1 00 = 0.1 x β 19 0 = 1.3 x β 59 0 = 0.9 x β 143 = 1 .0 x β 214 = 1.2
7 6 β 1 00 = 1 .0 x β 38 0 = 5 .0 x β 111 = 5 .0 x β 122 = 9 .0 x β 123 = 9 .0 x β 147 = 9 .0 x β 151 = 11 .
8 6 β 1 00 = 0.1 x β 38 0 = 1.5 x β 111 = 1.5 x β 122 = 1.9 x β 123 = 1.9 x β 147 = 1.9 x β 151 = 2.8

Comparison with Lin’s Design.

Case #Active Design Size μ E ( s 2 ) #Balance AFDR AMR MSE EME
1 2 Proposed 12 × 27 0.33 9.70 24 0.49 0.00 0.64 7.22
LIN 14 × 27 0.43 7.84 26 0.59 0.00 0.93 9.18
2 2 Proposed 12 × 27 0.33 9.70 24 0.51 0.02 0.74 9.15
LIN 14 × 27 0.43 7.84 26 0.59 0.08 1.13 12.51
3 4 Proposed 12 × 27 0.33 9.70 24 0.25 0.18 73.58 537.69
LIN 14 × 27 0.43 7.84 26 0.40 0.28 230.66 794.73
4 4 Proposed 12 × 27 0.33 9.70 24 0.49 0.29 2.15 18.52
LIN 14 × 27 0.43 7.84 26 0.64 0.53 4.19 19.79

For each design in our simulations, we will show its size, coherence, and number of balanced columns. Although E ( s 2 ) is not our focused criterion, we will also show the E ( s 2 ) of each design for reference. We use the following four criteria to compare variable selection and model fitting accuracy of each design:

Average False Discovery Rate (AFDR) ∑ i = 1 N FDR i / N, where FDR i = the number of falsely discovered coefficients the number of discovered coefficients if the discovered model is not a null model, and FDR i = 1 if the discovered model is a null model.

Average Miss Rate (AMR) ∑ i = 1 N MR i / N, where MR i = the number of undiscovered active coefficients the number of active coefficients .

Mean Squared Error (MSE) ∑ i = 1 N ‖ β − β ˆ ( i ) ‖ 2 / N to estimate E ( ‖ β − β ˆ ‖ 2 ).

Expected Model Error (EME) ∑ i = 1 N ‖ X β − X β ˆ ( i ) ‖ 2 / N to estimate E ( ‖ X β − X β ˆ ‖ 2 ).

 Smaller values of AFDR, AMR, MSE and EME are desirable. Note that AFDR and AMR are computed based on the number of discovered or active coefficients, not the number of discovered or active factors.

Comparison with Lin’s Design

We obtain the 12 × 27 proposed design by selecting the intercept, the 24 balanced columns and the first two unbalanced columns of the 12 × 64 design in Example 1. We use Lin’s 14 × 27 design directly. The result is shown in Table 6, where and in other tables below, the smaller values of the four criteria are in boldface. Table 6 indicates that the proposed design outperforms Lin’s design in terms of variable selection and parameter estimation with the Lasso. In addition, the proposed design has fewer runs than Lin’s design.

To verify that this result is not due to some better property of the selected active factors in one design versus another, we also conducted a follow-up comparison with randomly selected active factors. More specifically, for a given setting in Table 5, we randomly select an #Active number of active factors with the same values given in the setting, and repeat 10 times to get 10 sets of active factors. Then, for each set of randomly selected active factors, we repeat the above data generation and variable selection procedure N = 30 times for the proposed design and Lin’s design. The final criteria are averaged over the criteria from the 10 sets of active factors, and they can be found in Table 7. According to Table 7, the proposed design still generally outperforms Lin’s design, even with randomly selected active factors. To save computational cost and keep our paper concise, we will omit such a verification for the following comparisons.

Comparison with Lin’s Design with Randomly Selected Active Factors.

Case #Active Design Size μ E ( s 2 ) #Balance AFDR AMR MSE EME
1 2 Proposed 12 × 27 0.33 9.70 24 0.61 0.00 5.01 9.02
LIN 14 × 27 0.43 7.84 26 0.61 0.00 1.59 11.58
2 2 Proposed 12 × 27 0.33 9.70 24 0.58 0.09 1.34 10.83
LIN 14 × 27 0.43 7.84 26 0.59 0.14 1.36 12.93
3 4 Proposed 12 × 27 0.33 9.70 24 0.47 0.00 51.62 27.23
LIN 14 × 27 0.43 7.84 26 0.41 0.04 33.89 135.42
4 4 Proposed 12 × 27 0.33 9.70 24 0.55 0.15 2.11 13.23
LIN 14 × 27 0.43 7.84 26 0.45 0.27 2.13 18.56

Comparison with Wu’s Designs.

Case #Active Design Size μ E ( s 2 ) #Balance AFDR AMR MSE EME
1 2 Proposed 12 × 64 0.33 9.90 24 0.65 0.00 0.95 8.41
WU 12 × 64 0.33 11.06 63 0.54 0.00 0.89 9.03
2 2 Proposed 12 × 64 0.33 9.90 24 0.66 0.05 1.20 12.46
WU 12 × 64 0.33 11.06 63 0.56 0.04 0.98 11.22
__________________ _____________________ _ ________________________ ________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________ ________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________ ________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________ ________________________ ___________________ _________________ ________________ _________________
3 4 Proposed 24 × 256 0.33 21.84 168 0.41 0.00 0.96 19.96
WU 24 × 256 0.33 23.00 255 0.65 0.00 2.08 23.25
4 4 Proposed 24 × 256 0.33 21.84 168 0.67 0.04 0.94 18.20
WU 24 × 256 0.33 23.00 255 0.68 0.14 1.76 24.90
__________________ _______________________ __________________________ __________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________ __________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________ __________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________ __________________________ _____________________ ___________________ __________________ ___________________
7 6 Proposed 48 × 300 0.33 93.05 299 0.62 0.00 1.76 27.15
WU 48 × 300 0.33 41.99 299 0.75 0.00 2.52 37.94
8 6 Proposed 48 × 300 0.33 93.05 299 0.71 0.00 1.74 21.65
WU 48 × 300 0.33 41.99 299 0.77 0.01 2.58 44.09
Comparison with Wu’s Designs

We use the 12 × 64, 24 × 256 proposed designs in Example 1 directly. We obtain the 48 × 300 proposed design by selecting the intercept and the first 299 balanced columns of the 48 × 8192 design in Example 2. We obtain Wu’s 12 × 64, 24 × 256 and 48 × 300 designs by selecting the first 64 columns of Wu’s 12 × 67 design, the first 256 columns of Wu’s 24 × 277 design and the first 300 columns of Wu’s 48 × 1129 design, respectively. The result shown in Table 8 indicates that the proposed designs are generally better than Wu’s designs in terms of variable selection and parameter estimation with the Lasso. In almost all cases, the proposed designs have better performance while in the first two cases, Wu’s designs perform relatively better.

 Comparison with the Bayesian D-Optimal Supersaturated Designs

We use the 12 × 64 proposed design in Example 1 directly. We obtain the 24 × 253 proposed design by selecting the intercept, the 168 balanced columns and the first 84 unbalanced columns of the 24 × 256 design in Example 1. We first search a 12 × 64 and a 24 × 256 Bayesian D-optimal supersaturated design by the JMPő software [13] with 500 random starting designs. We use the 12 × 64 design directly. We remove the 22th, 134th and 209th columns of the 24 × 256 design to obtain the 24 × 253 design used in the simulations since the original searched 24 × 256 Bayesian D-optimal supersaturated design has μ = 0.83. The result shown in Table 9 indicates that the proposed designs generally outperform the Bayesian D-optimal supersaturated designs in terms of variable selection and parameter estimation with the Lasso.

Comparison with the Bayesian D-Optimal Supersaturated Designs.

Case #Active Design Size μ E ( s 2 ) #Balance AFDR AMR MSE EME
1 2 Proposed 12 × 64 0.33 9.90 24 0.65 0.00 0.95 8.41
BAYES 12 × 64 0.67 11.04 63 0.76 0.06 45.02 76.31
2 2 Proposed 12 × 64 0.33 9.90 24 0.66 0.05 1.20 12.46
BAYES 12 × 64 0.67 11.04 63 0.78 0.45 3.20 16.31
3 4 Proposed 12 × 64 0.33 9.90 24 0.30 0.00 1.83 18.04
BAYES 12 × 64 0.67 11.04 63 0.66 0.16 240.96 535.21
4 4 Proposed 12 × 64 0.33 9.90 24 0.65 0.06 1.48 12.77
BAYES 12 × 64 0.67 11.04 63 0.63 0.46 4.61 30.40
__________________ _______________________ __________________________ __________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________ __________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________ __________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________ __________________________ _____________________ ___________________ _____________________ _____________________
5 4 Proposed 24 × 253 0.33 21.84 168 0.60 0.00 1.61 21.59
BAYES 24 × 253 0.67 22.91 252 0.67 0.00 3.62 27.99
6 4 Proposed 24 × 253 0.33 21.84 168 0.70 0.05 1.85 29.70
BAYES 24 × 253 0.67 22.91 252 0.61 0.28 3.36 50.61
7 6 Proposed 24 × 253 0.33 21.84 168 0.54 0.00 4.24 38.37
BAYES 24 × 253 0.67 22.91 252 0.67 0.28 162.70 934.85
8 6 Proposed 24 × 253 0.33 21.84 168 0.72 0.01 3.49 25.34
BAYES 24 × 253 0.67 22.91 252 0.70 0.44 16.57 158.04

Comparison with the U E ( s 2 )-Optimal Designs.

Case #Active Design Size μ E ( s 2 ) #Balance MAFDR MAMR MMSE MEME
1 2 Proposed 12 × 64 0.33 9.90 24 0.67 0.00 0.87 7.99
JM 12 × 64 0.83 9.90 16 0.61 0.00 2.53 11.02
2 2 Proposed 12 × 64 0.33 9.90 24 0.67 0.03 1.04 10.38
JM 12 × 64 0.83 9.90 16 0.66 0.26 2.35 15.26
3 4 Proposed 12 × 64 0.33 9.90 24 0.30 0.00 2.03 21.28
JM 12 × 64 0.83 9.90 16 0.62 0.45 351.16 956.38
4 4 Proposed 12 × 64 0.33 9.90 24 0.59 0.09 1.44 13.18
JM 12 × 64 0.83 9.90 16 0.66 0.55 4.47 24.95
__________________ _______________________ __________________________ __________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________ __________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________ __________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________ __________________________ _________________________ _______________________ ______________________ ________________________
5 4 Proposed 24 × 256 0.33 21.84 168 0.58 0.00 1.68 26.33
JM 24 × 256 0.58 21.84 44 0.69 0.00 66.86 61.37
6 4 Proposed 24 × 256 0.33 21.84 168 0.72 0.02 1.42 20.00
JM 24 × 256 0.58 21.84 44 0.73 0.31 3.73 41.41
7 6 Proposed 24 × 256 0.33 21.84 168 0.71 0.30 260.42 398.54
JM 24 × 256 0.58 21.84 44 0.72 0.32 294.96 1336.26
8 6 Proposed 24 × 256 0.33 21.84 168 0.62 0.35 13.62 110.92
JM 24 × 256 0.58 21.84 44 0.72 0.43 18.54 147.44
Comparison with the <inline-formula id="j_nejsds34_ineq_269"><alternatives><mml:math> <mml:mi mathvariant="italic">U</mml:mi> <mml:mi mathvariant="italic">E</mml:mi> <mml:mo mathvariant="normal" fence="true" stretchy="false">(</mml:mo> <mml:msup> <mml:mrow> <mml:mi mathvariant="italic">s</mml:mi> </mml:mrow> <mml:mrow> <mml:mn>2</mml:mn> </mml:mrow> </mml:msup> <mml:mo mathvariant="normal" fence="true" stretchy="false">)</mml:mo></mml:math><tex-math><![CDATA[$UE({s^{2}})$]]></tex-math></alternatives></inline-formula>-Optimal Designs

According to Jones and Majumdar [15], for given n and p, there are a large number of U E ( s 2 )-optimal designs. In practical scenarios, people may randomly generate only one U E ( s 2 )-optimal design and use it to screen important factors. In our comparison, we want to mimic such a scenario, i.e., we want to know whether people will have a higher chance ( > 0.5) of getting a worse design (compared with our design) if they generate a U E ( s 2 )-optimal design in this way. Therefore, we randomly generate 500 U E ( s 2 )-optimal designs and for each of them, repeat the simulation procedure N = 30 times (rather than N = 300 in comparisons with other competing designs). We take the median values for the four criteria introduced in the beginning of the section over the 500 randomly generated designs and compare these median values with the corresponding criteria of the proposed design. If the AFDR of a proposed design is less than the median AFDR (MAFDR) of all the 500 U E ( s 2 )-optimal designs, it means the AFDR of the proposed design is less than that of the U E ( s 2 )-optimal design with more than half probability. The coherence, E ( s 2 ) and number of balanced columns in the table are also the corresponding medians. In addition, it is possible that a generated U E ( s 2 )-optimal design has μ = 1 and we will avoid using these designs in the simulations, because in practice if people get a bad design with μ = 1, they will probably regenerate another one.

We use the 12 × 64 and 24 × 256 proposed designs in Example 1 directly. We obtain the U E ( s 2 )-optimal designs by randomly sampling 12 rows of a 64 × 64 Hadamard matrix 500 times and 24 rows of a 256 × 256 Hadamard matrix 500 times, respectively. The result shown in Table 10 indicates that with more than half probability, the proposed designs outperform the U E ( s 2 )-optimal designs in terms of variable selection and parameter estimation with the Lasso. Furthermore, one may want to see the distributions of the criteria of 500 U E ( s 2 )-optimal designs. We take Case 8 for an example and show the distributions of the criteria in Figure 1.

Comparison between the U E ( s 2 )-optimal design and the proposed design in Case 8. The solid curves are the kernel density estimates of the probability density functions of U E ( s 2 )-optimal design’s criteria and the solid lines are the corresponding median values. The dashed lines are the values of the proposed design’s criteria.

 Summary of Comparisons

To help people better understand how the proposed design performs compared to competing designs with the same size and under the same active coefficient setting, we collect the comparisons between designs with exactly the same size and active coefficient setting all together in Table 11. Lin’s design is not of size 12 × 64, so it is not included in the table (it is actually impossible to construct a design of this size with Lin’s method).

Comparison Between Designs with the Same Size and Active Coefficients.

Case #Active Design Size μ E ( s 2 ) #Balance AFDR AMR MSE EME
1 2 Proposed 12 × 64 0.33 9.90 24 0.65 0.00 0.95 8.41
WU 12 × 64 0.33 11.06 63 0.54 0.00 0.89 9.03
BAYES 12 × 64 0.67 11.04 63 0.76 0.06 45.02 76.31
JM 12 × 64 0.83 9.90 16 0.61 0.00 2.53 11.02
2 2 Proposed 12 × 64 0.33 9.90 24 0.66 0.05 1.20 12.46
WU 12 × 64 0.33 11.06 63 0.56 0.04 0.98 11.22
BAYES 12 × 64 0.67 11.04 63 0.78 0.45 3.20 16.31
JM 12 × 64 0.83 9.90 16 0.66 0.26 2.35 15.26

According to Table 11, the proposed design significantly outperforms the Bayesian D-optimal design and the U E ( s 2 )-optimal design. In this scenario, Wu’s design performs slightly better than the proposed design, but there are also scenarios where the proposed design performs better than Wu’s design (see Table 8).

In summary, for Wu’s design, the proposed design is generally comparable to it, because of the same coherence they have. For Lin’s design, the Bayesian D-optimal design and the U E ( s 2 )-optimal design, the proposed design generally outperforms them, except that it sometimes tends to select more factors and gives a larger AFDR. However, the proposed design always gives a smaller AMR when compared to these three designs. In practice, a small AMR is more important than a small AFDR, because missing active factors is a more serious problem than falsely discovering inactive factors. This is because we can conduct follow-up experiments to screen out spurious factors, but we can never detect active factors that were already removed in the screening phase [21].

 Conclusion and Discussion

We have proposed a method for constructing supersaturated designs with small coherence. The constructed designs are allowed to be unbalanced to achieve more flexible sample sizes. The proposed method uses direct constructions and it entails no intensive computing even for large p. Since the proposed method can efficiently construct a supersaturated design with a large number of columns, it can be applied to high-dimensional variable selection problems in marketing, biology, engineering and other areas.

Here are possible directions for future work. First, the proposed method can expand a 6 m × p design D 0 with μ ≤ 1 / 3 to a larger design with μ ≤ 1 / 3. One may be interested in finding a special D 0 such that the expanded design has coherence strictly smaller than 1 / 3 or using a different construction method to reduce the coherence upper bound. Second, beyond the Lasso, the proposed design can be applicable to other variable selection and high-dimensional problems where controlling coherence of the design matrix is important. Moreover, many new powerful variable selection techniques have emerged. For example, Fan and Lv [10] and Fan, Feng and Song [11] proposed the sure independence screening and nonparametric independence screening for ultrahigh dimensional variable selection. Another example is the best subset regression. Shen, Pan, Zhu and Zhou [25] provided theoretical support for the use of best subset regression. Bertsimas, King and Mazumder [1] proposed a mixed integer optimization approach for the best subset selection problem, where they showed that coherence plays a role in parameter specification of their approach. In future work, our proposed designs will be applied to other penalized variable selection methods with results reported elsewhere.

 Proof for the Upper Bound of Coherence

For any matrix with levels − and +, if the signs of all the entries in a row are changed simultaneously, then the dot product of every two columns remains the same. Thus, the coherence of the matrix remains the same.

The coherence of D 1 constructed in Section 2.2 is no greater than 1 / 3.

Throughout, for two columns d i and d j of a two-level design D with n rows, we call | d i ⊤ d j | / n their absolute column correlation. Since D 1 is obtained by two copies of D 0 , the absolute column correlations of any two columns of D 1 are the same as those of D 0 . Thus, the coherence of D 1 equals that of D 0 , which is no greater than 1 / 3.  □

The coherence of D 2 constructed in Section 2.2 is no greater than 1 / 3.

Let d ( i ) denote the ith column of D 2 and d j ( i ) denote the jth entry of d ( i ) . Put

where C 1 and C 2 are 12 m × p matrices.

Pick two arbitrary columns of D 2 .

Suppose that both of them are from C 1 or C 2 . By Lemma 2, μ ( C 1 ) is ≤ 1 / 3. By Lemmas 1 and 2, μ ( C 2 ) is ≤ 1 / 3. Thus, the absolute column correlation of the two columns is ≤ 1 / 3.

Suppose that one of the two columns is the ith column d ( i ) of C 1 and the other column is the jth column d ( p + j ) of C 2 with 1 ≤ i , j ≤ p. Partition the two columns as d ( i ) d ( p + j ) = A B , where A and B have 8 m and 4 m rows, respectively. By Step 2 in Section 2.2,

Note that A consists of 2 m blocks with 4 rows and 2 columns − d 2 l − 1 ( i ) − d 2 l − 1 ( j ) d 2 l ( i ) d 2 l ( j ) − d 2 l − 1 ( i ) − d 2 l − 1 ( j ) d 2 l ( i ) d 2 l ( j ) , l = 1 , … , 2 m .

Since the absolute dot products of the two columns in these blocks are ≤ 2, the absolute dot product of the two columns in A is ≤ 4 m. Similarly, B consists of m blocks with 4 rows and 2 columns − d 8 m + 2 l − 1 ( i ) − d 8 m + 2 l − 1 ( j ) d 8 m + 2 l ( i ) − d 8 m + 2 l ( j ) − 1 − d 8 m + 2 l − 1 ( i ) − d 8 m + 2 l − 1 ( j ) d 8 m + 2 l ( i ) − d 8 m + 2 l ( j ) − 1 , l = 1 , … , m . Since the absolute dot products of the two columns in these blocks are 0, the absolute dot product of the two columns in B is 0. Thus, the absolute dot product of d ( i ) and d ( p + j ) is ≤ 4 m and their absolute column correlation is ≤ 4 m / 12 m = 1 / 3.

 Combining ( i ) and ( ii ) proves that μ ( D 2 ) is ≤ 1 / 3.  □

The coherence of D 3 constructed in Section 2.2 is no greater than 1 / 3.

Let d ( i ) denote the ith column of D 3 and d j ( i ) denote the jth entry of d ( i ) . Put

where C 1 , C 2 , C 3 and C 4 are 12 m × p matrices.

Pick two arbitrary columns of D 3 .

Suppose that both of them are from [ C 1 , C 2 ] or [ C 3 , C 4 ]. By Lemma 3, the coherence of [ C 1 , C 2 ] is ≤ 1 / 3. By Lemmas 1 and 3, the coherence of [ C 3 , C 4 ] is ≤ 1 / 3. Thus, the absolute column correlation of these two columns is ≤ 1 / 3.

Suppose that one of the two columns is the ith column d ( i ) of C 1 and the other column is the jth column d ( 2 p + j ) of C 3 with 1 ≤ i , j ≤ p. Partition these columns as d ( i ) d ( 2 p + j ) = A B , where A and B have 8 m and 4 m rows, respectively. By Step 3 in Section 2.2,

Note that A consists of 2 m blocks with 4 rows and 2 columns − d 2 l − 1 ( i ) − d 2 l − 1 ( j ) d 2 l ( i ) − d 2 l ( j ) − − d 2 l − 1 ( i ) − d 2 l − 1 ( j ) d 2 l ( i ) − d 2 l ( j ) − , l = 1 , … , 2 m . Since the absolute dot products of the two columns in these blocks are ≤ 2, the absolute dot product of the two columns in A is ≤ 4 m. Similarly, B consists of m blocks with 4 rows and 2 columns − d 8 m + 2 l − 1 ( i ) − d 8 m + 2 l − 1 ( j ) d 8 m + 2 l ( i ) − d 8 m + 2 l ( j ) − − d 8 m + 2 l − 1 ( i ) − . d 8 m + 2 l − 1 ( j ) d 8 m + 2 l ( i ) − d 8 m + 2 l ( j ) − , l = 1 , … , m . Since the absolute dot products of the two columns in these blocks are 0, the absolute dot product of the two columns in B is 0. Thus, the absolute dot product of d ( i ) and d ( 2 p + j ) is ≤ 4 m and their absolute column correlation is ≤ 4 m / 12 m = 1 / 3.

Suppose that one of the two columns is the ith column d ( p + i ) of C 2 and the other column is the jth column d ( 3 p + j ) of C 4 with 1 ≤ i , j ≤ p. This case can be proved similarly to ( ii ).

Suppose that one of the two columns is the ith column d ( i ) of C 1 and the other column is the jth column d ( 3 p + j ) of C 4 with 1 ≤ i , j ≤ p. Partition these columns as d ( i ) d ( 3 p + j ) = A B , where A and B have 8 m and 4 m rows, respectively. By Step 3 in Section 2.2,

Note that A consists of 2 m blocks with 4 rows and 2 columns − d 2 l − 1 ( i ) − d 2 l − 1 ( j ) d 2 l ( i ) − d 2 l ( j ) − − d 2 l − 1 ( i ) − d 2 l − 1 ( j ) d 2 l ( i ) − d 2 l ( j ) − , l = 1 , … , 2 m . Since the absolute dot products of the two columns in these blocks are 0, the absolute dot product of the two columns in A is 0. Similarly, B consists of m blocks with 4 rows and 2 columns − d 8 m + 2 l − 1 ( i ) − d 8 m + 2 l − 1 ( j ) d 8 m + 2 l ( i ) − d 8 m + 2 l ( j ) − − d 8 m + 2 l − 1 ( i ) − d 8 m + 2 l − 1 ( j ) d 8 m + 2 l ( i ) − d 8 m + 2 l ( j ) − , l = 1 , … , m . The absolute dot product of the two columns in the lth block is 2 | − d 8 m + 2 l − 1 ( i ) d 8 m + 2 l − 1 ( j ) + d 8 m + 2 l ( i ) d 8 m + 2 l ( j ) | , which is either 0 or 4. This implies that the absolute dot product of the two columns in B is ≤ 4 m. Thus, the absolute dot product of d ( i ) and d ( 3 p + j ) is ≤ 4 m and their absolute column correlation is ≤ 4 m / 12 m = 1 / 3.

Suppose that one of the two columns is the ith column d ( p + i ) of C 2 and the other column is the jth column d ( 2 p + j ) of C 3 with 1 ≤ i , j ≤ p. This case can be proved similarly to ( iv ).

 Combining ( i ), ( ii ), ( iii ), ( iv ) and ( v ) proves that μ ( D 3 ) is ≤ 1 / 3.  □

The coherence of D 2 constructed in Section 2.4 is no greater than t / n.

To prove this theorem, we use the same notation from the proof of Lemma 3, and pick two arbitrary columns of D 2 .

Suppose that both of the two columns are from C 1 or C 2 . By a similar argument to Lemma 2, the coherence of C 1 is ≤ t / n. Further, by Lemma 1, the coherence of C 2 is ≤ t / n. Thus, the absolute column correlation of these two columns is ≤ t / n.

Suppose that one of the two columns is from C 1 and the other column is from C 2 . Since the matrices A and B have 4 t and 2 n− 4 t rows respectively, both of which are multiples of 4, a similar argument to Lemma 3 (ii) still works. Hence, the absolute column correlation of these two columns is ≤ ( 4 t / 4 × 2 ) / 2 n = t / n.

 Combining (i) and ( ii) proves that the coherence of D 2 is ≤ t / n.  □

 Coherence and Model Selection Consistency of the Lasso

Zhao and Yu [34] studied the model selection consistency of the Lasso, i.e., the ability of the Lasso to exactly identify all the active factors from a large number of factors as the sample size n gets large. Assuming the model in (1.1), they defined the Strong Irrepresentable Condition, which requires the existence of a positive constant vector η such that | C 21 C 11 − 1 sign ( β ( 1 ) ) | ≤ 1 − η , where C 21 = X ( 2 ) ⊤ X ( 1 ) / n, C 11 = X ( 1 ) ⊤ X ( 1 ) / n is invertible, X ( 1 ) and X ( 2 ) are the submatrices consisting of the columns of X corresponding to the active and inactive factors respectively, the element-wise function sign ( · ) maps zero to zero, positive numbers to + 1 and negative numbers to − 1, β ( 1 ) is the vector of coefficients of active factors, 1 is a vector of + 1’s, and the inequality holds element-wise. They showed that for fixed p and β, under regularity conditions lim n → + ∞ 1 n X ⊤ X = C , C is a positive definite matrix , lim n → + ∞ 1 n max 1 ≤ i ≤ n x ( i ) ⊤ x ( i ) = 0 , where x ( i ) is the ith row of X, if the Strong Irrepresentable Condition holds, then there exists a λ n = f ( n ) such that lim n → + ∞ P ( sign ( β ˆ ( λ n ) ) = sign ( β ) ) = 1 , where β ˆ ( λ n ) is the Lasso estimate. Further, in Corollary 2 of their paper, they showed that if the worst case column correlation is no greater than c / ( 2 ‖ β ‖ 0 − 1 ) for a constant 0 ≤ c < 1, then the Strong Irrepresentable Condition holds. In sum, a small worst case column correlation implies the Strong Irrepresentable Condition, which in turn guarantees the model selection consistency of the Lasso.

 References Bertsimas, D., King, A. and Mazumder, R. Best subset selection via a modern optimization lens. The Annals of Statistics 44(2) 813–852 (2016). https://doi.org/10.1214/15-AOS1388. MR3476618 Booth, K. H. and Cox, D. R. Some systematic supersaturated designs. Technometrics 4(4) 489–495 (1962). https://doi.org/10.2307/1266285. MR0184369 Box, G. E. and Meyer, R. D. An analysis for unreplicated fractional factorials. Technometrics 28(1) 11–18 (1986). https://doi.org/10.2307/1269599. MR0824728 Bulutoglu, D. A. and Cheng, C. S. Hidden projection properties of some nonregular fractional factorial designs and their applications. The Annals of Statistics 31(3) 1012–1026 (2003). https://doi.org/10.1214/aos/1056562472. MR1994740 Chai, F. S., Das, A., Singh, R. and Stufken, J. Discriminating between superior U E ( s 2 )-optimal supersaturated designs. Statistics and Applications 18(2) (2020). Chen, J. and Lin, D. K. On the identifiability of a supersaturated design. Journal of Statistical Planning and Inference 72(1–2) 99–107 (1998). https://doi.org/10.1016/S0378-3758(98)00025-1. MR1655185 Davenport, M. A., Duarte, M. F., Eldar, Y. C. and Kutyniok, G. Introduction to Compressed Sensing. Citeseer (2012). https://doi.org/10.1017/CBO9780511794308.002. MR2963166 Deng, L. Y. and Lin, D. K. Criteria for supersaturated designs. In Proceedings of the Section on Physical and Engineering Sciences 124–128. American Statistical Association (1994). DuMouchel, W. and Jones, B. A simple Bayesian modification of D-optimal designs to reduce dependence on an assumed model. Technometrics 36(1) 37–47 (1994). Fan, J. and Lv, J. Sure independence screening for ultrahigh dimensional feature space. Journal of the Royal Statistical Society, Series B 70(5) 849–911 (2008). https://doi.org/10.1111/j.1467-9868.2008.00674.x. MR2530322 Fan, J., Feng, Y. and Song, R. Nonparametric independence screening in sparse ultra-high-dimensional additive models. Journal of the American Statistical Association 106(494) 544–557 (2011). https://doi.org/10.1198/jasa.2011.tm09779. MR2847969 Friedman, J., Hastie, T. and Tibshirani, R. Regularization paths for generalized linear models via coordinate descent. Journal of Statistical Software 33(1) 1–22 (2010). JMP (version 11). SAS Institute Inc., Cary, NC (1989–2007). Jones, B., Lin, D. K. and Nachtsheim, C. J. Bayesian D-optimal supersaturated designs. Journal of Statistical Planning and Inference 138(1) 86–92 (2008). https://doi.org/10.1016/j.jspi.2007.05.021. MR2369616 Jones, B. and Majumdar, D. Optimal supersaturated designs. Journal of the American Statistical Association 109(508) 1592–1600 (2014). https://doi.org/10.1080/01621459.2014.938810. MR3293612 Jones, B., Lekivetz, R., Majumdar, D., Nachtsheim, C. J. and Stallrich, J. W. Construction, properties, and analysis of group-orthogonal supersaturated designs. Technometrics 62(3) 403–414 (2020). https://doi.org/10.1080/00401706.2019.1654926. MR4125505 Li, R. and Lin, D. K. Data analysis in supersaturated designs. Statistics & Probability Letters 59(2) 135–144 (2002). https://doi.org/10.1016/S0167-7152(02)00140-2. MR1927525 Lin, D. K. A new class of supersaturated designs. Technometrics 35(1) 28–31 (1993). Lin, D. K. Generating systematic supersaturated designs. Technometrics 37(2) 213–225 (1995). https://doi.org/10.2307/1269765. MR1365722 Liu, Y., Ruan, S. and Dean, A. M. Construction and analysis of E ( s 2 ) efficient supersaturated designs. Journal of Statistical Planning and Inference 137(5) 1516–1529 (2007). https://doi.org/10.1016/j.jspi.2006.09.001. MR2303773 Marley, C. J. and Woods, D. C. A comparison of design and model selection methods for supersaturated experiments. Computational Statistics & Data Analysis 54(12) 3158–3167 (2010). https://doi.org/10.1016/j.csda.2010.02.017. MR2727742 Osborne, M. R., Presnell, B. and Turlach, B. A. Knot selection for regression splines via the lasso. Computing Science and Statistics. 44–49 (1998). Plackett, R. L. and Burman, J. P. The design of optimum multifactorial experiments. Biometrika 33(4) 305–325 (1946). https://doi.org/10.1093/biomet/33.4.305. MR0016624 R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria (2017). https://www.R-project.org/. Shen, X., Pan, W., Zhu, Y. and Zhou, H. On constrained and regularized high-dimensional regression. Annals of the Institute of Statistical Mathematics 65(5) 807–832 (2013). https://doi.org/10.1007/s10463-012-0396-3. MR3105798 Singh, R. and Stufken, J. Selection of two-level supersaturated designs for main effects models. Technometrics 65 1–9 (2022). https://doi.org/10.1080/00401706.2022.2102080. MR4543063 Tang, B. and Wu, C. F. J. A method for constructing supersaturated designs and its E ( s 2 ) optimality. Canadian Journal of Statistics 25(2) 191–201 (1997). https://doi.org/10.2307/3315731. MR1463319 Tibshirani, R. Regression shrinkage and selection via the Lasso. Journal of the Royal Statistical Society, Series B 58 267–288 (1996). MR1379242 Weese, M. L., Smucker, B. J. and Edwards, D. J. Searching for powerful supersaturated designs. Journal of Quality Technology 47(1) 66–84 (2015). Weese, M. L., Stallrich, J. W., Smucker, B. J. and Edwards, D. J. Strategies for supersaturated screening: Group orthogonal and constrained var (s) designs. Technometrics 63(4) 443–455 (2021). https://doi.org/10.1080/00401706.2020.1850529. MR4331445 Welch, L. R. Lower bounds on the maximum cross correlation of signals (corresp.). IEEE Transactions on Information Theory 20(3) 397–399 (1974). Wu, C. F. J. Construction of supersaturated designs through partially aliased interactions. Biometrika 80(3) 661–669 (1993). https://doi.org/10.1093/biomet/80.3.661. MR1248029 Wu, C. F. J. and Hamada, M. S. Experiments: planning, analysis, and optimization 2nd ed. John Wiley & Sons, Hoboken, NJ (2009). MR2583259 Zhao, P. and Yu, B. On model selection consistency of Lasso. Journal of Machine Learning Research 7 2541–2563 (2006). MR2274449