Orthogonal Arrays: A Review

C. Devon Lin1 and John Stufken2

1
Department of Mathematics and Statistics, Queen’s University
2
Department of Statistics, College of Engineering and Computing, George
arXiv:2505.15032v2 [[Link]] 5 Jun 2025

Mason University

Abstract
Orthogonal arrays are arguably one of the most fascinating and important statisti-
cal tools for efficient data collection. They have a simple, natural definition, desirable
properties when used as fractional factorials, and a rich and beautiful mathematical the-
ory. Their connections with combinatorics, finite fields, geometry, and error-correcting
codes are profound. Orthogonal arrays have been widely used in agriculture, engi-
neering, manufacturing, and high-technology industries for quality and productivity
improvement experiments. In recent years, they have drawn rapidly growing interest
from various fields such as computer experiments, integration, visualization, optimiza-
tion, big data, machine learning/artificial intelligence through successful applications
in those fields. We review the fundamental concepts and statistical properties and
report recent developments. Discussions of recent applications and connections with
various fields are presented.
Key Words: Big data, Computer experiment, Error-correcting code, Factorial
experiment, Hadamard matrix, Subsampling

1 Definition and History

Orthogonal arrays have applications in many areas and have proven to be a fascinating
and rich subject for research. Statisticians, mathematicians and other researchers have
studied orthogonal arrays since their introduction by C. R. Rao in a series of seminar papers
(Rao, 1946, 1947, 1949). In the field of design of experiments they are used to determine the
settings of factors for conducting experiments, so that effects of these factors on a response
variable can be explored simultaneously. The possible settings for these factors are called
levels. Let S be a set of s levels, conventionally denoted by 0, 1, . . . , s − 1. Formally, we can
define an orthogonal array as follows.

Definition 1.1. An N × k array A with entries from S is said to be an orthogonal array

with s levels, strength t, and index λ if every N × t subarray of A contains each t-tuple based
on S exactly λ times as a row.

1
 We use OA(N, sk , t) to denote such an orthogonal array. Note that λ is not included
in the notation because it can be derived as λ = N/st . An OA(N, sk , t) can be used to
determine N input settings or level combinations for k factors each with s levels. Each row
of an OA(N, sk , t) is called a run while each column corresponds to the setting of a factor or
input variable.

Example 1.1. Table 1 lists an OA(8, 24 , 3), an orthogonal array with two levels, strength
three, and index unity. It has eight runs and it is for four factors. It has the property that
every subarray consisting of its three distinct columns contains all the eight level combinations
exactly once.

Table 1: An OA(8, 24 , 3)

0 0 0 0
0 0 1 1
0 1 0 1
0 1 1 0
1 0 0 1
1 0 1 0
1 1 0 0
1 1 1 1

When all factors have the same number of levels, the arrays are often referred to as
“fixed-level” or “pure-level” or symmetrical orthogonal arrays. In some experiments, different
factors are allowed to have different numbers of levels, and this leads to a class of arrays
named mixed (or asymmetrical) orthogonal arrays. They began to receive more attention
in the early 1960s, especially through the consideration of asymmetrical orthogonal main
effects plans in Addelman and Kempthorne (1962). Rao (1973) discussed, among others, the
notion of asymmetrical orthogonal arrays. Significant advances have been made during the
past decades on their constructions (Wu, 1991; Wang and Wu, 1991; Hedayat et al., 1992;
Sitter, 1993; DeCock and Stufken, 2000; Suen et al., 2001; Pang et al., 2021), resulting in
wider applications. Formally, mixed (or asymmetrical) arrays can be defined as follows.

Definition 1.2. A mixed orthogonal array OA(N, sk11 sk22 · · · skvv , t) is an N × k array where
k = k1 + k2 + · · · + kv is the total number of factors, in which the first k1 factors have s1
levels, the next k2 factors have s2 levels, and so on, with the property that in any N × t
subarray every possible t-tuple occurs an equal number of times as a row.

Example 1.2. Table 2 lists an OA(12, 24 31 , 2), a mixed orthogonal array of strength two,
with 12 runs, with the first four factors at two levels and the fifth at three levels.

The term orthogonal arrays was first introduced by Bush (1950), although the concept
was previously described in Rao’s 1943 master thesis and three papers (Rao, 1946, 1947,
1949), first for the special case of what was called a hypercube of strength d (Rao, 1946)

2
 Table 2: A mixed orthogonal array OA(12, 24 31 , 2) (transposed)

0 0 1 1 0 0 1 1 0 0 1 1
0 1 0 1 0 1 0 1 0 1 0 1
0 0 1 1 1 1 0 0 1 0 0 1
0 1 0 1 1 0 0 1 1 0 1 0
0 0 0 0 1 1 1 1 2 2 2 2

and then for general orthogonal arrays. Orthogonal arrays have since found applications in
various fields, including as fractional factorials for agricultural, medical, industrial, and other
experiments. They are particularly useful in scenarios where the number of input variables
is relatively small but too large to permit exhaustive testing of all possible combinations. In
software testing, orthogonal arrays are effective in identifying errors related to faulty logic.
In quality control, they are closely associated with the Taguchi methods, developed in the
early 1950s and widely adopted by industries in Japan and later the United States. These
methods have been instrumental in improving the quality of manufactured goods, reducing
costs, and accelerating the time to market. Orthogonal arrays have also been applied in
engineering, biotechnology, marketing, advertising, and many other modern technological
fields.

2 Orthogonal Arrays as Fractional Factorial Experi-

ments
Orthogonal arrays were introduced for their use in fractional factorial experiments, and
this remains their most important application in statistics. An N × k orthogonal array
can be used to perform a fractional factorial experiment with N runs and k factors, where
the number of levels for a factor is equal to the number of symbols for the corresponding
column. Typically, the number of runs N would be much smaller than the number of all
possible level combinations for the k factors, and if all runs of the orthogonal array are
distinct, in which case it is called a simple orthogonal array, then its runs form a subset
of those in a full factorial which consists of all possible level combinations for the factors.
However, the definition of an orthogonal array does not require distinct runs, and we consider
it to represent a fractional factorial as long as N is smaller than the number of all possible
level combinations.
A well-known class of fractional factorials in which all factors have the same number of
levels s consists of the regular fractional factorials (cf. Box and Hunter, 1961, for s = 2).
Regular fractional factorials are a subclass of orthogonal arrays, and are closely related
to linear orthogonal arrays (cf. Hedayat et al., 2012, Chapter 4). A fractional factorial is
said to be regular if N is a power of s, say N = sp , there are p factors so that every
possible level combination appears once for these p factors, and the columns for the other
k − p factors can be computed as explained below from the columns of the initial p factors.
The relationships that exist between different columns in a regular fractional factorial are

3
called alias relationships. For example, consider the OA(8, 24 , 3) in Table 1. Any two or
three columns of this orthogonal array are orthogonal in the sense that every possible level
combination appears twice (for two columns) or once (for three columns). Selecting any
three columns, the levels for the fourth column can then be obtained by observing that
d1 + d2 + d3 + d4 = 0 (modulo 2), where the four columns are denoted as d1 , d2 , d3 , d4 .
Because of this relationship, the four-factor interaction is part of the defining relation of this
fractional factorial. It also implies that, with all computations modulo 2, d1 = d2 + d3 + d4 ,
d2 = d1 + d3 + d4 , d3 = d1 + d2 + d4 and d4 = d1 + d2 + d3 . This is to be interpreted as
the main effect of the first factor being aliased with (i.e., indistinguishable from) the three-
factor interaction of factors 2, 3, and 4; the main effect of the second factor being aliased
with the three-factor interaction of factors 1, 3, and 4; and so on. Every regular fractional
factorial has a defining relation. The relationships between the factors are captured by the
alias structure. For example, if we wish to construct a regular fractional factorial of 16 runs
with 7 two-level factors, we can start with a full factorial of 16 runs and 4 factors. Denoting
these factors by A, B, C, D, we can define the additional 3 factors using these four columns.
Rather than using levels 0 and 1, for 2-level fractional factorials, it is common practice to
use levels −1 and 1, which we will follow in this example. One possible choice for defining
three remaining factors is E = AB, F = BCD, G = AD. As such, we obtain the defining
relation I = ABE = BCDF = ADG = ACDEF = BDEG = ABCF G = CEF G, where
I represents a column of all ones. Here, ABE, BCDF , ADG, ACDEF , BDEG, ABCF G,
and CEF G are called words, and correspond to interactions in the defining relation that are
equal to I. Multiplying each term in the defining relation by A, and using that A2 = 1, we
obtain that A = BE = ABCDF = DG = CDEF = ABDEG = BCF G = ACEF G. This
is part of the alias structure, implying that if we try to estimate the main effect of factor A,
we are really estimating A+BE +ABCDF +DG+CDEF +ABDEG+BCF G+ACEF G.
We may interpret this as the main effect of A if the interactions in this sum can be assumed
to be negligible. The complete alias structure of this fractional factorial is,
A = BE = ABCDF = DG = CDEF = ABDEG = BCF G = ACEF G B = AE = CDF = ABDG = ABCDEF = DEG = ACF G = BCEF G
C = ABCE = BDF = ACDG = ADEF = BCDEG = ABF G = EF G D = ABDE = BCF = AG = ACEF = BEG = ABCDF G = CDEF G
E = AB = BCDEF = ADEG = ACDF = BDG = ABCEF G = CF G F = ABEF = BCD = ADF G = ACDE = BDEF G = ABCG = CEG
G = ABEG = BCDF G = AD = ACDEF G = BDE = ABCF = CEF AC = BCE = ABDF = CDG = DEF = ABCDEG = BF G = AEF G
AF = BEF = ABCD = DF G = CDE = ABDEF G = BCG = ACEG BC = ACE = DF = ABCDG = ABDEF = CDEG = AF G = BEF G
BD = ADE = CF = ABG = ABCEF = EG = ACDF G = BCF EF G BF = AEF = CD = ABDF G = ABCDE = DEF G = ACG = BCEG
BG = AEG = CDF G = AGD = ABCDEF G = DE = ACF = BCEF CE = ABC = BDEF = ACDEG = ADF = BCDG = ABEF G = F G
CG = ABCEG = BDF G = ACD = ADEF G = BCDE = ABF = EF

Together with the defining relation, these 15 equations establish the relationships between
the 27 = 128 effects including the grand mean effect, main effects, two-factor interaction
effects, and higher-order interaction effects. We refer to Chapters 4 and 5 of Wu and Hamada
(2011) and Chapter 9 of Cheng (2016) for in-depth discussions of regular fractional factorials.
Although the relationships between the columns of a regular fractional factorial are simple,
this simplicity comes at the cost of the run size having to be a power of the number of levels
s. For example, for a two-level regular fractional factorial, the run size must be a power of
2, i.e. 4, 8, 16, 32, 64, and so on, thereby creating an increasingly large gap between any
consecutive possible run sizes.
The general class of fixed-level orthogonal arrays does not have such a constraint on the
run size. For a two-level orthogonal array of strength two, for example, the run size can be

4
any multiple of 4, that is, 4, 8, 12, 16, 20, and so on. Such arrays are known as non-regular
fractional factorials if they are not regular. Their columns remain orthogonal (in the sense
that all level combinations appear equally often as a row for any set of up to t columns),
but factors are now partially aliased (in the sense that there are no alias relationships as for
regular factorials). For instance, Table 3 displays a two-level non-regular fractional factorial
of 12 runs, 11 factors, and strength two. Any two columns of are orthogonal. But, for
example, the levels for the factor corresponding to the third column cannot be computed
from those of the factors corresponding to the first two columns (or any other two columns),
because 00 for the first two factors can result in 0 or 1 for the third factor. Adding the
levels for the first three factors modulo 2 results eight times in a sum of 0 and four times
in a sum of 1, making the absolute correlation between the third factor and the interaction
of the first two 13 . Since this absolute correlation is not 0 or 1, we call this partial aliasing.
Similar observations hold for other sets of three columns.
Non-regular fractional factorials have a complex alias structure. Precisely because of
this, non-regular fractional factorials were traditionally not favored in choosing fractions.
However, with the advocacy of important work such as Hamada and Wu (1992); Lin and
Draper (1992); Cheng (1995); Deng and Tang (1999a); Hedayat et al. (2012); Dean et al.
(2015); Cheng and Tang (2025), the merits of non-regular fractional factorials are generally
well recognized and they are widely used in the design and analysis of experiments, and
beyond.

Table 3: An OA(12, 211 , 2)

0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 1 1 1 1 1
0 0 1 1 1 0 0 0 1 1 1
0 1 0 1 1 0 1 1 0 0 1
0 1 1 0 1 1 0 1 0 1 0
0 1 1 1 0 1 1 0 1 0 0
1 0 0 1 1 1 1 0 0 1 0
1 0 1 0 1 0 1 1 1 0 0
1 0 1 1 0 1 0 1 0 0 1
1 1 0 0 1 1 0 0 1 0 1
1 1 0 1 0 0 0 1 1 1 0
1 1 1 0 0 0 1 0 0 1 1

Definition 1.1 implies a key projection property of orthogonal arrays. Specifically, when
an orthogonal array of strength t is projected onto any subset of t or fewer factors, the
resulting subarray forms one or more replicates of the full factorial for those factors. This
projection property has an important statistical implication: if there are at most t active
(i.e., important) factors, then an orthogonal array of strength t allows for the estimation of
all factorial effects of these active factors. Notably, this holds regardless of which factors are
active.

5
 Deng and Tang (1999a) introduced the concept of J-characteristics for studying two-
level regular and non-regular fractional factorials. Let D = [d1 , d2 , . . . , dk ] = (dij ) denote
an N × k two-level factorial with levels 1 and −1. For 1 ≤ m ≤ k and any m-subset
S = {dj1 , dj2 , . . . , djm } of the columns in D, Deng and Tang (1999a) defined
N
X
Jm (S) = Jm (dj1 , dj2 , . . . , djm ) = | dij1 · · · dijm |. (1)
i=1

It can be shown that two-level orthogonal arrays of strength t have Jm (S) = 0 for all m-
subsets S and 1 ≤ m ≤ t. In addition, D is a regular fractional factorial if and only if
the J-characteristics of any subset is either 0 or N . Thus, if the J-characteristic of some
subset of a two-level fraction is strictly between 0 and N , then D is a non-regular fractional
factorial.
An immediate question when using an orthogonal array as a fractional factorial is how to
select the array (we will discuss this issue further in Section 7). In the pioneering work, Box
and Hunter (1961) introduced the concept of resolution to evaluate and compare two-level
fractional factorials. In the context of regular two-level fractional factorials, they defined a
fraction to be of resolution R if no c-factor effect is aliased with any other effect containing
less than R − c factors. For example, a fraction of resolution III does not alias main effects
with one another but allows main effects to be aliased with interactions of two or more
factors, and a fraction of resolution IV does not alias main effects with each other or with
two-factor interactions but does allow two-factor interactions to be aliased with each other.
In two-level regular fractional factorials, its resolution is the length of the shortest word in the
defining relation and is one more than the strength of the corresponding orthogonal array.
Deng and Tang (1999a) introduced generalized resolution to assess non-regular fractional
factorials. Let r be the smallest integer such that max|S|=r Jr (S) > 0, where Jr (S) is as
defined in (1) and the maximization is over all r-subsets S of r distinct columns of D. Deng
and Tang (1999a) defined the generalized resolution of design D to be

R(D) = r + [1 − max Jr (S)/N ]. (2)

|S|=r

As the values of Jr (S) are between 0 and N and the max in (2) is positive, we have r ≤
R(D) < r + 1. A regular or nonregular fractional factorial of resolution R is an orthogonal
array of strength ⌊R⌋ -1 where ⌊R⌋ is the largest integer that does not exceed R. In view of
this relationship, orthogonal arrays offer many attractive statistical properties when used as
fractional factorials. For example, an orthogonal array of strength t is universally optimal
for any model consisting of factorial effects involving at most ⌊t/2⌋ factors (Cheng, 1980).

3 Selected Applications of Orthogonal Arrays

Orthogonal arrays play a role in many applications. In Section 2, we have already ex-
plained the connection of orthogonal arrays to fractional factorials. Fractional factorials were
introduced by Finney (1945) and have found application in many fields, as already pointed
out in Section 1. Experiments using fractional factorials are especially prevalent in industry

6
for product development and improvement and quality control. In this section, we highlight
a few other applications with selected references for further reading.

3.1 Application in Numerical Integration

Orthogonal arrays are an important tool for numerical integration. This is due to the
structure imposed by their definition. Considering the rows of an s-level orthogonal array
of strength two as points in k-dimensional space, the array has the properties that (a) when
projecting the points on any coordinate axis, the projections are uniformly distributed over
s equally spaced points; and (b) when projecting the points on any coordinate plane, the
projections are uniformly distributed over s × s grids. In numerical integration, property
(a) helps filter out the main effects, while property (b) is useful for filtering out two-factor
interactions as well as main effects. This will become clear as we define main effects and
interaction effects below.
Let X = (X 1 , X 2 , . . . , X k ) be uniformly distributed in [0, 1]k and Y = f (X) ∈ R where
f is a known, but expensive to compute, function. In numerical integration, R a fundamental
problem is how to select N inputs X1 , X2 , . . . , XN to approximate f (X)dX accurately with
a minimal sample size N . Often, f (X)dX is estimated by Ȳ = N
R P
i=1 f (Xi )/N where f (Xi )
is the value of the function evaluated at the input Xi . The main effect of the input variable
X j is defined as fj (X j ) = E[f (X)|X j ] − E(Y ) and the interaction effect between variable
X i and variable X j is fij (X i , X j ) = E[f (X)|X i , X j ] − E(Y ) − fi (X i ) − fj (X j ), for 1 ≤
i ̸= j ≤ k. The simplest way to select N inputs is random sampling and the corresponding
variance of Ȳ is N −1 var[f (X)]. McKay et al. (1979) introduced Latin hypercube sampling,
which is based on a Latin hypercube that is an array that, by definition, satisfies that each
column is a permutation of 1, 2, . . . , N , where N is the run size of the Latin hypercube.
McKay et al. (1979) showed that Latin hypercube sampling achieves a smaller variance
for Ȳ than random sampling or stratified sampling when f satisfies certain monotonicity
conditions. Stein (1987) furtherPderived that the variance of Ȳ under Latin hypercube
sampling is N −1 var[f (X)] − N −1 ki=1 var[fi (X i )] + o(N −1 ) which is asymptotically smaller
than the variance of Ȳ under random sampling. Rows of different Latin hypercubes for the
same value of N can have very different projections on two-dimensional coordinate planes.
Owen (1992) and Tang (1993) independently took the idea in McKay et al. (1979) further
and introduced sampling methods based on orthogonal arrays to achieve further variance
reduction for numerical integration. Their methods are based on randomized orthogonal
arrays and orthogonal array-based Latin hypercubes, respectively. Tang (1993) proved that
the variance of Ȳ under Pk strength two orthogonal
Pk array-based Latin hypercube sampling is
−1 −1 i −1 −1
N var[f (X)] − N i=1 var[fi (X )] − N i<j var[fij ] + o(N ). Here we briefly explain
the procedure for obtaining an orthogonal array-based Latin hypercube from an orthogonal
array D = OA(N, sk , t). For each column of D, replace the N/s positions with entry u,
u = 0, 1, ..., s − 1, by a permutation of uN/s + 1, uN/s + 2, . . . , uN/s + N/s = (u + 1)N/s.
The resulting array has the property that each column is a permutation of 1, 2, . . . , N and
thus a Latin hypercube (McKay et al., 1979). Example 3.1 illustrates this procedure using
an OA(9, 34 , 2). Properties of sampling by using such an orthogonal array-based Latin
hypercube are discussed in Tang (1993), who refers to this procedure as U sampling. One
key conclusion is that if the underlying function f is additive, then U sampling gives a smaller

7
variance of the integral approximation than the sampling in Owen (1992).

Example 3.1. Consider D, an OA(9, 34 , 2), shown below. For each column of D, we replace
level 0 by a random permutation of 1, 2, 3, level 1 by a random permutation of 4, 5, 6, and
level 2 by a random permutation of 7, 8, 9. The orthogonal array-based Latin hypercube L
shown below is one of the arrays that we might obtain.
   
0 0 0 0 1 2 3 3
 0 1 1 2   2 6 4 9 
   
 0 2 2 1   3 8 7 6 
   
 1 0 1 1   4 3 6 4 
   
D=  1 1 2 0 , L =  5 5 9 1 .
  
 1 2 0 2   6 7 1 7 
   
 2 0 2 2   7 1 8 8 
   
 2 1 0 1   8 4 2 5 
2 2 1 0 9 9 5 2

3.2 Application in Computer Experiments

With the exponential growth of computing power, researchers are increasingly using com-
puter experiments to simulate real-world phenomena and complex systems through mathe-
matical models. These models are solved using numerical methods such as computational
fluid dynamics and finite element analysis to gain deeper insights into the systems being
studied (Santner et al., 2003; Gramacy, 2020). The underlying mechanisms of these com-
puter experiments are represented and executed through computer codes. To choose inputs
to run computer codes, a widely used approach is to use space-filling designs which aim to
spread out the design points evenly over the entire design space. Unlike traditional physical
experiments, computer experiments often involve a larger number of input variables and
require more runs. However, only a subset of these input variables is typically considered of
primary importance. To identify the key input variables, space-filling designs with desirable
low-dimensional projection properties are commonly used. These designs can be generated
using orthogonal array-based Latin hypercubes, which were described in Subsection 3.1. In
addition, Chen and Tang (2022) justified orthogonal array-based designs under a broad class
of space-filling criteria including commonly used distance-, orthogonality- and discrepancy-
based measures.
Orthogonal arrays are not only directly used in orthogonal array-based Latin hypercubes
but also serve as a cornerstone for various constructions of space-filling designs. For instance,
Steinberg and Lin (2006), Lin et al. (2009), Pang et al. (2009), Lin and Kang (2016), Sun and
Tang (2017a) and Sun and Tang (2017b) constructed another class of space-filling designs
known as orthogonal Latin hypercubes, which are Latin hypercubes with the property that
every two distinct columns have zero correlation. More specifically, for example, let us
consider how orthogonal arrays are used in Lin et al. (2009). For any positive integer u, let
gu be the u × 1 vector with ith element i − (u + 1)/2, 1 ≤ i ≤ u, and Γu be the set of the
u! vectors generated by permuting the elements of gu . Let B = (bij ) be an n × p matrix
with columns from Γn . Suppose that an orthogonal array OA(n2 , n2f , 2), say A, with n2

8
rows, 2f columns, n symbols, and strength two is available. Denote the symbols in A by
0, 1, 2, . . . , n − 1. Lin et al. (2009) proposed the following construction steps:

I. For 1 ≤ j ≤ p, obtain an n2 × (2f ) matrix Aj by replacing the symbols 0, 1, 2, . . . , n − 1

in A by b1j , b2j , . . . , bnj respectively, and then partition Aj as Aj = [Aj1 , . . . , Ajf ], where
each of Aj1 , . . . , Ajf has two columns.

II. For 1 ≤ j ≤ p, obtain the n2 × (2f ) matrix Mj = [Aj1 V, . . . , Ajf V ] , where

1 −n
h i
V = n 1 .

III. Finally, obtain the matrix M = [M1 , . . . , Mp ], of order N × q, where N = n2 and

q = 2pf .

Step II of this construction applies a rotation to pairs of columns of the orthogonal array
Aj . Lin et al. (2009) showed that the resulting matrix M is a Latin hypercube design in
which the correlation between any two distinct columns depends only on correlations of
columns of the matrix B. Thus, this construction can be used to obtain orthogonal Latin
hypercubes of larger run sizes from those with smaller run sizes.
In Section 8, we will discuss selected recent developments for orthogonal arrays. Some of
these are motivated by applications in computer experiments. For example, sliced orthogonal
arrays, introduced by Qian and Wu (2009), are a special class of orthogonal arrays with the
property that runs can be partitioned into smaller orthogonal arrays, and these arrays are
needed in computer experiments with both quantitative and qualitative inputs (Qian et al.,
2008).

3.3 Application in Subsampling of Big Data

With the advancement of technology, data generation continues to grow exponentially,
potentially resulting in huge datasets. As a result, many fields, including statistical science,
face unique challenges and unprecedented opportunities. One such challenge and opportunity
is the development of subsampling methods for efficiently selecting subdata (i.e., a subset of a
large dataset) with minimal loss of information. For example, Wang et al. (2019) proposed a
novel approach, termed information-based optimal subdata selection (IBOSS), in the context
of big data linear regression problems, and proposed a computationally efficient algorithm
for approximating the optimal subdata via the IBOSS method. Wang et al. (2021) also
considered a subsampling approach, but used orthogonal arrays to select optimal subdata.
Consider the linear regression model with p predictors or features,

y = β0 + β1 x1 + β2 x2 + . . . + βp xp + ϵ,

where y is the response and the random errors ϵ are independent and identically distributed
with mean 0 and variance σ 2 . The least squares estimate of β = (β0 , β1 , . . . , βp )T is β̂ =
(X̃ T X̃)−1 X̃ T y where y is the response vector, X̃ = (1, X) and X is the N × p matrix of
feature values for the N observations. Now take a subsample of size n from the full dataset

9
(X, y) and let (Xs , ys ) denote this subsample. The least squares estimate based on the
subsample is
β̂ s = (X̃sT X̃s )−1 X̃sT ys ,
where X̃s = (1, Xs ). The covariance matrix of β̂ s is σ 2 (X̃sT X̃s )−1 . Information-based sub-
sampling approaches aim to find subdata that, in some way, minimize the variance of β̂ s .
Using an optimality function, say ϕ, the optimal subdata minimizes ϕ((X̃sT X̃s )−1 ), i.e.,

Xs∗ = arg min ϕ((X̃sT X̃s )−1 ).

Xs ⊂X

Finding an exact solution for N >> n is too expensive, so that algorithms seek a highly
efficient solution for Xs∗ .
Common choices for ϕ are the determinant and trace, which correspond to the criteria
of D- and A-optimality, respectively, in optimal design of experiments. The orthogonal
subsampling proposed by Wang et al. (2021) selects the subsample Xs such that if all the
covariates are scaled to [-1,1], Xs mimics a two-level orthogonal array of strength two as
closely as possible. The method is inspired by the optimality of orthogonal arrays for linear
regression models. For example, Cheng (1980) showed that an orthogonal array of strength
two with s levels is universally optimal for a main-effects model. That is, such an array is
optimal under a wide range of criteria that include D- and A-optimality, among all s-level
factorial designs.
A subsample Xs for which the rows form exactly a two-level orthogonal array will gen-
erally not exist in the full dataset from which the sample is taken. Wang et al. (2021)
introduced a discrepancy function that is to be minimized in order to sequentially select
subdata that aims for the simultaneous attainment of two features. These two features are:
(i) select points with extreme values of the features: selected points are located near the
corners of the feature space and have a large distance from the center, and (ii) aim for or-
thogonality of columns corresponding to any two features. Wang et al. (2021) derived a lower
bound for their discrepancy function and proposed an efficient algorithm to sequentially se-
lect points for inclusion in the subdata that minimize the discrepancy function. Interested
readers are referred to their article for the details of this method and algorithm.
Zhang et al. (2024) and Zhu et al. (2024) extended the idea of the orthogonal subsampling
to independence-encouraging subsampling for nonparametric additive models and group-
orthogonal subsampling for hierarchical data based on linear mixed models, respectively. In
both cases, orthogonal arrays serve as an essential tool for the proposed subsamplings.

4 Orthogonal Arrays and Error-Correcting Codes

This section focuses on connections between orthogonal arrays and error-correcting codes.
These two concepts have deep connections as first observed in Bose (1961) for linear codes
and linear orthogonal arrays. The latter are orthogonal arrays with a defining relation (cf.
Hedayat et al., 2012, Section 11.5). An even deeper connection is based on the work by
Delsarte (1973). For further discussion, we begin by reviewing fundamental concepts of an
error-correcting code. Excellent references on this topic are MacWilliams and Sloane (1977),

10
Stinson (2006) and Hedayat et al. (2012). Error-correcting codes are used to detect and
correct errors that occur during data transmission over noisy communication channels. With
a set of symbols S of size s, called the alphabet, an error-correcting code is any collection C
of vectors from S k , the set of all sk vectors of length k based on the alphabet S. The vectors
in C are called the codewords. For codes it is not common that codewords are repeated (i.e.,
all vectors in C are typically distinct), but since this is not a requirement for orthogonal
arrays, we allow repetition of codewords for codes.
An important concept for a code C is its minimal distance, which is defined as

d = min dist(u, v),

u,v∈C
u̸=v

where dist(u, v) is the number of positions where vectors u and v differ, referred to as the
Hamming distance between u and v. A code with minimal distance d can correct ⌊(d − 1)/2⌋
errors by associating a received signal with the word in the code that is closest to the signal
in Hamming distance. If C contains N codewords, it is a code of length k, size N , and
minimal distance d over an alphabet of size s, denoted as a (k, N, d)s code. Table 4 provides
a (7, 8, 4)2 code.
When S corresponds to a Galois field, we define a linear code C of length k as a code with
distinct codewords that form a vector subspace of S k . This definition implies that C has size
N = sn for some nonnegative integer n, 0 ≤ n ≤ k, where k is now called the dimension of
the code. A linear code may be characterized by an n × k generator matrix G. The rows of G
form a basis for the code, so that all codewords can be obtained by taking all possible linear
combinations of the rows of G. For any linear code C, there is another linear code called its
dual, and denoted by C ⊥ . This consists of all vectors v ∈ S k such that uv T = 0 for all u ∈ C.
For a (k, sn , d)s linear code C, its dual code C ⊥ is a (k, sk−n , d⊥ )s code, where d⊥ is called
the dual distance of C. Example 4.1 illustrates these concepts. Common examples of linear
codes include Hamming codes, Bose-Chaudhuri-Hocquenghem (BCH) codes, Reed-Solomon
codes, cyclic codes, Golay codes, and Reed-Muller codes. Hedayat et al. (2012) devoted a
chapter discussing the construction of orthogonal arrays using these codes.
Nonlinear codes have been investigated much less than linear codes, but there are families
of nonlinear codes that tend to have better encoding properties than linear codes of the same
size. Nordstrom and Robinson (1967) provided the first nonlinear code now known as the
Nordstrom-Robinson code. It is a (16, 256, 6)2 code with the property that it has dual
distance 6 and offers the advantage over linear codes in that any binary linear code of length
16 with minimal distance 6 can contain at most 128 codewords. For development of families
of nonlinear codes that generalize the Nordstrom-Robinson code we refer to MacWilliams
and Sloane (1977).

Example 4.1. Table 4 lists a (7, 8, 4)2 code. Each row corresponds to a codeword. The
generator matrix of this code is
 
1 1 1 0 1 0 0
G= 0 1 1 1 0 1 0 .
0 0 1 1 1 0 1

11
Its dual code is a (7, 16, 3)2 code which has the generator matrix,
 
1 0 1 1 0 0 0
 0 1 0 1 1 0 0 
 .
 0 0 1 0 1 1 0 
0 0 0 1 0 1 1

Table 4: A (7, 8, 4)2 code

0 0 0 0 0 0 0
1 1 1 0 1 0 0
0 1 1 1 0 1 0
0 0 1 1 1 0 1
1 0 0 1 1 1 0
0 1 0 0 1 1 1
1 0 1 0 0 1 1
1 1 0 1 0 0 1

4.1 Basic Relationship Between Orthogonal Arrays and Codes

We define an orthogonal array OA(N, sk , t) based on the Galois field S = GF (s) to
be linear if its runs are distinct and form a vector space over S k . This is identical to the
definition for a linear code. Consequently, a linear orthogonal array OA(N, sk , t) is a linear
code (k, N, d)s for some d, and vice versa, a linear code (k, N, d)s is a linear orthogonal array
OA(N, sk , t) for some strength t. It turns out that this relationship connects the strength
of a linear orthogonal array to the dual distance of a code. The relationship is formulated
precisely in Theorem 4.1 below, which was first stated by Bose (1961) (see also Hedayat
et al., 2012, for a statement and proof). Through this relationship it is possible to obtain
linear orthogonal arrays from linear codes, and vice versa. Existence results can also be
translated from orthogonal arrays to codes and vice versa.

Theorem 4.1. If C is a (k, N, d)s linear code over S = GF (s) with dual distance d⊥ ,
then the codewords of C form the rows of an OA(N, sk , d⊥ − 1) with entries from GF (s).
Conversely, the rows of a linear OA(N, sk , t) over GF (s) form a (k, N, d)s linear code over
GF (s) with dual distance d⊥ ≥ t + 1. If the orthogonal array has strength t but not t + 1,
then d⊥ = t + 1.

It should be noted that a linear orthogonal array corresponds to a regular fractional

factorial and that the words in the defining relation correspond to the rows in the dual of
this linear orthogonal array.

12
4.2 Rao’s Bound and the Linear Programming Bound
Two very basic and related problems associated with the existence of orthogonal arrays
are the following:
1. For given values of k, s and t, what is the smallest number of runs N for which an
OA(N, sk , t) exists?
2. For given values of N , s and t, what is the largest number of factors k for which an
OA(N, sk , t) exists?
It can be seen that a complete answer to the second question implies a complete answer to
the first question (cf. Hedayat et al., 1992, Chapter 2).
Unfortunately, exact answers to these questions are often unknown. Rao (1947) provides
a lower bound for N in terms of k, s and t that applies to any orthogonal array. The bound
is now also known as Rao’s bound. Implicitly, Rao’s bound also provides an upper bound
for the value of k for given N , s and t. Rao’s bound states that for an OA(N, sk , t) it must
hold that
u  
X k
N≥ (s − 1)i , if t = 2u, and
i=0
i
u    
X k
i k−1
N≥ (s − 1) + (s − 1)u+1 , if t = 2u + 1,
i=0
i u

for u ≥ 0. This result can be understood by counting degrees of freedom for main effects and
interaction effects that can be estimated orthogonally when the orthogonal array is used in
a fractional factorial experiment.
While improvements on Rao’s bound were found for special cases, it wasn’t until the
seminal work by Delsarte (1973) that another general bound was established, which is known
as the linear programming bound. Unlike Rao’s bound, the linear programming bound
provides a bound through computation (by linear programming) and does not provide an
explicit lower bound for N . Theorem 4.2 presents the linear programming bound.
Theorem 4.2. Let NLP (k, d⊥ ) be the solution to the following linear programming problem:
find real numbers A0 , A1 , . . . , Ak to minimize

A0 + A1 + · · · + Ak ,

subject to the constraints

A0 ≥ 1, Ai ≥ 0, 1 ≤ i ≤ k
B0 = 1, Bi ≥ 0, 1 ≤ i ≤ k
B1 = . . . = Bt = 0,

where Bi = kj=0 Aj Pi (j), 0 ≤ i ≤ k, the Pi (j) are the Krawtchouck polynomials Pi (j) =
P
Pi i−r j k−j
, and t = d⊥ − 1. Then, in an OA(N, sk , t), it holds that
r



r=0 (−1) (s − 1) r i−r

N ≥ NLP (k, d⊥ ).

13
 Delsarte (1973) demonstrated that the Rao’s bound follows as a consequence of The-
orem 4.2, establishing that the linear programming bound is always at least as strong as
the Rao’s bound. In many cases, the linear programming bound is significantly stronger.
Hedayat et al. (2012) provided a comparative table highlighting the differences between the
Rao’s and linear programming bounds for binary orthogonal arrays of strength 4 with k
factors. Building on the work of Delsarte (1973), Sloane and Stufken (1996) extended these
results and formulated the linear programming bound for mixed-level orthogonal arrays.
Notably, Table 9.7 in Hedayat et al. (2012) presents the linear programming bounds for
OA(N, 2k1 3k2 , t)’s, demonstrating significant improvements over the Rao’s bound extended
to mixed orthogonal arrays.

5 Connections Between Orthogonal Arrays and Other

Combinatorial Structures
In this section, we will briefly explore connections between orthogonal arrays and other
combinatorial structures, such as mutually orthogonal Latin squares, Hadamard matrices,
incomplete block designs and difference schemes. Interested readers are referred to Chapters
6, 7 and 8 in Hedayat et al. (2012).

5.1 Mutually Orthogonal Latin Squares

Interest by mathematicians in Latin squares and mutually orthogonal Latin squares dates
back to the early 1700s. Their use in statistics was discovered in the 1930s largely due to the
influential work of R.A. Fisher. In particular, Latin squares and mutually orthogonal Latin
squares can play a role in experiments for comparing different treatments in the presence
of multiple blocking variables. A Latin square of order s is an s × s array with entries of a
set of S of size s such that each element of S appears once in every row and column. For
example, the following three arrays are Latin squares of order 4:

0 2 3 1 0 2 3 1 0 2 3 1
3 1 0 2 1 3 2 0 2 0 1 3
1 3 2 0 2 0 1 3 3 1 0 2
2 0 1 3 3 1 0 2 1 3 2 0
Two Latin squares of order s are said to be orthogonal to each other if, when superimposed
on each other, each of the s2 pairs (i, j) appears in exactly one cell, for 1 ≤ i, j ≤ s. It can
be verified that any two of the above three Latin squares of order 4 are orthogonal to each
other. A set of mutually orthogonal Latin squares is a collection of Latin squares of order s in
which any pair is orthogonal. The existence and construction of mutually orthogonal Latin
squares are discussed in Hedayat et al. (2012). The primary connection between mutually
orthogonal Latin squares and orthogonal arrays is summarized in Theorem 5.1. The proof
can be found in Hedayat et al. (2012) and Cheng (2016).

Theorem 5.1. A set of k mutually orthogonal s × s Latin squares is equivalent to an

OA(s2 , sk+2 , 2).

14
 For example, the three mutually orthogonal Latin squares of order 4 are equivalent to
an OA(16, 45 , 2). It is known that an upper bound on the number of mutually orthogonal
Latin squares of order s is s − 1 and this upper bound can be achieved when s is a prime or
prime power. Even with today’s computing power, the problem of obtaining the maximum
possible number of mutually orthogonal Latin squares for other values of s remains, with
few exceptions, a challenging undertaking. The study of this problem was pioneered by Bose
et al. (1960) and Wilson (1974), and pursued by many others. Interested readers are referred
to the brief survey provided by Colbourn and Dinitz (2001). We note that although the
connection between mutually orthogonal Latin squares and orthogonal arrays is fascinating,
the use of mutually orthogonal Latin squares for constructing new orthogonal arrays may be
limited. Hedayat et al. (2012) described a number of interesting related research problems
to be addressed.

5.2 Hadamard Matrices

The concept of a Hadamard matrix originated in 1893 from the work of the French
mathematician Jacques Salomon Hadamard. The study of Hadamard matrices expanded
significantly in the 20th century, particularly with their applications in coding theory, signal
processing, and design of experiments. Formally, a Hadamard matrix is a square matrix
whose entries are either 1 or −1 and whose rows (and, hence, columns) are mutually orthog-
onal. Mathematically, a Hadamard matrix H satisfies,

HH T = N IN

where N is the order of the matrix and IN is the N × N identity matrix. Theorem 5.2
establishes a connection between Hadamard matrices and orthogonal arrays. Multiplication
of entire rows or columns of a Hadamard matrix by −1 will again result in a Hadamard
matrix, so that for every order N for which a Hadamard matrix exists, there is one with all
entries in the first row or column equal to 1.

Theorem 5.2. Suppose H is a Hadamard matrix of order N > 2 such that all entries
in its first column are 1. If the first column of H is removed, the resulting matrix is an
OA(N, 2N −1 , 2). Conversely, appending a column of all 1’s to an OA(N, 2N −1 , 2) produces a
Hadamard matrix of order N .

The proof of Theorem 5.2 can be found in Hedayat et al. (2012) and Cheng (2016). Since
the existence of an OA(N, 2N −1 , 2) is equivalent to the existence of an OA(2N, 2N , 3), the
existence of the latter orthogonal array is also equivalent to the existence of a Hadamard
matrix of order N (cf. Hedayat and Wallis, 1978; Hedayat et al., 2012). If there exists a
Hadamard matrix of order N > 2, then N must be a multiple of 4. According to the
Hadamard conjecture, Hadamard matrices exist for all orders that are multiples of 4, but
their existence has not been proven for all such orders. The smallest multiple of 4 for which
no Hadamard matrix has been found is 668, despite extensive computational searches. Given
the equivalence between Hadamard matrices and orthogonal arrays, it is clear that two-level
orthogonal arrays can be constructed by selecting specific columns from Hadamard matrices.
Orthogonal arrays obtained from Hadamard matrices, also called Hadamard designs, are not

15
regular when the order N of the Hadamard matrix is not a power of 2; when N is a power
of 2, they may be either regular or non-regular. This flexibility allows Hadamard designs
to accommodate a wider range of run sizes compared to regular designs. In their seminal
work, Plackett and Burman (1946) introduced the use of Hadamard designs in factorial
experiments. The Hadamard designs described in their paper are now commonly known as
Plackett-Burman designs.
Several well-established methods have been introduced to construct Hadamard matrices.
These include the Sylvester construction, Paley construction, Williamson construction, con-
ference matrices, and algebraic approaches involving group theory or combinatorial designs
(Sylvester, 1867; Paley, 1933; Williamson, 1944; Hall, 1998; Georgiou et al., 2003; Kline,
2019). For small orders, Hadamard matrices can be found through exhaustive search or
optimization algorithms. Several researchers have conducted theoretical investigations into
identifying suitable Hadamard matrices for constructing two-level non-regular designs (Shi
and Tang, 2018; Chen et al., 2023). In addition, Hadamard matrices have been used to
construct mixed orthogonal arrays through clever use of properties of Hadamard matrices
(Dey and Ramakrishna, 1977; Chacko et al., 1979; Chacko and Dey, 1981; Agrawal and Dey,
1982; Cheng, 1989; Wang, 1990). See Section 4.3 of Dey and Mukerjee (2009) for detailed
construction methods and results.

5.3 Incomplete Block Designs

Blocking is one of the three fundamental principles in the design of experiments. It is
an effective strategy for explaining response variability by controlling for known nuisance
factors, thereby increasing the precision of treatment effect estimation (Wu and Hamada,
2011). The simplest and most frequently used block designs are randomized complete block
designs. However, blocking should be based on differences in experimental units, and it
may not always be feasible to include all treatments within a single block. Incomplete block
designs may then be a sensible choice. Yates introduced balanced incomplete block designs, for
which v treatments are arranged in b blocks of k experimental units, k < v, each treatment
occurring in r blocks, and any two treatments occurring together in λ blocks. Each balanced
incomplete block design has therefore the parameters: v (the number of treatments); b (the
number of blocks); r (the number of blocks in which each treatment appears); k (the number
of treatments in each block); λ (the number of blocks in which any pair of treatments appears
together). Necessary conditions for the existence of a balanced incomplete block design are
that the parameters satisfy bk = vr, λ(v − 1) = r(k − 1) and b ≥ v. If b = v, then r = k and
the design is called a symmetric balanced incomplete block design (cf. Lander, 1983).
The connection between orthogonal arrays and incomplete block designs is partly due to
both being connected to Hadamard matrices (Hedayat and Wallis, 1978). One connection
between orthogonal arrays and incomplete block designs is stated in the following theorem.
Theorem 5.3. The existence of an orthogonal array OA(N, 2N −1 , 2) implies the existence
of a symmetric balanced incomplete block design with v = b = N − 1, r = k = N/2 − 1, and
λ = N/4 − 1.
To see the validity of Theorem 5.3, we represent an incomplete block design through
its incidence matrix, say M = (mij ). This is a v × b matrix with entries 0 and 1, and

16
with mij = 1 if and only if the ith treatment appears in the jth block. Starting from the
orthogonal array, we obtain a Hadamard matrix of order N (see Theorem 5.2), and convert it
to one in which all entries in both the first column and row are 1. The matrix that remains
after deleting the first column and row is the incidence matrix for a symmetric balanced
incomplete block design with the parameters as in Theorem 5.3. The converse also holds, so
that the existence of this symmetric balanced incomplete block design is equivalent to the
existence of an orthogonal array OA(N, 2N −1 , 2).
There are multiple other block designs that are related to the incomplete block design in
Theorem 5.3. For example, by replacing each block by its complement in the treatment set
we obtain a symmetric balanced incomplete block design with v = b = N − 1, r = k = N/2,
and λ = N/4. As another example, starting from the symmetric block design in Theorem 5.3,
upon deleting an entire block and deleting the treatments in that block also from all of the
other blocks, we obtain the so-called residual design, which is a balanced incomplete block
design with v = N/2, b = N − 2, r = N/2 − 1, k = N/4, and λ = N/4 − 1. Additional
connections can be found in Hedayat and Wallis (1978) and Raghavarao (1988).

5.4 Difference Schemes

Difference schemes were first defined by Bose and Bush (1952). We review their con-
cept and their connection with orthogonal arrays. More importantly, they are a simple yet
powerful tool for constructing orthogonal arrays as we will see in Section 6.
An r × c array with entries from a finite Abelian group containing s entries is called a
difference scheme if each vector difference between any two distinct columns of the array
contains every element from the group equally often (Bose and Bush, 1952). Such an array
is denoted by D(r, c, s). Table 5 displays a D(9, 9, 3). Hedayat et al. (2012) discussed the
existence, construction, properties, and generalizations of difference schemes. Statistical
Analysis Software (SAS) provides a library of difference schemes with a large number of
different values for r and c and 3 ≤ s ≤ 22. The library can be accessed at the webpage
SAS (2019). In this definition, difference schemes are of strength two. Hedayat et al. (1996)
introduced difference schemes of strength t for the integer t ≥ 2.
Every orthogonal array OA(N, sk , 2) based on an Abelian group is a difference scheme
D(N, k, s), but typically not a very interesting one. When looking for a difference scheme
D(r, c, s) one is typically interested in obtaining an array with the largest possible value of
c for given values of r and s, and the additional structure required by an orthogonal array
will fail to provide an array with the maximum possible value for c.
For s = 2, using the multiplicative Abelian group consisting of 1 and −1, a Hadamard
matrix of order N is a difference scheme D(N, N, 2) and, conversely, such a difference scheme
can be converted to a Hadamard matrix of order N . Thus, through Hadamard matrices,
this provides a connection between orthogonal arrays and difference schemes. While this
is a simple and interesting connection, what is perhaps more interesting is that difference
schemes can be used to construct larger arrays that are orthogonal arrays. For example, as we
will see in Subsection 6.1, the existence of a difference scheme D(λs, k, s) implies that of an
orthogonal array OA(λs2 , sk+1 , 2) and a mixed orthogonal array OA(λs2 , (λs)1 sk , 2). These
orthogonal arrays can therefore be presented via the smaller difference schemes, and exploring
the construction of the difference schemes, whether through combinatorial or computational

17
methods, can be easier than direct constructions of orthogonal arrays. On the flip side, not
every orthogonal array can be constructed through a difference scheme.
If s is a prime or prime power, then D(s, s, s) can be obtained based on the Galois
field GF (s), so that can we obtain an OA(s2 , ss+1 , 2), which provides equality in the Rao’s
bound. Hedayat et al. (2012) studied an extension to difference schemes of strength t,
t ≥ 3, established the existence of such arrays, and showed how they could be used for the
construction of orthogonal arrays of strength t. Difference schemes have also been used for
the construction of mixed orthogonal arrays (cf. Wang and Wu, 1991).

Table 5: A difference scheme D(9, 9, 3)

0 0 0 0 0 0 0 0 0
0 1 2 0 1 2 0 1 2
0 2 1 0 2 1 0 2 1
0 0 0 2 2 2 1 1 1
0 1 2 2 0 1 1 2 0
0 2 1 2 1 0 1 0 2
0 0 0 1 1 1 2 2 2
0 1 2 1 2 0 2 0 1
0 2 1 1 0 2 2 1 0

6 Selected Construction Methods for Orthogonal Ar-

rays
This section provides a concise review of selected construction methods for orthogonal
arrays. The literature on methods for constructing orthogonal arrays is extensive, encom-
passing a wide range of approaches. These include direct construction methods and meth-
ods based on other mathematical structures such as Hadamard matrices, finite fields, Latin
squares, difference sets, cyclic groups, and projective planes. Recursive methods, such as
those utilizing Kronecker products or Kronecker sums, are also widely used. In addition,
intelligent computer-based searches have become an important tool for generating orthog-
onal arrays tailored to specific requirements. Here, we review several general and widely
recognized methods.

6.1 Constructions Based on Difference Schemes

Theorem 6.1 indicates that a difference scheme can be converted into an orthogonal array
straightforwardly.

18
Theorem 6.1. If D is a difference scheme D(r, c, s) based on an Abelian group A =
{σ0 , σ1 , . . . , σs−1 } with binary operation +, then
 
D0
 D1 
A =  .. 
 
 . 
Ds−1

is an OA(rs, sc , 2), where Di is obtained from D by adding σi to each of its entries.

The orthogonal array in Theorem 6.1 can be extended by at least one column. For this
we use the Kronecker product of an n1 × m1 matrix A = (aij ) and an n2 × m2 matrix B that
are based on the same group with binary operation +. Their Kronecker product is denoted
by A⊗B, and is defined as the (n1 n2 )×(m1 m2 ) block matrix with the block in location (i, j)
equal to aij + B, 1 ≤ i ≤ n1 , 1 ≤ j ≤ m1 . If σ0 denotes the identity element of the group,
then the column that can be added to the array in Theorem 6.1 is σ0 1r/s ⊗ (σ0 , ..., σs−1 )T ,
where σ0 1r/s denotes the r/s × 1 vector with every entry equal to σ0 and (σ0 , ..., σs−1 )T is the
s×1 vector that contains every element of A. Thus, a difference scheme D(r, c, s) implies the
existence of an orthogonal array OA(rs, sc+1 , 2). This explains the claim in Subsection 6.1
that a difference scheme D(s, s, s) results in an OA(s2 , ss+1 , 2).
Kronecker products are a workhorse in constructing orthogonal arrays from difference
schemes. Another use is outlined in Theorem 6.2. This theorem states that a difference
scheme, combined with an existing orthogonal array, can be used to generate an orthogonal
array with a larger run size and more columns.
Theorem 6.2. If D is a difference scheme D(r, c, s) and B is an OA(N, sk , 2), both based
on the same Abelian group, then the array

A=B⊗D

is an orthogonal array OA(N r, skc , 2).

Additional methods for using difference schemes in the construction of orthogonal arrays
can be found in Chapter 6 of Hedayat et al. (2012). The concept was generalized by Seiden
(1954), who introduced difference schemes of strength t. Subsequently, Hedayat et al. (1996)
examined their existence and construction, providing methods to construct orthogonal arrays
of strength t from these schemes for t > 2.

6.2 Constructions Based on Hadamard Matrices

In Theorem 5.2 we already saw how a Hadamard matrix of order N can be used to
construct an orthogonal array OA(N, 2N −1 , 2). In addition, if H is a Hadamard matrix,
then the array  
H
−H
is an orthogonal array OA(2N, 2N , 3). Both arrays give equality in the Rao’s bound. Thus,
the construction of these orthogonal arrays boils down to the construction of Hadamard

19
matrices. For the same order, different Hadamard matrices can result in orthogonal arrays
with different properties. For example, if N is a power of 2, we could obtain a regular or
non-regular fractional factorial depending on the Hadamard matrix that we use.
The simplest construction method for Hadamard matrices of order N = 2m , m ≥ 2, is
the Sylvester method, which is based on repeated use of the Kronecker product (with the
binary operation being multiplication), as defined in Subsection 6.1. This result is formally
stated in Theorem 6.3.
Theorem 6.3. Let H be an array obtained by the (m − 1)-fold Kronecker product, m ≥ 2,
   
1 −1 1 −1
H= ⊗ ··· ⊗ .
1 1 1 1

Then H is a Hadamard matrix of order 2m .

More generally, for Hadamard matrices H1 and H2 of order N1 and N2 , their Kronecker
product H1 ⊗ H2 is a Hadamard matrix of order N1 N2 .
Based on Galois fields, Paley (1933) developed two methods of construction for infinite
families of Hadamard matrices. The first Paley construction works when q = N − 1 is an
odd prime power, where N , a multiple of 4, is the order of the Hadamard matrix to be
constructed. Denote the elements of a Galois field GF (q) by α1 = 0, α2 , . . . , αq , and define
a function χ : GF (q) → {0, 1, −1} as

1,
 if β = x2 for some x ∈ GF (q);
χ(β) = 0, if β = 0;

−1, otherwise.


Let A = (aij ) be the q × q matrix with aij = χ(αi − αj ) for i, j = 1, 2, . . . , q, and define
 
1 −1Tq
H= (3)
1q A + Iq

where 1q is a column of q ones and Iq is the q × q identity matrix. The matrix H in (3) is a
Hadamard matrix of order N (cf. Hedayat et al., 2012; Cheng, 2016). The OA(N, 2N −1 , 2)
constructed from H in (3) is also known as the Paley design of order N .
The second Paley construction is for Hadamard matrices of order N = 2q + 2, where q is
an odd prime power and q − 1 is a multiple of 4. Using the matrix A as defined for the first
Paley construction, define

1Tq 1Tq
 
1 −1
 1q A + Iq 1q A − Iq 
H= T
. (4)
−1 1q −1 −1Tq 
1q A − 1q −1q −A − Iq

Then H in (4) is a Hadamard matrix of order N . Multiplying the (q + 2)nd row of H in (4)
by −1 gives a Hadamard matrix with first column equal to 1, which can be used to construct
an orthogonal array OA(N, 2N −1 , 2). Chen et al. (2023) conducted a comprehensive study of

20
orthogonal arrays constructed from Hadamard matrices obtained by the Paley constructions
in terms of generalized resolution, projectivity, and hidden projection property.

6.3 The Rao-Hamming Construction

Multiple methods of construction for orthogonal arrays make use of Galois fields or finite
geometries. One of the earlier and simpler methods to describe was introduced independently
by Rao (1947, 1949) and Hamming (1950), and the corresponding orthogonal arrays were
named Rao-Hamming arrays by Hedayat et al. (2012). There are actually multiple methods
of construction for the Rao-Hamming arrays (see Hedayat et al., 2012, for three methods),
and we present one of these methods here. Rao-Hamming arrays are orthogonal arrays
OA(sn , sk , 2) where s is a prime power, k = (sn − 1)/(s − 1), and n ≥ 2.
Construction of Rao-Hamming Arrays: Construct an sn × n array with the rows con-
sisting of all n-tuples based on GF (s). Denote the columns of this array as C1 , . . . , Cn . The
columns of the orthogonal array are obtained as

z1 C1 + · · · + zn Cn = [C1 . . . Cn ]z,

using all n-tuples z = (z1 , . . . , zn )T based on GF (s), with zi ̸= 0 for at least one i and the first
non-zero zi equal to the unit element 1. There are (sn − 1)/(s − 1) such n-tuples, resulting
in the required number of columns. Theorem 6.4 summarizes the result of this construction.

Theorem 6.4. If s is a prime power, then an OA(sn , sk , 2) with k = (sn − 1)/(s − 1) exists
whenever n ≥ 2.

The Rao-Hamming arrays achieve equality in the Rao’s bound and can also be presented
as linear orthogonal arrays. The number (sn −1)/(s−1) corresponds to the number of points
in the finite projective geometry P G(n − 1, s), and the n-tuples used in the proof can be
thought of as the points in this geometry.

6.4 Recursive Constructions

Recursive methods are essential tools for constructing orthogonal arrays. In Subsec-
tion 6.1, we presented a recursive approach based on difference schemes. This subsection
examines a recent recursive technique proposed by He et al. (2022) for constructing larger
orthogonal arrays from smaller ones with a large number of factors.
We first define an operator based on the Kronecker product. For an n1 × m1 matrix A
and n2 × m2 matrix B, n1 ≤ n2 , and the rows of B partitioned into submatrices B1 , ..., Bn1 ,
we define the generalized Kronecker product as
   
B1 a1 ⊗ B1
A⃝ ∗ B =A⃝ ∗  ...  = [ai ⊗ Bi ]1≤i≤n1 =  ..
, (5)
   
.
Bn1 an1 ⊗ Bn1

where ai denotes the ith row of A and the operator ⊗ represents the usual Kronecker product
defined in Subsection 6.1. Note that A ⃝∗ B depends on the partition of the rows of B, that

21
A⃝∗ B has the same number of runs as B, and that the number of factors in A ⃝
∗ B is m1 m2 .
T
For example, with addition modulo 3, let A = (0, 1, 2) and
   
0 0 0 0 0 0 0 0
 0 1 1 2   0 1 1 2 
   
 0 2 2 1   0 2 2 1 
     
B1  1 0 1 1 
 
 2 1 2 2 
 
B =  B2  =   1 1 2 0 , then we have A ⃝
 ∗ B=  2 2 0 1 .

B3  1 2 0 2 
 
 2 0 1 0 
 
 2 0 2 2   1 2 1 1 
   
 2 1 0 1   1 0 2 0 
2 2 1 0 1 1 0 2
The construction in He et al. (2022) uses an n1 × m1 matrix A with rows a1 , . . . , an1 and
n1 matrices B1 , . . . , Bn1 each of size (n2 /n1 ) × m2 that partition the rows of B, with entries
for all matrices from the Galois field GF (s) = {α0 = 0, α1 , . . . , αs−1 } for a prime power s.
They define s+1 arrays D1 , D2 , . . . , Ds+1 as follows, where the Kronecker products are based
on addition in GF (s):
(i) For αg ∈ GF (s) and g = 1, . . . , s − 1, define
 
a1 ⊗ (αg ∗ B1 )
∗ (αg ∗ B) = 
Dg = A ⃝ ..
,
 
.
an1 ⊗ (αg ∗ Bn1 )

where ∗ represents multiplication in GF (s).

(ii) For g = s,  
0 ⊗ B1
Ds = 0n1 ∗ B=
⃝ ..
,
 
.
0 ⊗ Bn1
where 0 denotes the zero element in GF (s) and 0n1 is the n1 × 1 vector of 0’s.
(iii) For g = s + 1,
Ds+1 = A ⊗ 0(n2 /n1 ) ,
where 0(n2 /n1 ) is the (n2 /n1 ) × 1 vector of 0’s.
Now define the array
E = [D1 , D2 , . . . , Ds+1 ] . (6)
For 1 ≤ g ≤ s − 1, Dg is an n2 × (m1 m2 ), while Ds and Ds+1 are n2 × m2 and n2 × m1 arrays,
respectively. Consequently, E is an n2 × [(s − 1)m1 m2 + m1 + m2 ] array.
The key result in He et al. (2022) is summarized in the following theorem.
Theorem 6.5. For a prime power s, if A is an OA(n1 , s1 , 1) or an OA(n1 , sm1 , 2) for
m1 > 1, and each Bi , i = 1, . . . , n1 , is an OA((n2 /n1 ), sm2 , 2), then the array E in (6) is an
OA(n2 , sk , 2), where k = (s − 1)m1 m2 + m1 + m2 .

22
 Thus, for a prime power s, if A and the Bi ’s in this construction are s-level orthogonal
arrays of strength two, then so is the resulting array E. This construction uncovers the hidden
structure of many existing fixed-level orthogonal arrays, often producing arrays with more
factors then previously known for a given run size. In addition, He et al. (2022) explored
how the construction in (6) can be leveraged in the construction of orthogonal arrays of
strength three, as well as other types of orthogonal arrays, such as resolvable orthogonal
arrays, balanced sliced orthogonal arrays, and nested orthogonal arrays. These variants of
orthogonal arrays, which will be further discussed in Section 8, offer diverse applications in
experimental design and enhance flexibility for specific design goals.

7 Distinguishing Between Orthogonal Arrays With the

Same Parameters
Two orthogonal arrays OA(N, sk , t) with the same values for the parameters N , k, s and
t may have different properties with respect to statistically meaningful criteria when using
the orthogonal arrays as fractional factorials. Often, if there is an orthogonal array for a set
of parameters, there are many such arrays, leading to the question of which one to use in an
experiment. Some are isomorphic and others are non-isomorphic. Two orthogonal arrays are
(algebraically) isomorphic if one can be obtained from the other by row permutations, column
permutations, and/or relabeling the levels of one or more factors. Two isomorphic orthogonal
arrays have similar properties with respect to some statistical criteria, but not with respect
to others. Nevertheless, studying isomorphism classes is valuable for the enumeration of all
possible orthogonal arrays for a given set of parameters. This is a challenging combinatorial
problem due to the huge number of possible arrays. For any given array, there are after all
N ! × k! × (s!)k permutations to consider, although this number can be reduced since many
of these permutations will result in the same array. With increasing values for N and k,
determining the isomorphism classes is however an NP-hard problem.
The study of isomorphism of orthogonal arrays has resulted in an extensive literature.
Draper and Mitchell (1967) compared word-length patterns of designs to determine their
isomorphism. However, two designs with the same word-length pattern could be non-
isomorphic. Draper and Mitchell (1970) proposed a more sensitive test for isomorphism
using a “letter pattern comparison.” However, Chen and Lin (1991) gave two non-isomorphic
231−15 designs with identical letter pattern matrices and thus demonstrated that the letter
pattern also does not uniquely determine a fractional factorial. Lin and Sitter (2008) subse-
quently showed that there are many such cases. Chen (1992) discussed the isomorphism of
2m−p fractional factorials in terms of the existence of a relabeling map between two frequency
vectors together with an appropriately defined matrix. Using this frequency representation,
Chen (1992) proved that the word-length pattern uniquely determines any 2m−p fractional
factorial with p = 1 or two and any minimum aberration 2m−p fractional factorial when
p = 3 or 4. Clark and Dean (2001) introduced a method of determining isomorphism of any
two (regular and non-regular) fractional factorials by examining the Hamming distance ma-
trices of their projection designs. An algorithm was provided for checking the isomorphism
of fractional factorials when all the factors have two levels which saves considerable time

23
for detecting non-isomorphic arrays. Other work on determining isomorphism of orthogonal
arrays focuses on using different properties of the arrays such as the centered L2 discrepancy
(Ma et al., 2001), minimal column base (Sun et al., 2002), singular value decomposition
(Katsaounis et al., 2013), and degree of isomorphism (Weng et al., 2023). Katsaounis and
Dean (2008) provided a survey and evaluation of methods for determination of isomorphism
of fractional factorials.
While the criteria mentioned in the previous paragraph can be helpful in deciding whether
two arrays are isomorphic, they don’t solve the enumeration problem of non-isomorphic or-
thogonal arrays for given parameters. Nonetheless, significant progress has been made on
this enumeration problem. The foundational work of Chen et al. (1993) provided a compre-
hensive catalog of all possible 2m−p fractional factorials of size 16 and 32 and all resolution
four (or higher) fractions of size 64. Block and Mee (2005) extended the enumeration to
two-level 128-run resolution IV designs, offering an expanded set of designs for experimen-
tal applications requiring a larger number of runs. Xu (2005) proposed methods based on
coding theory to efficiently classify and rank fractional factorials. This approach facilitated
the enumeration of fractional factorials with 27, 81, and 243 runs at resolution IV or higher,
as well as 729-run arrays at resolution V or higher. Stufken and Tang (2007) provided a
complete solution to enumerating non-isomorphic two-level orthogonal arrays of strength d
with d + 2 constraints for any d and any run size N = λ2d . Their work represents a signifi-
cant milestone, as it systematically classified all structurally unique designs within this class.
Schoen et al. (2010) introduced a minimum complete set algorithm for generating catalogs
of non-isomorphic symmetric and mixed-level orthogonal arrays with specified strength, run
sizes, number of factors, and number of factor levels. The cases include all mixed-level
strength two orthogonal arrays with the run size N ≤ 28, all symmetric strength two or-
thogonal arrays with N ≤ 27, all OA(28, 2a , 2) with a ≤ 7, all strength three orthogonal
arrays with N ≤ 64 except OA(56, 2a , 3), OA(64, 2a , 3), and OA(64, 41 2a , 3), and all strength
four orthogonal arrays with N ≤ 168 except OA(160, 2a , 4). For OA(56, 2a , 3), OA(64, 2a , 3),
OA(64, 41 2a , 3) and OA(160, 2a , 4), the catalog of non-isomorphic designs are obtained for
a ≤ 8, 7, 6, 8, respectively. The accompanying Python-based software provides an exten-
sive collection of catalogs containing non-isomorphic orthogonal arrays covering more cases
than the results provided by Schoen et al. (2010). More recent studies, such as Vazquez
and Xu (2019); Bohyn et al. (2023); Eendebak et al. (2023), have focused on constructing
non-isomorphic orthogonal arrays for specific levels and larger run sizes.
Another critical problem in the study of orthogonal arrays is ranking arrays with specific
parameters to identify the optimal ones according to a meaningful criterion. Early founda-
tional works include Box and Hunter (1961) and Fries and Hunter (1980), which introduced
the optimality criteria of resolution and minimum aberration, respectively, for regular frac-
tional factorial designs. This line of research gained momentum in the late 20th century, and
numerous studies have been published since that time. Proposed optimality criteria include
minimum G-aberration (Deng and Tang, 1999a), minimum G2 -aberration (Deng and Tang,
1999b), minimum moment aberration (Xu, 2003), various uniformity measures (Fang and
Mukerjee, 2000), and estimation capacity (Cheng et al., 1999). For a comprehensive review,
see Xu et al. (2009) and Cheng (2016).

24
8 Selected Recent Developments
In recent years, practical applications in science, engineering, and technology have led to
multiple innovations in research related to orthogonal arrays. This section highlights some of
the most notable developments that have drawn particular interest within the experimental
design community. These include sliced orthogonal arrays, nested orthogonal arrays, strong
orthogonal arrays, and grouped orthogonal arrays. The following subsections provide an
overview of these concepts and key advancements in their development.

8.1 Sliced Orthogonal Arrays

Sliced orthogonal arrays were developed to address the needs of computer experiments in-
volving both quantitative and qualitative input variables, a scenario that frequently occurred
in real-world applications. For example, a data center experiment might include qualitative
factors such as diffuser height and the location of hot-air return vents (Schmidt et al., 2005).
Similarly, computer experiments in marketing and social sciences often involve qualitative
factors such as education level, race, and social background. Qian and Wu (2009) was
the first to introduce the concept of sliced orthogonal arrays to tackle experimental design
challenges in such mixed-input computer experiments.
To define sliced orthogonal arrays, we first define the concept of a level-collapsing pro-
jection. A mapping δ(·) is called a level-collapsing projection if the mapping is from a set S
with s elements to a set of s0 elements and satisfies (a) S can be partitioned into s0 subsets
S1 , . . . , Ss0 , with each having s/s0 elements, and (b) for any two elements x ∈ Si , y ∈ Sj ,
δ(x) = δ(y) for i = j, and δ(x) ̸= δ(y) otherwise. Table 6 provides an example of the
level-collapsing projection δ(0) = δ(1) = 0 and δ(2) = δ(3) = 1.
Definition 8.1. An OA(N, sk , 2) D is called a sliced orthogonal array if the N rows of D can
be partitioned into v subarrays D1 , D2 , . . . , Dv such that each Di becomes an OA(N0 , sk0 , 2)
with N0 = N/v after the s levels in each column of D are collapsed to s0 levels according to
some level-collapsing projection.
We denote such an array by SOA(N, sk , 2; v, s0 ). Furthermore, if each column in each
slice Di of a sliced orthogonal array is balanced, that is, has an equal occurrence of each
of the s levels, it is called a balanced sliced orthogonal array (Ai et al., 2014). We use
BSOA(N, sk , 2; v, s0 ) to denote such an array. The sliced orthogonal array D in Table 6 is
a BSOA(16, 43 , 2; 4, 2), where rows 1–4, 5–8, 9–12, 13–16 correspond to D1 , D2 , D3 and D4 ,
respectively.
A sliced orthogonal array can be used to construct sliced space-filling designs, which
are used to choose inputs for computer experiments with both qualitative and quantitative
inputs. Suppose we have v level combinations of the qualitative factors, and k quantitative
input variables in the experiment. We will use a sliced orthogonal array with v slices and
k factors. The procedure of constructing sliced space-filling designs consists of two main
steps. First, for the quantitative factors, a Latin hypercube design is generated using a
sliced orthogonal array in the same way as done for generating an orthogonal array-based
Latin hypercube (see Subsection 3.1), with the rows of the resulting Latin hypercube designs
being partitioned into v slices corresponding to the partition in the sliced orthogonal array.

25
Table 6: An orthogonal array D = OA(16, 43 , 2) and the level-collapsing projection δ: D is
also a BSOA(16, 43 , 2; 4, 2) and an N OA(16, 43 , 2; 4, 2)

D δ(D)
0 0 0 0 0 0
2 1 3 1 0 1
1 3 2 0 1 1
3 2 1 1 1 0
2 2 2 1 1 1
0 3 1 0 1 0
3 1 0 1 0 0
1 0 3 0 0 1
1 1 1 0 0 0
3 0 2 1 0 1
0 2 3 0 1 1
2 3 0 1 1 0
3 3 3 1 1 1
1 2 0 0 1 0
2 0 1 1 0 0
0 1 2 0 0 1

Second, each of these v slices is then associated with a level combination of the qualitative
factors. This approach ensures an efficient and balanced design for experiments with mixed
inputs. For example. the sliced orthogonal array in Table 6 can be used to generate sliced
space-filling designs for a computer experiment with three quantitative input variables and
four level combinations for qualitative factors, where the four level combinations could form
a full factorial for two qualitative factors or a fractional factorial for three qualitative factors.
Qian and Wu (2009) introduced several methods for constructing sliced orthogonal arrays,
including based on the Rao-Hamming method. Building on this pioneering work, additional
construction methods have been developed, as detailed in, among others, Xu et al. (2011);
Ai et al. (2014); Li et al. (2015); Hwang et al. (2016); Zhang et al. (2018); He (2019); Tsai
(2022); He et al. (2022); Pang and Zhu (2024), to construct sliced orthogonal arrays with
more flexible parameters.

8.2 Nested Orthogonal Arrays

Qian et al. (2009) introduced the concept of nested orthogonal arrays to construct nested
space-filling designs, which are used to select inputs for computer experiments at two levels
of accuracy or fidelity. Such experiments are commonly seen in practice when computer
models can be run with varying degrees of sophistication, resulting in different computational
times. For instance, Qian and Wu (2008) studied computer simulations for a heat exchanger
in an electronic cooling application. In this case, two computer codes were employed to
simulate linear cellular alloys for electronic cooling systems. One code used finite-element
analysis, providing high accuracy but requiring more computational time, while the other

26
relied on the finite difference method, yielding a faster but less precise approximation. The
differences in numerical methods and grid resolution created a trade-off between accuracy
and computational efficiency, highlighting the need for nested designs.
Definition 8.2. An OA(N, sk , 2) D is called a nested orthogonal array if D contains a
subarray D0 that, after applying a specific level-collapsing projection to each column of D,
becomes an OA(N0 , sk0 , 2).
We denote such an array by N OA(N, sk , 2; N0 , s0 ). For illustration, consider the BSOA
(16, 3, 4, 2; 4, 2) D in Table 6. This array is also an N OA(16, 43 , 2; 4, 2) where any Di for
i = 1, 2, 3, 4, can serve as the subarray D0 (rows 1–4, 5–8, 9–12, 13–16 correspond to D1 ,
D2 , D3 and D4 , respectively). In fact, a sliced orthogonal array is a nested orthogonal array,
however the reverse does not hold.
Mukerjee et al. (2008) and Wang and Yin (2013) explored the existence of nested or-
thogonal arrays. Building on the work of Qian et al. (2009), many studies have investigated
methods for constructing nested orthogonal arrays, including Qian et al. (2009); Dey (2010);
He and Qian (2011); Dey (2012); Wang and Li (2013); Sun et al. (2014); Zhang et al. (2018,
2019); He et al. (2022); Pang and Zhu (2024). In general, nested orthogonal arrays can
be constructed through various direct methods, such as the Rao-Hamming method, Bush’s
method, or by leveraging structures like nested difference matrices, Hadamard matrices, re-
solvable orthogonal arrays, zero-sum arrays, and operators like the Kronecker product and
subgroup projection. Recursive methods are also used, offering an iterative approach to
generate larger nested orthogonal arrays from smaller ones. These construction methods
aim to provide greater flexibility in terms of run sizes and the number of levels for each
factor, enabling the development of nested designs tailored to a wide range of experimental
requirements.

8.3 Strong Orthogonal Arrays

He and Tang (2013) introduced the concept of strong orthogonal arrays, which are used
to generate space-filling designs with enhanced projection and space-filling properties. These
arrays have broad applications in optimization, prediction, and sensitivity analysis of com-
plex systems, where effective exploration of the design space is critical. The introduction of
strong orthogonal arrays by He and Tang (2013) has received considerable interest from both
researchers and practitioners, leading to advancements in both design theory and practical
applications.
Definition 8.3. An N × k array with entries from {0, 1, . . . , st − 1} is called a strong or-
thogonal array of size N , k factors, st levels, and strength t if any subarray of g columns
for any g with 1 ≤ g ≤ t can be collapsed into an OA(N, su1 su2 · · · sug , g) for any positive
integers u1 , . . . , ug with u1 + u2 + · · · + ug = t where collapsing st levels into suj levels is done
through the map a → [a/st−uj ], with [x] denoting the largest integer not exceeding x.
Note that in Definition 8.3, the notation su1 stands for a single factor at su1 levels rather
than for u1 factors at s levels. Similarly for the other g − 1 factors. We use SOA(N, (st )k , t)
to denote such a strong orthogonal array. Note that, despite sharing the same acronym as
sliced orthogonal arrays, the notation differs in both the type and number of parameters.
Table 7 presents an SOA(8, 83 , 3). This array has the following properties:
27
 (i) The array becomes an OA(8, 23 , 3) after the eight levels are collapsed into two levels
according to [a/4] = 0 for a = 0, 1, 2, 3 and [a/4] = 1 for a = 4, 5, 6, 7, where [x] denotes
the largest integer not exceeding x;

(ii) Any two columns of the array become an OA(8, 21 41 , 2) or OA(8, 41 21 , 2) after the
eight levels of one column are collapsed into two levels by [a/4] and the eight levels of
the other column are collapsed into four levels by [a/2];

(iii) Any column of the array is an OA(8, 81 , 1).

Table 7: An SOA(8, 83 , 3)

0 0 0
2 3 6
3 6 2
1 5 4
6 2 3
4 1 5
5 4 1
7 7 7

Clearly, N = λst must hold for some integer λ. As was the case for an orthogonal array,
λ is called the index of the strong orthogonal array. He and Tang (2013) noted that the
definition of a strong orthogonal array is motivated by nets for quasi-Monte Carlo point sets,
and nets are a special case of strong orthogonal arrays. Their connection is documented in
He and Tang (2013). If λ = sw for some integer w, then the existence of an SOA(λst , (st )k , t)
is equivalent to that of a (w, m, k)-net with base s where m = w + t. Nets are defined with
the restriction that the index is a power of s while strong orthogonal arrays do not have such
a restriction.
Looking into the definition of strong orthogonal arrays, if we focus on t = 2, we can
see that an SOA(N, (s2 )k , 2) of strength two is an orthogonal array of strength one itself
and becomes an orthogonal array of strength two if its s2 levels are collapsed into s levels
according to [a/s]. This implies that an SOA(N, (s2 )k , 2) has the same projection property as
an orthogonal array of strength two. An SOA(N, (s3 )k , 3) however, offers better stratification
and projection property than an s-level orthogonal array of strength three. This is because
an SOA(N, (s3 )k , 3) achieves stratification on s2 × s and s × s2 grids in two-dimensions
and s × s × s grids in three-dimensions while an s-level orthogonal array of strength three
can only promise stratification on s × s grids in two-dimensions and s × s × s grids in
three-dimensions. Similar examinations reveal that to enjoy the benefits of better space-
filling properties, when compared to ordinary orthogonal arrays, strong orthogonal arrays
need to have strength three or higher which may require run sizes that are too large for
experimenters to afford in practice. He et al. (2018) introduced a new class of arrays, called
strong orthogonal arrays of strength two plus.

28
Definition 8.4. An N × k array with entries from {0, 1, . . . , s2 − 1} is called a strong
orthogonal array of size N , k factors, s2 levels, and strength 2+ if any subarray of two
columns can be collapsed into an OA(N, (s2 )1 s1 , 2) and an OA(N, s1 (s2 )1 , 2).

Since their introduction, several research topics on strong orthogonal arrays have been
explored, leading to significant advancements in their theory and applications. First, a num-
ber of studies have focused on developing construction methods for strong orthogonal arrays,
particularly those of strength three and strength two plus. Notable contributions in this area
include He et al. (2018) and Shi and Tang (2020). Second, selecting an optimal strong or-
thogonal array from the broader class based on specific design criteria remains a fundamental
problem. Representative works addressing this challenge include Shi and Tang (2019) and
Chen and Tang (2024). Understanding and characterizing strong orthogonal arrays is an-
other critical research direction. He and Tang (2014) made significant contributions in this
area, providing deeper insights into the properties of strong orthogonal arrays. In addition,
substantial advancements have been achieved by imposing additional structures on strong
orthogonal arrays. For example, Li et al. (2021) and Zhou and Tang (2019) considered
column-orthogonal strong orthogonal arrays; Liu and Liu (2015) studies sliced strong or-
thogonal arrays; Zheng et al. (2024) introduced nested strong orthogonal arrays; Wang et al.
(2022) investigated strong group-orthogonal arrays. Shi et al. (2023) empirically showed the
advantage of strong orthogonal arrays in hyperparameter tuning in deep neural networks,
providing an example of the use of strong orthogonal arrays in machine learning.

8.4 Grouped Orthogonal Arrays

Chen et al. (2025) introduced the concept of grouped orthogonal arrays, motivated by
applications in computer experiments where interactions between factors only arise from
disjoint groups of variables. In such experiments, the true response surface or the preferred
surrogate model is additive, with each component being a function of one specific group of
variables. A practical example is provided in the engine block and head joint sealing exper-
iment discussed by Joseph et al. (2008), where eight factors were involved. Their analysis
identified significant linear, quadratic, and interaction effects among the first, second, and
sixth factors. The eight factors can be divided into two groups: one comprising the first,
second, and sixth factors, and the other containing the remaining factors. In such cases, a
design with superior projection and space-filling properties for the first group is more de-
sirable than traditional space-filling designs that lack this feature because it enables more
accurate estimation of the main and interaction effects within the first group.

Definition 8.5. An orthogonal array is called an s-level grouped orthogonal array with N
runs, g groups and strength t if its factors can be partitioned into g groups where the ith
group has ki factors and is of strength ti , where ti ≥ t, for i = 1, . . . , g.

We use GOA(N, (k1 , k2 , . . . , kg ), (t1 , t2 , . . . , tg ), s, t) to denote such an array. An an ex-

ample, Table 8 displays a GOA(27; (4; 3; 3); (3, 3, 3); 3; 2), for which the whole array is a
three-level orthogonal array of strength two, but where the columns are partitioned into
three groups, denoted by D1 , D2 and D3 in Table 8, that each form an orthogonal array of
strength three.

29
 Table 8: The design matrix D = (D1 , D2 , D3 ) for GOA(27, (4, 3, 3), 3 × 3, 3, 2).

D1 D2 D3
0000 000 000
1110 111 111
2220 222 222
0120 012 012
1200 120 120
2010 201 201
0210 021 021
1020 102 102
2100 210 210
0111 122 200
1221 200 011
2001 011 122
0201 101 212
1011 212 020
2121 020 101
0021 110 221
1101 221 002
2211 002 110
0222 211 100
1002 022 211
2112 100 022
0012 220 112
1122 001 220
2202 112 001
0102 202 121
1212 010 202
2022 121 010

The concept of grouped orthogonal arrays is not new. In addressing the experimental
design issue in physical experiments and applications in combinatorics, Lin (2012) introduced
designs of variable resolution and Raaphorst et al. (2014) coins variable strength orthogonal
arrays. Lin (2012) and Lekivetz and Lin (2016) provided several constructions for designs of
variable resolution but the focus was on two-level designs. The variable strength orthogonal
arrays obtained by Raaphorst et al. (2014) have groups of three factors and their run sizes
are limited to s3 for a prime power s. Zhang et al. (2023) constructed variable strength
orthogonal arrays with strength l containing a subarray with strength greater than l, where
l ≥ 2. In addition to the drawback that the designs constructed have only one group with
larger strength, the resulting designs have very restrictive run sizes st for a prime power s
and an integer t ≥ 4. Chen et al. (2025) proposed several construction methods to generate
many more designs with flexible run sizes and better within-group projection properties for

30
any prime power number of levels.
The research on grouped orthogonal arrays is quite new. Several important directions are
called for. First, current constructions mostly produce regular designs and thus construction
methods for non-regular grouped orthogonal arrays are needed. Second, constructions on
grouped orthogonal arrays with differing group sizes and mixed-level grouped orthogonal
arrays are worth exploring. Another important topic is the use of group orthogonal arrays
in analysis of computer experiments and beyond such as models with blocked additive ker-
nels, as done in Lin and Morrill (2014), which showed the advantages of designs of variable
resolution in model selection of linear models.

9 Summary
This review is based on the experiences and interests of the authors, and therefore in-
complete and selective. Yet, even this selective review demonstrates the enormous impact
that the introduction of orthogonal arrays by Rao (1946, 1947, 1949) has had on applications
and research. Orthogonal arrays are a simple yet elegant mathematical and statistical tool
with a rich theoretical foundation and diverse applications across many fields. This review
highlights fundamental results and recent advancements, hoping to interest more readers in
research on orthogonal arrays and related structures or to use orthogonal arrays in novel
applications. In an era where, in some fields, large datasets are regularly collected, there
continue to be emerging roles for orthogonal arrays, as we have pointed out in this review.
This includes applications in data subsampling and machine learning. There is no doubt
that orthogonal arrays will remain a prominent and versatile tool that will continue to drive
innovation, including in the rapidly evolving realm of artificial intelligence.

Conflicts of Interest Statement

The authors have no conflicts of interest to disclose.

References
Addelman, S. and O. Kempthorne (1962). Orthogonal main-effect plans. Aeronautical Re-
search Laboratory, Office of Aerospace Research US Air Force.

Agrawal, V. and A. Dey (1982). A note on orthogonal main effect plans for asymmetrical
factorials. Sankhyā: The Indian Journal of Statistics, Series B 44, 278–282.

Ai, M., B. Jiang, and K. Li (2014). Construction of sliced space-filling designs based on
balanced sliced orthogonal arrays. Statistica Sinica 24, 1685–1702.

Block, R. M. and R. W. Mee (2005). Resolution IV designs with 128 runs. Journal of Quality
Technology 37 (4), 282–293.

31
Bohyn, A., E. D. Schoen, and P. Goos (2023). Enumeration of regular fractional factorial
designs with four-level and two-level factors. Journal of the Royal Statistical Society Series
C: Applied Statistics 72 (3), 750–769.

Bose, R. (1961). On some connections between the design of experiments and information
theory. Bull. Inst. Inter. Statist. 38, 257–271.

Bose, R. C. and K. A. Bush (1952). Orthogonal arrays of strength two and three. The
Annals of Mathematical Statistics 23 (4), 508–524.

Bose, R. C., S. S. Shrikhande, and E. T. Parker (1960). Further results on the construction of
mutually orthogonal Latin squares and the falsity of Euler’s conjecture. Canadian Journal
of Mathematics 12, 189–203.

Box, G. E. and J. S. Hunter (1961). The 2k−p fractional factorial designs. Technometrics 3 (3),
311–351.

Bush, K. A. (1950). Orthogonal arrays. PhD thesis, North Carolina State University.

Chacko, A. and A. Dey (1981). Some orthogonal main effect plans for asymmetrical factorials.
Sankhyā: The Indian Journal of Statistics, Series B 43, 384–391.

Chacko, A., A. Dey, and G. Ramakrishna (1979). Orthogonal main-effect plans for asym-
metrical factorials. Technometrics 21 (2), 269–270.

Chen, G., Y. He, C. D. Lin, and F. Sun (2025). Grouped orthogonal arrays and their
constructions. Statistica Sinica Accepted.

Chen, G., C. Shi, and B. Tang (2023). Nonregular designs from Paley’s hadamard matrices:
Generalized resolution, projectivity and hidden projection property. Electronic Journal of
Statistics 17 (2), 2120–2138.

Chen, G. and B. Tang (2022). A study of orthogonal array-based designs under a broad
class of space-filling criteria. The Annals of Statistics 50 (5), 2925–2949.

Chen, G. and B. Tang (2024). Selecting strong orthogonal arrays by linear allowable level
permutations. Electronic Journal of Statistics 18 (2), 3573–3589.

Chen, J. (1992). Some results on 2n−k fractional factorial designs and search for minimum
aberration designs. The Annals of Statistics, 2124–2141.

Chen, J. and D. K. Lin (1991). On the identity relationships of 2 − p designs. Journal of

Statistical Planning and Inference 28 (1), 95–98.

Chen, J., D. Sun, and C. Wu (1993). A catalogue of two-level and three-level fractional
factorial designs with small runs. International Statistical Review/Revue Internationale
de Statistique 61, 131–145.

Cheng, C.-S. (1980). Orthogonal arrays with variable numbers of symbols. The Annals of
Statistics 8 (2), 447–453.

32
Cheng, C.-S. (1989). Shorter communication: some orthogonal main-effect plans for asym-
metrical factorials. Technometrics 31 (4), 475–477.

Cheng, C.-S. (1995). Some projection properties of orthogonal arrays. The Annals of Statis-
tics 23 (4), 1223–1233.

Cheng, C.-S. (2016). Theory of factorial design. Chapman and Hall/CRC Boca Raton, FL,
USA.

Cheng, C.-S., D. M. Steinberg, and D. X. Sun (1999). Minimum aberration and model
robustness for two-level fractional factorial designs. Journal of the Royal Statistical Society:
Series B (Statistical Methodology) 61 (1), 85–93.

Cheng, C.-S. and B. Tang (2025). Theory of nonregular factorial designs. Chapman and
Hall/CRC.

Clark, J. B. and A. M. Dean (2001). Equivalence of fractional factorial designs. Statistica

Sinica 11, 537–547.

Colbourn, C. J. and J. H. Dinitz (2001). Mutually orthogonal Latin squares: a brief survey
of constructions. Journal of Statistical Planning and Inference 95 (1-2), 9–48.

Dean, A. M., M. Morris, J. Stufken, and D. Bingham (2015). Handbook of design and analysis
of experiments, Volume 7. CRC Press Boca Raton.

DeCock, D. and J. Stufken (2000). On finding mixed orthogonal arrays of strength 2 with
many 2-level factors. Statistics and Probability Letters 50 (4), 383–388.

Delsarte, P. (1973). An algebraic approach to the association schemes of coding theory.

Philips Res. Rep. Suppl. 10.

Deng, L.-Y. and B. Tang (1999a). Generalized resolution and minimum aberration criteria
for Plackett-Burman and other nonregular factorial designs. Statistica Sinica 9, 1071–1082.

Deng, L.-Y. and B. Tang (1999b). Minimum G2 -aberration for nonregular fractional factorial
designs. the Annals of Statistics 27 (6), 1914–1926.

Dey, A. (2010). Construction of nested orthogonal arrays. Discrete Mathematics 310 (21),
2831–2834.

Dey, A. (2012). On the construction of nested orthogonal arrays. Australas. J Comb. 54,
37–48.

Dey, A. and R. Mukerjee (2009). Fractional factorial plans. John Wiley & Sons.

Dey, A. and G. Ramakrishna (1977). A note on orthogonal main-effect plans. Technomet-

rics 19 (4), 511–512.

Draper, N. R. and T. J. Mitchell (1967). The construction of saturated 2k−p R designs. The
Annals of Mathematical Statistics 38, 1110–1126.

33
Draper, N. R. and T. J. Mitchell (1970). Construction of a set of 512-run designs of resolution
≥ 5 and a set of even 1024-run designs of resolution ≥ 6. The Annals of Mathematical
Statistics 41 (3), 876–887.

Eendebak, P. T., E. D. Schoen, A. R. Vazquez, and P. Goos (2023). Systematic enumeration

of two-level even-odd designs of strength 3. Computational Statistics and Data Analy-
sis 180, 107678.

Fang, K.-T. and R. Mukerjee (2000). Miscellanea. a connection between uniformity and
aberration in regular fractions of two-level factorials. Biometrika 87 (1), 193–198.

Finney, D. J. (1945). The fractional replication of factorial arrangements. Annals of Eugen-

ics 12, 291–331.

Fries, A. and W. G. Hunter (1980). Minimum aberration 2k−p designs. Technometrics 22 (4),
601–608.

Georgiou, S., C. Koukouvinos, and J. Seberry (2003). Hadamard matrices, orthogonal designs
and construction algorithms. DESIGNS 2002: Further Computational and Constructive
Design Theory, 133–205.

Gramacy, R. B. (2020). Surrogates: Gaussian process modeling, design, and optimization

for the applied sciences. Chapman and Hall/CRC.

Hall, M. (1998). Combinatorial theory, Volume 71. John Wiley & Sons.

Hamada, M. and C. J. Wu (1992). Analysis of designed experiments with complex aliasing.

Journal of Quality Technology 24 (3), 130–137.

Hamming, R. W. (1950). Error detecting and error correcting codes. The Bell System
Technical Journal 29 (2), 147–160.

He, X. (2019). Sliced rotated sphere packing designs. Technometrics 61 (1), 66–76.

He, X. and P. Z. Qian (2011). Nested orthogonal array-based Latin hypercube designs.
Biometrika 98 (3), 721–731.

He, Y., C.-S. Cheng, and B. Tang (2018). Strong orthogonal arrays of strength two plus.
The Annals of Statistics 46 (2), 457–468.

He, Y., C. D. Lin, and F. Sun (2022). A new and flexible design construction for orthogonal
arrays for modern applications. The Annals of Statistics 50 (3), 1473–1489.

He, Y. and B. Tang (2013). Strong orthogonal arrays and associated Latin hypercubes for
computer experiments. Biometrika 100 (1), 254–260.

He, Y. and B. Tang (2014). A characterization of strong orthogonal arrays of strength three.
The Annals of Statistics 42 (4), 1347–1360.

34
Hedayat, A., K. Pu, and J. Stufken (1992). On the construction of asymmetrical orthogonal
arrays. The Annals of Statistics 20, 2142–2152.

Hedayat, A., J. Stufken, and G. Su (1996). On difference schemes and orthogonal arrays of
strength t. Journal of Statistical Planning and Inference 56 (2), 307–324.

Hedayat, A. and W. Wallis (1978). Hadamard matrices and their applications. The Annals
of Statistics 6 (6), 1184–1238.

Hedayat, A. S., N. J. A. Sloane, and J. Stufken (2012). Orthogonal arrays: theory and
applications. Springer Science & Business Media.

Hwang, Y., X. He, and P. Z. Qian (2016). Sliced orthogonal array-based Latin hypercube
designs. Technometrics 58 (1), 50–61.

Joseph, V. R., Y. Hung, and A. Sudjianto (2008). Blind kriging: A new method for developing
metamodels. Journal of Mechanical Design, 031102.

Katsaounis, T. I. and A. M. Dean (2008). A survey and evaluation of methods for deter-
mination of combinatorial equivalence of factorial designs. Journal of Statistical Planning
and Inference 138 (1), 245–258.

Katsaounis, T. I., A. M. Dean, and B. Jones (2013). On equivalence of fractional facto-

rial designs based on singular value decomposition. Journal of Statistical Planning and
Inference 143 (11), 1950–1953.

Kline, J. (2019). Geometric search for hadamard matrices. Theoretical Computer Sci-
ence 778, 33–46.

Lander, E. (1983). Symmetric designs: An algebraic approach. Cambridge University Press.

Lekivetz, R. and C. D. Lin (2016). Designs of variable resolution robust to non-negligible

two-factor interactions. Statistica Sinica 26, 1269–1278.

Li, K., B. Jiang, and M. Ai (2015). Sliced space-filling designs with different levels of two-
dimensional uniformity. Journal of Statistical Planning and Inference 157, 90–99.

Li, W., M.-Q. Liu, and J.-F. Yang (2021). Column-orthogonal nearly strong orthogonal
arrays. Journal of Statistical Planning and Inference 215, 184–192.

Lin, C. D. (2012). Designs of variable resolution. Biometrika 99 (3), 748–754.

Lin, C. D. and L. Kang (2016). A general construction for space-filling Latin hypercubes.
Statistica Sinica 26, 675–690.

Lin, C. D. and S. Morrill (2014). Design of variable resolution for model selection. Journal
of Statistical Planning and Inference 155, 127–134.

Lin, C. D., R. Mukerjee, and B. Tang (2009). Construction of orthogonal and nearly orthog-
onal Latin hypercubes. Biometrika 96 (1), 243–247.

35
Lin, C. D. and R. Sitter (2008). An isomorphism check for two-level fractional factorial
designs. Journal of Statistical Planning and Inference 138 (4), 1085–1101.

Lin, D. K. and N. R. Draper (1992). Projection properties of Plackett and Burrnan designs.
Technometrics 34 (4), 423–428.

Liu, H. and M.-Q. Liu (2015). Column-orthogonal strong orthogonal arrays and sliced strong
orthogonal arrays. Statistica Sinica 25, 1713–1734.

Ma, C.-X., K.-T. Fang, and D. K. Lin (2001). On the isomorphism of fractional factorial
designs. Journal of Complexity 17 (1), 86–97.

MacWilliams, F. and N. J. A. Sloane (1977). The theory of error-correcting codes. Elsevier

Science Publishers.

McKay, M. D., R. J. Beckman, and W. J. Conover (1979). A comparison of three methods

for selecting values of input variables in the analysis of output from a computer code.
Technometrics 42 (1), 55–61.

Mukerjee, R., P. Z. Qian, and C. J. Wu (2008). On the existence of nested orthogonal arrays.
Discrete Mathematics 308 (20), 4635–4642.

Nordstrom, A. W. and J. P. Robinson (1967). An optimum nonlinear code. Information and

Control 11 (5-6), 613–616.

Owen, A. B. (1992). Orthogonal arrays for computer experiments, integration and visual-
ization. Statistica Sinica 2, 439–452.

Paley, R. E. (1933). On orthogonal matrices. Journal of Mathematics and Physics 12 (1-4),

311–320.

Pang, F., M.-Q. Liu, and D. K. Lin (2009). A construction method for orthogonal Latin
hypercube designs with prime power levels. Statistica Sinica, 1721–1728.

Pang, S., J. Wang, D. K. Lin, and M.-Q. Liu (2021). Construction of mixed orthogonal
arrays with high strength. The Annals of Statistics 49 (5), 2870–2884.

Pang, S. and Y. Zhu (2024). On the construction of nested orthogonal arrays with the
adjacent numbers of levels. Stat 13 (2), 1–15.

Plackett, R. L. and J. P. Burman (1946). The design of optimum multifactorial experiments.

Biometrika 33 (4), 305–325.

Qian, P. Z., M. Ai, and C. J. Wu (2009). Construction of nested space-filling designs. The
Annals of Statistics 37 (3), 3616–3643.

Qian, P. Z., B. Tang, and C. J. Wu (2009). Nested space-filling designs for computer exper-
iments with two levels of accuracy. Statistica Sinica 19, 287–300.

36
Qian, P. Z. and C. J. Wu (2008). Bayesian hierarchical modeling for integrating low-accuracy
and high-accuracy experiments. Technometrics 50 (2), 192–204.

Qian, P. Z. and C. J. Wu (2009). Sliced space-filling designs. Biometrika 96 (4), 945–956.

Qian, P. Z. G., H. Wu, and C. J. Wu (2008). Gaussian process models for computer experi-
ments with qualitative and quantitative factors. Technometrics 50 (3), 383–396.

Raaphorst, S., L. Moura, and B. Stevens (2014). A construction for strength-3 covering arrays
from linear feedback shift register sequences. Designs, Codes and Cryptography 73 (3),
949–968.

Raghavarao, D. (1988). Constructions and combinatorial problems in design of experiments.

Dover Publications.

Rao, C. R. (1946). Hypercubes of strength’d’leading to confounded designs in factorial

experiments. Bull. Calcutta Math. Soc. 38, 67–78.

Rao, C. R. (1947). Factorial experiments derivable from combinatorial arrangements of

arrays. Supplement to the Journal of the Royal Statistical Society 9 (1), 128–139.

Rao, C. R. (1949). On a class of arrangements. Proceedings of the Edinburgh Mathematical

Society 8 (3), 119–125.

Rao, C. R. (1973). Some combinatorial problems of arrays and applications to design of

experiments. In J. N. Srivastava (Ed.), A Survey of Combinatorial Theory, Chapter 29,
pp. 349–359. Amsterdam: North Holland.

Santner, T. J., B. J. Williams, W. I. Notz, and B. J. Williams (2003). The design and
analysis of computer experiments, Volume 1. Springer.

SAS (2019). A library of difference schemes. [Link]

technote/ts723_DifferenceSchemes.txt.

Schmidt, R. R., E. E. Cruz, and M. Iyengar (2005). Challenges of data center thermal
management. IBM Journal of Research and Development 49 (4.5), 709–723.

Schoen, E. D., P. T. Eendebak, and M. V. Nguyen (2010). Complete enumeration of pure-

level and mixed-level orthogonal arrays. Journal of Combinatorial Designs 18 (2), 123–140.

Seiden, E. (1954). On the problem of construction of orthogonal arrays. The Annals of

Mathematical Statistics 25 (1), 151–156.

Shi, C., A. K. Chiu, and H. Xu (2023). Evaluating designs for hyperparameter tuning in deep
neural networks. The New England Journal of Statistics in Data Science 1 (3), 334–341.

Shi, C. and B. Tang (2018). Designs from good hadamard matrices. Bernoulli 24 (1),
661–671.

37
Shi, C. and B. Tang (2019). Design selection for strong orthogonal arrays. Canadian Journal
of Statistics 47 (2), 302–314.

Shi, C. and B. Tang (2020). Construction results for strong orthogonal arrays of strength
three. Bernoulli 26 (1), 418–431.

Sitter, R. R. (1993). Balanced repeated replications based on orthogonal multi-arrays.

Biometrika 80 (1), 211–221.

Sloane, N. J. A. and J. Stufken (1996). A linear programming bound for orthogonal arrays
with mixed levels. Journal of Statistical Planning and Inference 56 (2), 295–305.

Stein, M. (1987). Large sample properties of simulations using latin hypercube sampling.
Technometrics 29 (2), 143–151.

Steinberg, D. M. and D. K. Lin (2006). A construction method for orthogonal Latin hyper-
cube designs. Biometrika 93 (2), 279–288.

Stinson, D. R. (2006). A tutorial on orthogonal arrays: Constructions, bounds and links to

error-correcting codes.

Stufken, J. and B. Tang (2007). Complete enumeration of two-level orthogonal arrays of

strength d with d + 2 constraints. The Annals of Statistics 35, 793—-814.

Suen, C.-Y., A. Das, and A. Dey (2001). On the construction of asymmetric orthogonal
arrays. Statistica Sinica 11, 241–260.

Sun, D. X., W. Li, and K. Q. Ye (2002). An algorithm for sequentially constructing noniso-
morphic orthogonal designs and its applications. preprint.

Sun, F., M.-Q. Liu, and P. Z. Qian (2014). On the construction of nested space-filling designs.
The Annals of Statistics 42 (4), 1394–1425.

Sun, F. and B. Tang (2017a). A general rotation method for orthogonal Latin hypercubes.
Biometrika 104 (2), 465–472.

Sun, F. and B. Tang (2017b). A method of constructing space-filling orthogonal designs.

Journal of the American Statistical Association 112 (518), 683–689.

Sylvester, J. J. (1867). Lx. thoughts on inverse orthogonal matrices, simultaneous signsucces-

sions, and tessellated pavements in two or more colours, with applications to newton’s rule,
ornamental tile-work, and the theory of numbers. The London, Edinburgh, and Dublin
Philosophical Magazine and Journal of Science 34 (232), 461–475.

Tang, B. (1993). Orthogonal array-based Latin hypercubes. Journal of the American Sta-
tistical Association 88 (424), 1392–1397.

Tsai, S.-F. (2022). Some results on the construction of sliced orthogonal arrays of parallel-
flats type. Communications in Statistics-Theory and Methods 51 (3), 569–580.

38
Vazquez, A. R. and H. Xu (2019). Construction of two-level nonregular designs of strength
three with large run sizes. Technometrics 61 (3), 341–353.

Wang, C., J. Yang, and M.-Q. Liu (2022). Construction of strong group-orthogonal arrays.
Statistica Sinica 32 (3), 1225–1243.

Wang, H., M. Yang, and J. Stufken (2019). Information-based optimal subdata selection
for big data linear regression. Journal of the American Statistical Association 114 (525),
393–405.

Wang, J. and C. Wu (1991). An approach to the construction of asymmetrical orthogonal

arrays. Journal of the American Statistical Association 86 (414), 450–456.

Wang, K. and Y. Li (2013). Constructions of nested orthogonal arrays. Journal of Combi-

natorial Designs 21 (10), 464–477.

Wang, K. and J. Yin (2013). Further results on the existence of nested orthogonal arrays.
Designs, Codes and Cryptography 67, 233–243.

Wang, L., J. Elmstedt, W. K. Wong, and H. Xu (2021). Orthogonal subsampling for big
data linear regression. The Annals of Applied Statistics 15 (3), 1273–1290.

Wang, P. (1990). On the constructions of some orthogonal main-effect plans. Sankhyā: The
Indian Journal of Statistics, Series B 52, 319–323.

Weng, L.-C., K.-T. Fang, and A. Elsawah (2023). Degree of isomorphism: a novel criterion
for identifying and classifying orthogonal designs. Statistical Papers 64 (1), 93–116.

Williamson, J. (1944). Hadamard’s determinant theorem and the sum of four squares. Duke
Math. J. 11, 65–81.

Wilson, R. M. (1974). Concerning the number of mutually orthogonal Latin squares. Discrete
Mathematics 9 (2), 181–198.

Wu, C. (1991). Balanced repeated replications based on mixed orthogonal arrays.

Biometrika 78 (1), 181–188.

Wu, C. J. and M. S. Hamada (2011). Experiments: planning, analysis, and optimization.

John Wiley & Sons.

Xu, H. (2003). Minimum moment aberration for nonregular designs and supersaturated
designs. Statistica Sinica 13, 691–708.

Xu, H. (2005). A catalogue of three-level regular fractional factorial designs. Metrika 62 (2),
259–281.

Xu, H., F. K. Phoa, and W. K. Wong (2009). Recent developments in nonregular fractional
factorial designs. Statistics Surveys 3, 18–46.

39
Xu, X., B. Haaland, and P. Z. Qian (2011). Sudoku-based space-filling designs.
Biometrika 98 (3), 711–720.

Zhang, Q., S. Pang, and Y. Li (2023). On the construction of variable strength orthogonal
arrays. IEICE Trans. Foundamentals E106, 683–688.

Zhang, T., J. Jiang, and Z. Li (2018). Some families of asymmetric nested orthogonal arrays
and asymmetric sliced orthogonal arrays. Statistics and Probability Letters 141, 109–113.

Zhang, T.-f., G. Wu, and A. Dey (2019). Construction of some new families of nested
orthogonal arrays. Communications in Statistics-Theory and Methods 48 (3), 774–779.

Zhang, Y., L. Wang, X. Zhang, and H. Wang (2024). Independence-encouraging subsampling

for nonparametric additive models. Journal of Computational and Graphical Statistics 33,
1–10.

Zheng, C., W. Li, and J.-F. Yang (2024). Nested strong orthogonal arrays. Statistical
Papers 65, 1–26.

Zhou, Y. and B. Tang (2019). Column-orthogonal strong orthogonal arrays of strength two
plus and three minus. Biometrika 106 (4), 997–1004.

Zhu, J., L. Wang, and F. Sun (2024). Group-orthogonal subsampling for hierarchical data
based on linear mixed models. Journal of Computational and Graphical Statistics 33,
1–10.

40