Ascent with Quadratic Assistance for the Construction of Exact Experimental Designs
Lenka Filová Radoslav Harman
stat.CO 62K05

	Scalable Mutual Information Estimation using Dependence Graphs
    Christine Nattrass
    December 30, 2023
==============================================================


In the area of statistical planning, there is a large body of theoretical knowledge and computational experience concerning so-called optimal approximate designs of experiments. However, for an approximate design to be executed in practice, it must be converted into an exact, i.e., integer, design, which is usually done via rounding procedures. Although rapid, rounding procedures have many drawbacks; in particular, they often yield worse exact designs than heuristics that do not require approximate designs at all.
		
In this paper, we build on an alternative principle of utilizing optimal approximate designs for the computation of optimal, or nearly-optimal, exact designs. The principle, which we call ascent with quadratic assistance (AQuA), is an integer programming method based on the quadratic approximation of the design criterion in the neighborhood of the optimal approximate information matrix.
		
To this end, we present quadratic approximations of all Kiefer's criteria with an integer parameter, including D- and A-optimality and, by a model transformation, I-optimality. Importantly, we prove a low-rank property of the associated quadratic forms, which enables us to apply AQuA to large design spaces, for example via mixed integer conic quadratic solvers. We numerically demonstrate the robustness and superior performance of the proposed method for models under various types of constraints. More precisely, we compute optimal size-constrained exact designs for the model of spring-balance weighing, and optimal symmetric marginally restricted exact designs for the Scheffè mixture model. We also show how can iterative application of AQuA be used for a stratified information-based subsampling of large datasets under a lower bound on the quality and an upper bound on the cost of the subsample. 
		
Keywords: optimal design, rounding procedures, mixed integer conic quadratic programming, experimental constraints, subsampling
		






§ INTRODUCTION


Suppose that we intend to perform an experiment consisting of a set[Note that we implicitly assume that reordering of the trials does not influence the relevant properties of the experimental design.] of trials. Assume that the observed response in each trial depends on a design point chosen from a finite design space ={_1,…,_n}. For instance,  may be the set of all available combinations of levels of several discrete factors[In some experimental situations, the set of available design points can be modeled as a continuous domain. However, in many applications, the design space is finite. This is the case if each factor has - in principle or effectively - only a finite number of levels that the experimenter can select, or if the optimal design problem corresponds to data sub-selection (see the examples in Section <ref>). Moreover, the method proposed in this paper can also be useful for solving the problems with continuous design spaces; cf. Section <ref>.]. 


An “exact” design (ED) is a selection ξ of design points, not necessarily distinct, to be used for individual trials. We will formalize an ED ξ as a non-negative integer-valued vector (ξ_1,…,ξ_n)^T ∈ℕ_0^n,[The symbols , _+, , _0, and ^k × n denote the sets of real, non-negative real, natural, non-negative integer numbers, and the set of all k × n real matrices, respectively.] where ξ_i, called the i-th weight, represents the number of trials to be performed at the design point _i, i=1,…,n.[Therefore, we do not represent designs by normalized (probability) measures, as is frequently done in optimal design, but by non-normalized vectors of numbers of trials.] An “approximate” design (AD), ξ=(ξ_1,…,ξ_n)^T ∈_+^n, is allowed to have general non-negative components, which means that the weight ξ_i is a continuous relaxation of the integer number of trials to be performed at _i, i=1,…,n.[Approximate designs are sometimes also called “continuous” designs, which refers to the continuity of the space of designs, not the design space.] Thus, an AD must be converted into an ED prior to its application in a real experiment. 


Let Ξ^E_,={ξ∈_0^n: ξ≤}̱ be a non-empty set of permissible EDs, where ∈^k × n and ∈̱^k. In the classical situation[The symbols _n, _n, _n and _n denote the n-dimensional vector of ones, n-dimensional vector of zeros, the n × n unit matrix and the n × n matrix of ones, respectively.] =_n^T and =̱N ∈; in that case we only restrict the number N of trials, the so-called size of the experiment. Nevertheless, there are also many situations where  and $̱ are more complex. They can correspond to various time, budget, material, unbiasedness and safety restrictions, or requirements on the form of the design (see, e.g., <cit.> and Section <ref> of this paper). An important constraint necessary for applications to subsampling is the that each design point (i.e., data-point) can be used only once; formally=_nand=̱_n.[In actual computation using integer programming solvers this “without replication” constraint can be forced by setting the type of variables to binary.]

Suppose that the information gained from an experiment based onξ∈Ξ^E_,can be represented by a matrix(ξ). For instance,(ξ)may be proportional to the Fisher information matrix for the unknown parameters of an underlying statistical model. In optimal experimental design, it is usual to select a concave functionΦwith a target set∪{-∞}to quantify the information content of(ξ). Such an optimality criterion allows an experimenter to compare different designs and, in principle, to select aΦ-optimal ED, i.e., a design that maximizes[Alternatively, it is possible to select a convex criterion Φ such that Φ((ξ)) can be interpreted as a loss from the experiment that depends on the design ξ. In this case, the optimal design would minimize Φ((·)) over Ξ^E_,. Note that some criteria do not depend on the design via its information matrix; we will not discuss them in this paper.]Φ((·))over the discrete setΞ^E_,. In many cases, the information matrix can be consistently extended to ADs. Then, maximizingΦ((·))over the convex setΞ^A_,={ ξ∈_+^n: ξ≤}̱results in the so-called optimal AD.


The construction of optimal EDs is typically a difficult problem of discrete optimization. There are two general approaches to computing optimal or nearly-optimal EDs (see, e.g., <cit.> for a survey):

  *  Convex computational methods, in which an optimal AD is first determined and a process called “rounding” is then used to obtain an ED;
  
  *  Computational methods of discrete optimization, including complete or partial enumeration methods, as well as various specialized or general-purpose solvers and heuristics of mathematical programming.
 

It is usually much simpler to determine an optimal AD than an optimal ED, both theoretically and computationally. Therefore, a large part of the literature is concerned only with approximate designs (cf., <cit.>, <cit.>). Although ADs cannot be directly used for conducting experiments, relatively little attention has been paid to their conversion into efficient EDs. 


The standard methods for converting an AD into an ED are called rounding algorithms, developed for the classical, size-constrained problem. A rounding algorithm begins with an ADξ^*=(ξ_1^*,…,ξ_n^*)^T ∈Ξ^A__n^T,Nand extracts a vector𝒲=(ξ^*_i_1, …, ξ^*_i_s)^Tof positive weights, where{i_1,…,i_s}={i: ξ^*_i>0}is the support ofξ^*andsis the size of the support. Then, typically using simple rules, the algorithm converts𝒲into a vector(ξ^±_i_1,…,ξ^±_i_s)^T ⊂_0^ssuch that∑_j ξ^±_i_j=Nand, finally, transforms the vector of rounded weights into an ED that belongs toΞ^E__n^T,N.


The first notable rounding method was suggested by <cit.>, who formulated the rounding problem as the minimization of the maximum of the difference between the exact and approximate design weights. By using techniques similar to those applied in voting apportionment, <cit.> arrived at a criterion-independent rounding algorithm known as efficient rounding (ER). More recent proposals include randomized rounding heuristics, e.g., proportional and pipage rounding, as well as incremental rounding, and bounds on the approximation ratios of the resulting designs have been presented (see <cit.> and <cit.>). However, these methods are only applicable if the criterion function is submodular (e.g.,D-optimality). We are not aware of rounding procedures for generalΞ^E_,, but for specific classes of constraints, it is not difficult to mimic the existing rounding procedures originally developed for the size-constrained problem. 


ER and its variants, although prevalent to this day, have several major drawbacks. First, for any positive coordinate of the initial AD, the value of the corresponding coordinate of the resulting ED is forced to be at least1. This implies the restrictionN ≥s, which can completely prevent the application of ER if the support size of the AD is large. From the opposite perspective, if a design point is not present in the support of the AD, then ER cannot add a corresponding design point into the resulting ED. Moreover, ER does not account for any design criterion nor any underlying statistical model, although it is based on an optimal AD, which itself can strongly depend on the adopted criterion and model. In addition, for many statistical models, an infinite number of optimal ADs exist, and it is unclear which of them should be used for the rounding operation. All of these disadvantages generally make approach (ii) preferable to (i) in practice; see, e.g., the examples in <cit.>. 
<cit.> proposed a substantially different approach to the use of an optimal AD for ED construction, which overcomes many disadvantages of ER and similar methods. In particular, it does not depend on the choice of the optimal AD if the AD is not unique, it is not restricted to the support of the optimal AD, and the resulting EDs are usually significantly more efficient than the EDs computed by ER. The method is based on a second-order approximation of theD-criterion in the neighborhood of theD-optimal approximate information matrix, and to arrive at an ED, it employs rapid off-the-shelf solvers for integer quadratic programming (IQP).


In this paper, we view the idea of a quadratic criterion approximation based on an optimal or nearly-optimal AD as a broadly applicable principle in computational experimental design, and we call this principle AQuA (ascent with quadratic assistance). As we will show, AQuA can be realized by means of heuristics but also via solvers of IQP or mixed integer conic quadratic programming (MICQP) in situations with various budget and structural constraints on the design. AQuA can also be used sequentially, similarly to the sequential quadratic programming.


The new results of this paper demonstrate that AQuA can be applied to a wide range of criteria, including the important criteria ofA- andI-optimality, and, utilizing a low-rank property of key quadratic forms, to much larger design spaces than competing methods.


This paper is organized as follows: In Section <ref>, we present the general statistical model that we consider and two versions of Kiefer'sΦ_p-criteria. Subsequently, in Section <ref>, we demonstrate how to compute quadratic approximations of these criteria. We propose a low-rank method for the efficient application of AQuA in Section <ref>. This leads to the main result of this paper, a MICQP formulation of AQuA that can be practically used for large structured or unstructured design spaces. Section <ref> provides various remarks. Finally, Section <ref> presents examples of optimal designs that can be computed by the application of the AQuA approach.





§ THE MODEL AND KIEFER'S CRITERIA

 
For a trial in_i ∈,i=1,…,n, the observed responseY(_i)is anr-dimensional random vector that is assumed to satisfy the linear regression modelE(Y(_i))=^T_i β, whereβ∈^mis a vector of unknown parameters and_i ∈^m ×ris a known matrix. For different observations, the errors are assumed to be independent and identically distributed with a finite and non-zero variance. Note that we consider a linear regression model with homoscedastic errors only for simplicity. It is straightforward to use the results of this paper for the construction of locally optimal designs of non-linear regression models (it only requires a linearization in a nominal parameter of the model; see, e.g., <cit.>, Chap. 17) and also to heteroscedastic observations (by means of a proper transformation of the model; see <cit.>, Chap. 23).


The information matrix associated with a designξon, either exact or approximate, is

    (ξ)=∑_i=1^nξ_i _̋i,

where the_̋i=_i^T_i,i=1,…,n, are non-negative definite “elementary” information matrices with dimensions ofm×m[For brevity, we will henceforth use 𝒮^m, 𝒮^m_+, and 𝒮^m_++  to denote the sets of all symmetric, non-negative definite and positive definite m × m matrices, respectively.]. For the classical case with univariate observations,_i=(_i) ∈^m, i.e., them-dimensional regressor corresponding to_i ∈. The general form of the elementary information matrices may be useful for instance for problems with grouped or multivariate observations with possibly correlated components (see <cit.>, Sec. II.5.3.), optimal augmentation of a set of existing trials (as shown in <cit.>, Section 6) and elsewhere.


LetΦ: 𝒮^m_+ →∪{-∞}be a continuous optimality criterion that attains its smallest value for singular non-negative definite matrices. Note that aΦ-optimal ED exists becauseΞ^E_,is finite and non-empty. SinceΞ^A_,is a non-empty compact set, aΦ-optimal AD is also guaranteed to exist. Ifξ^*is aΦ-optimal exact (approximate) design, then(ξ^*)is called aΦ-optimal exact (approximate) information matrix. To make the optimal design problem non-trivial, we will suppose that there exists aξ∈Ξ^E_,such that(ξ)is non-singular, which implies that both the approximate and exactΦ-optimal information matrices are non-singular. We will also assume thatΦis twice differentiable in𝒮^m_++and that there exists a version[By two versions of a criterion, we mean two criteria that induce the same ordering on the set of information matrices.] ofΦthat is strictly concave in the optimal approximate information matrix_*. This assumption is satisfied for most models and standard optimality criteria, and it implies that theΦ-optimal approximate information matrix is unique. 


All properties stated above are satisfied by Kiefer's criteria, which are commonly used in practice. In the optimal design literature, several versions of Kiefer's criteria forΦ_p-optimality appear, and usually the choice of the particular version does not affect the strength of the theoretical or computational results. However, it turns out that in general, criterion-approximation methods do depend on the particular version of the criterion that is chosen. Therefore, we will consider two concave versions ofΦ_pcriteria, as follows.


The “positive” version (cf. <cit.>): Forp ∈and∈𝒮^m_++,
 let 

    Φ^+_p()=(1/m (^-p) )^-1/p

andΦ^+_p()=0for a singular matrix∈𝒮^m_+. In particular, forp=1, we obtain the criterionΦ^+_1ofA-optimality. The corresponding criterion ofD-optimality is defined asΦ^+_0()=(())^1/mfor all∈𝒮^m_+.


The “negative” version (cf. <cit.>, Section IV.2.7): Forp ∈and∈𝒮^m_++, let  

    Φ^-_p()=-(1/m(^-p) )^1/p

andΦ^-_p()= -∞for a singular matrix∈𝒮^m_+. In particular,Φ^-_1is a version of theAcriterion, and the correspondingDcriterion isΦ^-_0()=-(())^-1/mfor∈𝒮^m_++orΦ^-_0()= -∞for a singular∈𝒮^m_+. 


Another commonly used concave version of theD-optimality criterion isΦ^0_0()=log(())for all∈𝒮^m_+(log(0):=-∞); cf. <cit.>. This is the version ofD-optimality used in <cit.>.


Note that both positive and negative versions of the criterion are smooth on the set of positive definite matrices, and the gradients are

    ∇_Φ^±_p()=±Φ^±_p()/ (^-p)^-p-1.


It is customary to evaluate the quality of a design with respect to the optimal AD. Letξ^*be the optimal AD, and letΦbe a non-negative, positively homogeneous criterion that is not constantly equal to zero (these conditions are satisfied by the positive versionΦ^+_pof Kiefer's criteria). Then, theΦ-efficiency of a designξis defined aseff_Φ(ξ)=Φ((ξ))/Φ((ξ^*)), see <cit.>, Section 5.15. 





§ QUADRATIC APPROXIMATIONS OF KIEFER'S CRITERIA


Suppose that we have a quadratic approximationΦ_Q: 𝒮^m_+ →of a concave criterionΦin the neighborhood of_*. Our experience shows that in most optimal design problems, the ordering onΞ^E_,that is induced byΦ_Qlargely coincides with the ordering induced by the original criterionΦ. At the same time, the quadratic approximation criterionΦ_Qcan be evaluated (or updated) much more rapidly thanΦ, it has a simpler analytic properties and there are powerful available solvers that can maximizeΦ_Q.


Let_*∈𝒮^m_++denote theΦ-optimal approximate information matrix[Note that the optimal approximate information matrix _* with respect to Φ_p^+ and Φ_p^- is non-singular for any p ∈_0.] and letΦbe twice differentiable in𝒮^m_++. Then, a second-order Taylor approximation ofΦin terms of∈𝒮^m_+can be written as follows
(see, e.g., <cit.>, Appendix D) 

    Φ()    ≈   Φ(_*)+∂Φ(_*,-_*) 
           +1/2∂^2Φ(_*,-_*),

where∂Φ(,)denotes the directional derivative at the pointin the direction, i.e.,∂Φ(,)=(∇_Φ()), and∂^2Φ(,)denotes the second directional derivative at the pointin the direction, i.e.,∂^2Φ(,)=(∇_∂Φ(,)), with∇_Φ()denoting the gradient with respect to. 


ForΦ=Φ_p^+as defined in (<ref>) withp ∈_0, we obtain

    ∂Φ_p^+(_*,-_*)=Φ_p^+(_*)/(_*^-p)(_*^-p-1(-_*))

and

    ∂^2Φ_p^+(_*,-_*)=Φ_p^+(_*)/(_*^-p)[(p+1)^2(_*^-p-1)/(_*^-p).
    . -ℱ_p(_*,, )],

where 

    ℱ_p(_*,_1,_2)=∑_r=1^p+1(_*^-r_1_*^-p-2+r_2).

Note that in particular,

    ℱ_0(_*,, )   =   ([_*^-1]^2) 
    ℱ_1(_*,, )   =   2(_*^-2_*^-1).

According to (<ref>), we haveΦ_p^+()≈Φ_pQ^+(), whereΦ_pQ^+is the second-order approximation of the criterion ofΦ_p^+-optimality and is given by

    Φ_pQ^+()=Φ_p^+(_*)/(_*^-p)( (_*^-p-1)+ .
    .p+1/2^2(_*^-p-1)/(_*^-p)-1/2ℱ_p(_*,,)).


Similar computations can be performed forΦ=Φ_p^-,p ∈_0, as defined in (<ref>), leading to the quadratic approximation 

    Φ_pQ^-()=-3Φ_p^-(_*)/(_*^-p)((_*^-p-1)-(_*^-p).
    . +p-1/6^2(_*^-p-1)/(_*^-p)
    	-1/6ℱ_p(_*,, )).


We remark that it is also possible to computeΦ_pQ^-based on the formulas for the Hessian of a modified version ofΦ_p^-that were derived by <cit.> (for integer values ofp) and <cit.> (forp=0,1).


Because the mappingξ→(ξ)is linear,Φ_Q((·))is a quadratic function on^n_+; i.e.,Φ_Q((·))=ϕ_,(·)+cfor some∈^n,∈𝒮^n_+andc ∈, where

    ϕ_,(ξ)=^Tξ - ξ^Tξ,   ξ∈^n_+.


Then, the problem of optimal ED based on the AQuA approach can be expressed as the integer quadratic problem

    .[      max_ξ    ϕ_,(ξ),; subject to   ξ∈Ξ^E_,. ].


For a general criterion, there are several possible ways of constructing the appropriate vectorand matrix, for instance, through the use of standard numerical differentiation techniques. However, for Kiefer's criteria withp ∈_0, it is simple to derive analytical forms forand, as we show next.


Consider an EDξ=(ξ_1,…,ξ_n)^Twith the information matrix=(ξ). Clearly,(_*^-p-1)=∑_i ξ_i (_*^-p-1_̋i), and forr=1,…, p+1,

    ℱ_p(_*,,)=∑_i,j=1^nξ_i ξ_j ℱ_p(_*,_̋i,_̋j);

therefore, the maximization ofΦ_pQ^+((·))overΞ^E_,is equivalent to the integer quadratic optimization problem expressed in (<ref>), where=_p^+has the components(_p^+)_i=(_*^-p-1_̋i),i=1,…,n, and the matrix=_p^+has the elements[Note that the matrix _p^+ is symmetric, as is the matrix _p^- defined below, because (_1_̋1_2_̋2)=(_1_̋2_2_̋1) for the symmetric non-negative definite matrices _1, _2, _̋1, and _̋2.]
    (_p^+)_i,j=p+1/2(_p^+)_i(_p^+)_j/(_*^-p)-1/2ℱ_p(_*,_̋i,_̋j),i,j=1,…,n. Similarly, the maximization ofΦ_pQ^-((·))overΞ^E_,is equivalent to the integer quadratic optimization problem expressed in (<ref>), where=_p^-has the components(_p^-)_i=(_*^-p-1_̋i),i=1,…,n, and the matrix=_p^-has the elements

    (_p^-)_i,j=p-1/6(_p^-)_i(_p^-)_j/(_*^-p)-1/6ℱ_p(_*,_̋i,_̋j),i,j=1,…,n.





§ EFFICIENT COMPUTATIONAL APPROACH TO AQUA


 §.§ A low-rank property of the quadratic approximations of Kiefer's criteria


We can use quadratic approximation of criteria in combination with many algorithms for optimal ED (e.g., <cit.>, <cit.>, <cit.>). To do so, we must be able to compute the values of the quadratic functionϕ_,given in (<ref>) for designsξ, as required by the algorithm. We will show that for the quadratic approximation criteria resulting from the optimal ED problem based on the Kiefer's criteria, this computation can be performed rapidly, based on a low-rank property of the associated quadratic forms. As a key by-product, we will obtain a useful quadratic cone representation of the AQuA optimization problem.


The ability to efficiently numerically evaluate multivariate quadratic functions of the formϕ_,generally depends on various specifics of the problem at hand, the known theoretical properties ofand, the selected optimization algorithm, and the available hardware. Here, we will consider problems that are typical of optimal experimental design. In particular, we will assume thatmis a small number (usually less than10), whereasnis a much larger number, possibly ranging from the order of tens to hundreds of thousands.


Let the functionϕ_,be based on the quadratic approximation of a criterion defined on the set of information matrices. That is,ϕ_,(·)=Φ_Q((·)), whereΦ_Q()is a quadratic function of the elements of. For a designξ, the most problematic part of computingϕ_,(ξ)is the evaluation of the quadratic formξ^Tξfor then ×nmatrix, becausenis often large. However, as we will show, we can construct a matrixwith dimensions ofn ×t, wheret ≤s:=m(m+1)/2 ≪n, such that=^T. Importantly, we can constructwithout computing; i.e., we can completely avoid working with potentially enormous matrices.


To this end, let the functionΦ_Q: 𝒮^m_+ →be represented in the form   

    Φ_Q()=a(^T() - (())^T ()) + c,

where ∈^s, ∈𝒮^s_+, anda>0,care real numbers which do not influence the maximum. Let_m∈ℝ^m^2×sbe the duplication matrix that relates theandvecoperators[The symbols  and vec denote the vectorization and half-vectorization of a matrix, respectively.]; i.e.,vec()=_m (). Then, the versions of Kiefer's criteria defined in the previous section can be represented using Theorem 16.2.2. from <cit.> and the formulas

    ()   =   (vec())^T_m (), 
    ^2()   =   ()^T _m^T vec()
          (vec())^T _m (),
    (_1_2)   =   ()^T ^T_m (_2⊗_1)
          _m (),

which are valid for all, , _1, _2 ∈𝒮^m; thus, we obtain

    Φ_pQ^±()=a^±((_p^±)^T() - ()^T_p^± ()) + c^±,

where 

    a^+=Φ_p^+(_*)/(_*^-p),    a^-=-3Φ_p^-(_*)/(_*^-p),

    ^+_p=^-_p=_m^Tvec(_*^-p-1),

    ^+_p   =   _m^T[-1+p/2vec(_*^-p-1)(vec(_*^-p-1))^T/(_*^-p)+.
          .1/2∑_r=1^p+1_*^-p-2+r⊗_*^-r]_m,
    ^-_p   =   _m^T[1-p/6vec(_*^-p-1)(vec(_*^-p-1))^T/(_*^-p)+.
          .1/6∑_r=1^p+1_*^-p-2+r⊗_*^-r]_m,

andc^+=0,c^-= 3Φ_p^-(_*).


Next, we can construct a decomposition=^Tsuch that thes ×tmatrixis of rankt,[Note that  can be a singular non-negative definite matrix; therefore, t can be even smaller than s.] using, for instance, the Cholesky algorithm or the singular value decomposition. We have 

    ϕ_,(ξ)   =   Φ_Q(∑_i=1^nξ_i_̋i)=∑_i=1^nξ_i ^T(_̋i) 
          - ∑_i=1^n∑_j=1^nξ_iξ_j((_̋i))^T (_̋j)
       =   ξ^T (̃̋)^T - ∑_i=1^n∑_j=1^nξ_iξ_j(̃̋^T^̋T)_i,j
       =   ^Tξ - ^T ξ^2,

where=̋((_̋1),…,(_̋n))^Tis ann ×smatrix,=̃̋, and=̃̋. Equation (<ref>) allows us to computeϕ_,(ξ)without evaluating and storing.


An advantage of the previous expression is that with the use of,ϕ_,(ξ)can be rapidly evaluated; for instance, the exchange step in an exchange algorithm (<cit.>, Sec. 12.3) can be performed based on the equation

    ϕ_,(ξ+_l-_k)=
          ϕ_,(ξ)+_l-_k-2(^Tξ)^T[_l·- _k·]-
          _l·^2+2(_l·)^T_k·-_k·^2,

where_l,_kare thel-th andk-th standard unit vectors, and_l·,_k·are thelth and thekth rows of. Note that^Tξis updated as follows:^T(ξ+_l-_k)=^Tξ+_l·-_k·. If thenvalues of_k·^2,k=1,…,n, are precomputed and stored in memory, then each update involves only2t+2multiplications and4t+4subtractions or additions. An example of how these formulas can be utilized with heuristic exchange algorithm can be found in a preprint of the previous version of this paper; see <cit.>. Here we will focus on a more versatile application of the low-rank property, as detailed in the next section.




 §.§ Mixed integer conic quadratic programming formulation of AQuA


Once we proved the low-rank property=^T, whereis ann ×tmatrix,t ≪n, we can use a known trick to reformulate the problem of quadratic programming (e.g., <cit.>, Chapter 10). Introducing an auxiliary continuous variabler, the optimization problem <ref> of the AQuA approach can be written as

    .[         max_ξ, r          ^Tξ - r;             s.t. ξ≤,̱ ξ≥_n, ξ∈ℤ^n;                    r ≥ ||^Tξ||^2. ].

	
It is simple to verify that the last constraint in (<ref>) can be expressed as(1/2,r,ξ^T)^T ∈Q^2+t, whereQ^2+tis the second-order cone

    Q^2+t={(a,b,^̌T)^T: a ≥(b,^̌T)^T}

andis the orthogonal matrix

    =
      ( [  1/√(2)  1/√(2)    _t^T;  1/√(2) -1/√(2)    _t^T;      _t      _t      _t ]).

	
We thus obtained a mixed integer conic quadratic problem (MICQP) which can be formulated as follows: 
	
    .[        max_ξ,,̌r,a,b              ^Tξ - r;                 s.t.     ξ≤,̱ ξ≥_n, ξ∈ℤ^n;                             2√(2)a-2r =1,;                             2√(2)b+2r =1,;                                 ^Tξ-=̌_t,;                      (a,b,^̌T)^T ∈ Q^2+t. ].

	
  Note that the formulation (<ref>) has a linear objective function and does not require potentially hugen ×nmatrixat all; it only requires then ×tmatrixwhich is often much smaller in optimum design problems. Moreover, the number of variables of (<ref>),n+t+3, is only marginally larger than the numbernof variables in the direct integer quadratic formulation (<ref>). Indeed, fort ≪nour numerical studies prove that the formulation (<ref>) can be significantly more computationally efficient than (<ref>), as we will demonstrate in Section <ref>.





§ MISCELLANEOUS COMMENTS


 §.§ Continuous design spaces


In some applications, it is possible to use a continuous design space, instead of a finite one. This is typical of factor experiments under the theoretical assumption that the levels of some factor can be any real numbers in a given interval. In such cases, AQuA cannot be directly applied[Of course, the same is true for a multitude of other popular design algorithms which work only on finite spaces.].
However, a straightforward strategy is to first apply AQuA to a finite subset of, and then use its result as an initial design for any constrained continuous optimization method which adjusts the positions of the support points within. Note that the search for optimal positions of design points in a continuous space is generally a highly non-convex problem, and a good initial feasible solution provided by a finite-space method such as AQuA can make a crucial difference. See Subsection 5.1 in <cit.> that demonstrates this approach for the full quadratic model (m=15) with4continuous factors, i.e.,=[-1,1]^4. 




 §.§ Quadratic approximations of different versions of the same criterion


We can regard criteriaΦ^+_pandΦ^-_pas part of a larger class of concave criteria: forγ∈[-1,1]and forp ∈_0we can define 

    Φ^(γ)_p:=(1+γ)Φ^+_p/2+(1-γ)Φ^-_p/2,
 where we set0 ×∞=0.
Thus,Φ^+_p=Φ^(+1)_pandΦ^-_p=Φ^(-1)_pfor allp ∈_0. Clearly,Φ^(γ)_pis a concave version of the same criterion for allγ∈[-1,1]and its quadratic approximation is 

    Φ^(γ)_pQ=(1+γ)Φ^+_pQ/2+(1-γ)Φ^-_pQ/2.

Note that settingp=0andγ=γ_d=1-d^2/1+d^2, whered=(_*)^1/m, leads to the optimization problem of the form (<ref>) with

    _0=_m^Tvec(_*^-1)

and

    _0=1/4_m^T[_*^-1⊗_*^-1]_m.

It is straightforward to verify that this choice of_0and_0corresponds to the same quadratic approximation as the one that can be obtained from theD-optimality criterion in the formlog(()), used in <cit.>. Note that we always haveγ_d ∈(-1,1). That is, in the sense of the AQuA approach, thelogcriterion is always “between” the positive and the negative versions ofD-optimality. 


Different versions of the same criterion lead to different quadratic approximations. Nonetheless, our numerical observations suggest that the differences are minor (see Subsection <ref>).



 §.§ Generalization of I-optimality and its conversion into A-optimality

  
Recently, there has been much interest inI-optimality[This criterion is sometimes called called IV- or V-optimality (see Section 10.6 in <cit.>).], becauseI-optimality may be a more appropriate criterion thanD-optimality if we are interested in the estimation of the mean value of the response (see, e.g., <cit.>, <cit.>, and <cit.>).
  
The results forA-optimality can be easily adapted to computeI-optimal designs. StandardI-optimal designs are applied to models with one-dimensional observations (r=1), and they minimize the integral of the variances of the BLUEs of the response surface over a region𝔜with respect to some measure. We will generalize the notion ofI-optimal design to potentially multivariate observations and show thatI-optimal designs areA-optimal in a transformed model, giving us the possibility to use the theory and algorithms developed forA-optimality. 


 Let𝔜 ⊆ℝ^dbe a measurable set representing a region of prediction interest, and letηbe a measure on𝔜. Suppose that for each∈𝔜, there is a matrix() ∈𝒮^m_+such that(^-1())is a measure of variance of the response surface estimator in, provided that the information matrix for the parameters is∈𝒮^m_++. For a positive definite, we can define a (generalized)I-optimality criterion

    Φ_I()=-∫_∈𝔜(^-1())dη()
      =-(^-1Ł),

whereŁ=∫_∈𝔜()dη(), and for a singular, we can setΦ_I()=-∞. Suppose thatŁ=𝐒𝐒^T, where𝐒is non-singular. Then, clearly, a designξisI-optimal if and only if it isA-optimal in the model given by the elementary information matrices 

    𝐒^-1_̋1(𝐒^T)^-1,…,𝐒^-1_̋n(𝐒^T)^-1.
 
 The standard situation corresponds tor=1,𝔜=,()=()^T(), andηbeing a uniform measure on.


We demonstrate the computation ofI-optimal designs using AQuA in Subsections <ref> and <ref>.




 §.§ Iterative application of AQuA


The central idea of this paper is to apply integer quadratic programming to a problem constructed on the basis of the optimal approximate information matrix_*, which is often available, either theoretically or via an efficient algorithm of convex optimization. Note, however, that the approximation is quite precise even if the criterion is based on a matrix_*(henceforth called the 'anchor matrix') which is not perfectly optimal. Thus, in more difficult situations, in particular with a large design space and complex design restrictions, when_*may be difficult to compute, we suggest to apply the following heuristic iterative scheme, similar to the successive application of the Newton's method to sequential quadratic optimization: 

  *  Compute a rough estimate ^(0)_* of _* at a random subsample of  or neglecting some design constraints. Set j to 0.
  
  *  Use AQuA with the anchor matrix ^(j)_* instead of _*.[If this is not last iteration of the algorithm, we can use AQuA without the integer constraints on the design. Indeed this iterative approach can also be used for computing optimal approximate designs, but we do not explore this possibility here.] Set ^(j+1)_* to be the information matrix of the resulting design.
  
  *  If a stopping rule is not satisfied, increase j by one, and continue with the previous step.


The previous scheme uses a sequence of successive quadratic optimization problems, which, in some cases, can be solved via the conic formulation of AQuA, despite the fact that we cannot solve the original optimal approximate problem because of its size or complexity. In the last subsection of the next section we will demonstrate that this approach can indeed lead to efficient EDs for large design spaces.  



 §.§ Current limitations of AQuA


AQuA can be a valuable tool in the toolbox of computational methods of experimental design as numerically demonstrated in Section <ref>. However, it has currently no theoretical underpinnings in the sense of lower bounds on the efficiency of the resulting designs depending on general properties of the problem at hand[Note that after we already have a candidate exact design for a specific problem, we can compute a lower bound on its efficiency relative to the optimal approximate design. This often leads to a guarantee which is fully satisfactory for practical purposes. Moreover, many optimization heuristics which are eminently useful across sciences also lack theoretical bounds on the efficiency of the results that they generate.] Note that we have observed that AQuA sometimes produces significantly suboptimal designs for small design sizesN ≥m, in particular forN=m[See <cit.> for an example a strongly suboptimal result of AQuA for N=m.]. Moreover, we do not have a theoretical proof of convergence of the sequential approach outlined in Subsection <ref>.


With easily available hardware and software, the IQP formulation of AQuA can solve problems with middle sizen(up to thousands) and anym ≤n. On the other hand, the MICQP version of AQuA can solve “tall” problems with a largen(up to hundreds of thousands) and a relatively smallm ≪n. However, we do not know how to use AQuA to handle problems with bothnandmlarge. 





§ NUMERICAL STUDIES


The principle of AQuA can be applied to a wide spectrum of optimal design problems in various creative ways. Here we will choose several very different examples to inform the reader about general properties of AQuA, for instance:

  *  the degree of reliability in achieving the optimal ED and the robustness with respect to the anchor matrix;

  *  the possibility to efficiently construct solutions to optimal ED problems with complex constraints on the structure of the design;

  *  the possibility to apply the conic version of AQuA to specific problems with a large design space, in particular to the problem of an information-based sub-selection of “tall” datasets.


We will demonstrate the application and explore the performance of AQuA in the  computing environment (<cit.>) employing the packages  (<cit.>),  (<cit.>), and the mathematical programming solvers of  (<cit.>). Note that there are also several other professional solvers that can handle IQP and MIQCP problems, for instance  (<cit.>). The examples were computed on a 64-bit Windows 10 system with an Intel Core i5-5500U processor at 2.40 GHz and 8 GB of RAM. The codes and additional information can be found at

http://www.iam.fmph.uniba.sk/ospm/Harman/design/http://www.iam.fmph.uniba.sk/ospm/Harman/design/.

For the application of the provided R codes, the user only needs to create the model (the matrix of all possible regressors(x)), the constraints (in the form ofand$̱), and choose the criterion (D, A or I).




 §.§ Size-constrained D- and A-optimal exact designs for the model of spring balance weighing of 6 items


Consider the linear regression of the first degree without an intercept term on the vertices of the m-dimensional unit cube given by the formula

    E(Y())=x_1β _1+…+x_mβ_m,

where the components x_j of  are chosen to be either 0 or 1. In (<ref>), the measurement Y() can be interpreted as the result of the weighing of items with unknown weights β_1,…,β_m on a spring balance, where x_j denotes the presence or the absence of the item j. Here, the design space is the set of n=2^m vertices of the unit cube in ^m. For this example, we selected m=6 items, that is, n=64.
The AD theory for model (<ref>) with the standard constraint on the size of the experiment has been worked out in great detail: see, e.g., <cit.>, who used the equivalence theorem to find Φ_p-optimal ADs for all values of p. For the application of AQuA, we can use the well-known “neighbor vertex” D-optimal and A-optimal ADs as described in <cit.>, Sec. 14.10. For non-normalized ADs of size N, and for s∈[0,m], the neighbor vertex design is

    ξ_s=(1-(s-⌊ s ⌋))ζ_⌊ s ⌋ + (s-⌊ s ⌋)ζ_⌊ s ⌋+1,

where ζ_j is a j-vertex design, i.e., ζ_j assigns N / mj to the vertices of  having j components equal to 1 and m-j components equal to 0 and ⌊ s ⌋ denotes the largest integer not exceeding s. For our case of m=6, the design ξ_s with s=24/7 is D-optimal, its support size is 35 and its information matrix is ^*_D=2N/7_6+2N/7_6. Similarly, the design ξ_s with s=3 is A-optimal, its support size is 20 and its information matrix is ^*_A=3N/10_6+2N/10_6.


In this model, the optimal ADs are not unique; the designs from Tables <ref> and <ref> are evidently not neighbor vertex designs, yet they are D- and A-optimal, respectively, which can be directly verified. Notice that the D-optimal approximate design from Table <ref> is evidently a D-optimal exact of size N=7k, k ∈ℕ, and the A-optimal approximate design from Table <ref> is evidently an A-optimal exact design of size N=10k, k ∈ℕ.[We stress that it is not completely trivial to find these balanced small-support D-, and A-optimal ADs in class of all optimal ADs; in fact, we have found them using the integer programming capabilities of AQuA. In this respect, AQuA can be very useful also for the problem with a single size constraint.]


We remark that, according to our experience, for a problem of optimal ED constrained only by the experimental size, well implemented heuristics such as the KL-exchange algorithm (<cit.>, Section 12.6) will often outperform methods based on IP solvers, including AQuA, in terms of time required to achieve a practically optimal design. However, the existing heuristics and theoretical results for the selected size-constrained problem provide benchmarks that can be used to assess the properties of the AQuA method, as we show next.


For the numerical study of ED, we will use the experimental sizes of N=6,7,…,30. For m=6, the D-optimal EDs are theoretically known (see <cit.>). For A-optimality and m=6 items, we are not aware of any publication which provides optimal EDs; therefore, we have computed the A-optimal EDs using the KL heuristic. We tested the AQuA approach realized by the integer quadratic solver of gurobi against the exact optimal values. To anchor the quadratic approximations, we used either the theoretically known optimal approximate information matrix _*, or a perturbed information matrix _* that corresponds to a random design with efficiency 0.95. The sub-optimal anchor matrix allows us to assess the robustness of the AQuA approach for problems where precise optimal AD is unavailable.  


The results, visualized in Figures <ref> and <ref>, can be summarized as follows:


  * If _* is precise (see the top panels of Figs. <ref>, <ref>), AQuA usually provides not only good, but perfectly optimal EDs. Less efficient results tend to occur for smaller sizes of N, in particular for N=m. 

  * The time to compute the solution generally increases with N (see the right panels of Figs. <ref>, <ref>). However, if there is an optimal AD that coincides with an optimal ED of a given size, the computation tends to be rapid, in particular if _* is precisely computed.

  * AQuA is generally robust with respect to the choice of the anchor matrix (see the bottom panels of Figs. <ref>, <ref>). Even using a significantly sub-optimal anchor matrix _*, the resulting EDs are either perfectly optimal or reasonably efficient, without a significant increase of the computation time (except a few specific cases of N as discussed in the comment above).

  * There are some numerical differences between the two approximations of the D- and A-criteria, but they do not tend to be pronounced.


Note that the reported computation time corresponds to the moment at which the solver determines that its current design is good enough with respect to the quadratic criterion[We did not alter the default stopping rules and other options of the gurobi solver.]; the actual time that the solver first obtains the resulting design may be shorter.


It is also worth noting that the standard ER procedure cannot be applied to the neighbor vertex optimal ADs, for N < 35 in case of D-optimality, and for N<20 in case of A-optimality. The reason is that the neighbour vertex ADs have too many support points for ER to be applicable. Even in the remaining cases that ER can be applied, for instance if we used some auxiliary tools to obtain optimal ADs with a smaller support (such as those in Tables <ref> and <ref>), our computational experience suggests that the resulting EDs tends to be worse than those found by AQuA.













 §.§ Marginally constrained symmetric D- and I-optimal exact designs for the 3-component Scheffé mixture model


The most important applications of AQuA can be expected in those situations for which there are no specialized heuristics, such as for problems with complex constraints on the design weights[We would like to stress that here we do not focus on the constraints on the design region, which are trivial to incorporate (at least in the case of finite design spaces); we work with constraints on the design vector itself in the polyhedral set of designs in ^n.]. For models with general linear constraints, few options are available if one wishes to find informative EDs. Namely, <cit.> have shown that the D- and A-optimal EDs under general linear constraints on the design weights can be obtained by solving a specific mixed integer second-order cone programming problem (MISOCP). This approach, although, given enough time, it is guaranteed to find a perfectly optimal ED, is practically feasible only for problems that are small to medium in size (with currently common hardware, up to a thousand design points even with m<10). In this subsection, we will demonstrate that our approach can be superior to both the method of <cit.>, as well as the direct application of a quadratic approximation as suggested in <cit.>.


Consider a mixture of three components with the ratio of each varying between 0% and 100% in increments of 2.5%. The response can be modeled by a quadratic Scheffé mixture model given by

    E(Y())=∑_j=1^3 β_j x_j+∑_u<vβ_(uv)x_u x_v,

where =(x_1,x_2,x_3), x_j∈{0,0.025,0.05,…,1}, j=1,2,3. Hence, the model contains m=6 unknown parameters, and the dimensionality of the set of designs n=861 (for more details and applications of mixture designs see, e.g., <cit.> and <cit.>).


Suppose that, in addition to the size constraint, we are required to compute a design that fulfils a set of marginal constraints which require that each level of each factor can be used at most once, i.e., for all permissible designs ξ and all =(x̃_1,x̃_2,x̃_3) ∈ we have ∑_x_2,x_3ξ(x̃_1,x_2,x_3) ≤ 1, ∑_x_1,x_3ξ(x_1,x̃_2,x_3) ≤ 1 and ∑_x_1,x_2ξ(x_1,x_2,x̃_3) ≤ 1. These “non-collapsibility” constraints can be justified similarly as the Latin hypercube designs (<cit.>) and “bridge” designs (<cit.>) on cubes. In particular, they lead to designs without replications of design points, which is important for computer experiments. Additionally, we have imposed  constraints of the form ξ(x_1,x_2,x_3)=ξ(x_2,x_3,x_1)=ξ(x_3,x_1,x_2) that force the design to be symmetric. Therefore, we aim to find an optimal exact  design which combines properties of non-collapsibility of individual factor levels, symmetry, and efficiency of parameter estimation. 


In this setting, we computed the D-, and I-optimal exact designs with the MISOCP approach of <cit.> and with AQuA, realized by both the standard IQP solver and by the MICQP as proposed in Subsection <ref>. The results are depicted in Figures <ref> and <ref> and described in Tables <ref> and <ref>. 


We see that all three methods of computing EDs provide designs of similar efficiency, but the conic reformulation of AQuA can decrease the computation time by as much as two orders of magnitude.  


Note that in this case, ER method cannot be used at all to transform AD to ED. Besides the support of ADs in this model being very large, as can be seen in Figures <ref> and <ref>, the marginal and symmetry constraints cannot be incorporated into ER without its significant modification. 













 §.§ D- and I-optimal subsampling of a dataset under an upper constraint on budget and lower constraint on average quality


Lastly, we will show that the conic specification of AQuA can be used as a tool for computing EDs for large design spaces, in particular for a constrained information-based subsampling of “tall” datasets; see, e.g., <cit.> for a justification of this approach. Here, the purpose is to select a subsample for a screening with the quality based on a linear regression model. In contrast to the existing information-based subsampling methods, we can require a subsample that keeps limits on the numbers of selected objects within given strata, and, simultaneously, a lower constraint on the quality as well as an upper constraints on the price of the subsample. 


To this end we used the wine datafile <cit.> that contains data on approximately 150000 wine reviews from WineEnthusiast. The aim is to subsample this database for a survey, marketing or educational purposes.


After cleaning duplicities and incomplete entries, we were left with n=111534 wine reviews containing variables on country of origin, description, points, price, province, title, variety, and winery. Out of this dataset, we are to sample wines so that the upper bound on the combined price of the wines is 1000$, lower bound on the average quality points is 90 and there is exactly one wine from each of the 42 countries. To avoid selecting the same wine more than once, we added the upper bounds ξ_i ≤ 1 for all i=1,…,n. The model used was the linear regression with m=3 parameters with intercept, the quality points and the logarithm of the price as explanatory variables.


We will use the robustness of AQuA with respect to the selection of the anchor matrix and the iterative approach explained in Subsection <ref>. For the computation of the first AD, we used the SOCP formulation from <cit.> applied to a random sub-selection comprising 1500 data-points. Then, we sequentially applied AQuA until convergence.


For both criteria, we run the randomly initiated computation 5 times, and in every case it converged to the same solution in as few as 4 steps (including the first, SOCP computation), each taking less than 4 seconds. The resulting 42-element subsamples are visualized in Figures <ref> and <ref>. It turns out that for both resulting subsamples, the cumulative price of the 42 wines is exactly 1000$ and the average quality is exactly 90 points. To meet the restrictions, the subsample computed using the D-optimality criterion is automatically concentrated largely in the region of inexpensive wines of a good quality while still making the samples diverse enough to permit precise estimation of parameter of the linear model. The result based on I-optimality is similar, but since I-optimality minimizes the average variance throughout all points, the subsample is more concentrated in the area that is most densely populated.   


We also remark that, using a standard computer, AQuA based on the IQP solver (unlike the specific MICQP formulation of AQuA) cannot be applied to design spaces of size larger than a few thousands, because of the quadratic memory requirements. That is, here we again demonstrated the advantage of the proposed conic AQuA approach over the approach of AQuA from <cit.>, not only over methods directly based on a MISOCP formulation of the problem as in <cit.>.


 

 	





§ CONCLUSIONS


  We extended the quadratic approximation from a single version of the D-criterion used in <cit.> to two versions of all Kiefer's criteria with an integer parameter, including the criterion of A-optimality and, via transformation, to the criterion of I-optimality.
  Importantly, we also proved a low-rank property of the associated quadratic forms and used it to construct efficient conic formulation of the integer quadratic programming problem. The formulation permits using the method of AQuA in case of large design spaces that are out of reach of the previous methods. On the other hand, for smaller design spaces (provided that m ≪ n is still satisfied), the conic formulation of AQuA can significantly speed up the computation; in particular, it can rapidly identify an optimal exact design in case where one of the optimal approximate designs is also optimal exact. Moreover, using AQuA it is possible to obtain efficient exact designs for situations with simultaneous constraints on various characteristics of the experiment, e.g., its form, cost, and quality.
  The basic AQuA approach presumes the knowledge of the optimal approximate information matrix. However, because the algorithms for computing optimal approximate designs are well developed and fast, this is not considered to be a drawback. Moreover, there is a large body of literature that provides theorems that explicitly yield optimal approximate designs. Note that with rounding procedures alone, the practical value of optimal approximate designs is weaker since direct heuristic computational methods can often find better designs, entirely circumventing approximate design theory and computation. We prove that optimal approximate designs carry more useful information for the construction of exact designs than is utilized by rounding procedures.
  We also showed that the AQuA approach is generally robust with respect to the misspecification of the optimal information matrix, and can even be used sequentially, starting from an anchor matrix that is far from the approximate optimum. 


  Finally, the approach of AQuA can be extended to various criteria other than those analyzed in this paper; what is needed is only their quadratic approximation[For the application of the conic improvement, the quadratic forms must have low ranks.], which can be found either analytically or numerically. This opens up new possibilities for the computation of optimal experimental designs with respect to criteria that are difficult to evaluate.


Acknowledgments The work was supported by Grant No 1/0341/19 from the Slovak Scientific Grant Agency (VEGA).





50
	
[Anderson-Cook et al. (2009)] ABM Anderson-Cook CM, Borror CM, Montgomery DC: Response surface design evaluation and comparison. Journal of Statistical Planning and Inference 139(2), 629-641 (2009)
			
[Atkinson et al. (2007)]Atkinson Atkinson AC, Donev AN, Tobias RD: Optimum Experimental Designs, with SAS. Oxford University Press, Oxford (2007)

[Bouhtou et al. (2010)]Bouhtou Bouhtou M, Gaubert S, Sagnol G: Submodularity and randomized rounding techniques for optimal experimental design. Electronic Notes in Discrete Mathematics 36, 679-686 (2010)

[Cheng (1987)]Cheng Cheng C: An application of the Kiefer-Wolfowitz equivalence theorem to a problem in Hadamard transform optics. The Annals of Statistics, 15(4), 1593-1603  (1987)

[Cornell (2011)]CornellCornell JA: Experiments with mixtures: designs, models, and the analysis of mixture data (Vol. 403). John Wiley & Sons (2011)

[Dattorro (2008)]Dattorro Dattorro, J: Convex optimization & Euclidean distance geometry. Meboo Publishing (2008)

[Dykstra (1971)]Dykstra Dykstra O: The augmentation of experimental data to maximize |X'X|. Technometrics 13, 682-688 (1971)

[Filová and Harman (2018)]aqua Filová L, Harman R: Ascent with Quadratic Assistance for the Construction of Exact Experimental Designs. arXiv preprint arXiv:1801.09124v2 (2018)

[Goos and Jones (2011)]GoosJones Goos P, Jones B: Optimal design of experiments: a case study approach. John Wiley & Sons (2011)

[Goos et al. (2016)]GJS Goos P, Jones B, Syafitri U: I-Optimal Design of Mixture Experiments. Journal of the American Statistical Association 111,  899-911 (2016)

[Gurobi Optimization, LLC (2018)]gurobi Gurobi Optimization, LLC: Gurobi Optimizer Reference Manual. <http://www.gurobi.com> (2018)

[Haines (1987)]Haines Haines LM: The application of the annealing algorithm to the construction of exact optimal designs for linear-regression models. Technometrics 29, 439-447 (1987)
		
[Harman and Filová (2014)]HF Harman R, Filová L: Computing efficient exact designs of experiments using integer quadratic programming. Computational Statistics & Data Analysis 71, 1159-1167 (2014)
		
[Harman and Filová (2016)]RLIB Harman R, Filová L: Package 'OptimalDesign'.  <https://cran.r-project.org/web/packages/OptimalDesign/index.html> (2016)
		
[Harman et al. (2016)]HBF Harman R, Bachratá A, Filová L: Construction of efficient experimental designs under multiple resource constraints. Applied Stochastic Models in Business and Industry 32(1), 3-17 (2016)

[Harman et al. (2019)]REX Harman R, Filová L, Richtárik P: A Randomized Exchange Algorithm for Computing Optimal Approximate Designs of Experiments. Journal of the American Statistical Association, to appear (2019)

[Harman and Trnovská (2009)]HTHarman R, Trnovská M: Approximate D-optimal designs of experiments on the convex hull of a finite set of information matrices. Mathematica Slovaca, 59(6), 693-704 (2009)

[Harville (1997)]Harville Harville, DA: Matrix algebra from a statistician's perspective (Vol. 1). New York, Springer (1997)
		
[Jones et al. (2015)]bridge Jones B, Silvestrini RT, Montgomery DC, Steinberg DM: Bridge designs for modeling systems with low noise. Technometrics, 57(2), 155-163 (2015)
				
[Kiefer (1971)]kiefer Kiefer J: The role of symmetry and approximation in exact design optimality. Statistical Decision Theory and Related Topics, 109-118 (1971)
		
[Liu and Neudecker (1995)]LNLiu S, Neudecker H: A V-optimal design for Scheffé's polynomial model. Statistics & Probability Letters, 23(3), 253-258 (1995)
		
[Mandal et al. (2015)]MWY Mandal A, Wong WK, Yu Y: Algorithmic Searches for Optimal Designs. In Dean, A., Morris, M., Stufken, J., Bingham, D.:  Handbook of Design and Analysis of Experiments,  211-218, Chapman and Hall (2015)
		
[McKay et al. (1979)]LHD McKay MD, Beckman RJ, Conover WJ: A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 21(2), 239-245 (1979)

[Montgomery (2017)] Mont Montgomery DC: Design and analysis of experiments. John Wiley & Sons (2017)

[MOSEK ApS (2018a)] moseka Mosek modeling cookbook, Release 3.1., available online at <https://docs.mosek.com/MOSEKModelingCookbook-letter.pdf> (2018)

[MOSEK ApS (2018b)] mosekb The MOSEK optimization toolbox for MATLAB manual. Version 8.1., <http://docs.mosek.com/8.1/toolbox/index.html> (2018)

[Neubauer et al. (2000)]NWZ Neubauer MG, Watkins W, Zeitlin J: D-optimal weighing designs for six objects. Metrika, 52(3), 185-211 (2000)

[Novomestky (2012)]matrixcalc Novomestky F: matrixcalc: Collection of functions for matrix calculations. <https://cran.r-project.org/package=matrixcalc> (2012)
		
[Pázman (1986)]Pazman Pázman A: Foundations of optimum experimental design. Reidel, Dordrecht (1986)
		
[Pukelsheim (2006)]Puk Pukelsheim F: Optimal Design of Experiments (Classics in Applied Mathematics). SIAM (2006)
		
[Pukelsheim and Rieder (1992)]PukRieder Pukelsheim F, Rieder S: Efficient rounding of approximate designs. Biometrika, Vol. 79, 763-770 (1992)
		
[R Development Core Team (2011)]R R Development Core Team: R: A Language and Environment for Statistical Computing. Foundation for Statistical Computing, Vienna, Austria (2011)
		
[Sagnol (2013)]Sagnol Sagnol G: Approximation of a maximum-submodular-coverage problem involving spectral functions, with application to experimental designs. Discrete Applied Mathematics 161(1), 258-276 (2013)		

[Sagnol and Harman (2015)]SagnolHarman Sagnol G, Harman R: Computing exact D-optimal designs by mixed integer second-order cone programming. The Annals of Statistics 43(5),  2198-2224 (2015)

[Stufken and Yang (2012)]SY Stufken J, Yang M: On locally optimal designs for generalized linear models with group effects. Statistica Sinica 22(4), 1765-1786 (2012)

[Thoutt (2017)]wine Thoutt Z: Wine Reviews Data. <https://github.com/zackthoutt/wine-deep-learning> (2017)
		
[Yang et al. (2013)]YBT Yang M, Biedermann S, Tang E: On optimal designs for nonlinear models: a general and efficient algorithm. Journal of the American Statistical Association 108, 1411-1420 (2013)

[Yu (2011)]Yu Yu Y: D-optimal designs via a cocktail algorithm. Statistics and Computing 21(4), 475-481 (2011)
		
[Wang et al. (2018)]iboss Wang H, Yang M, Stufken J: Information-based optimal subdata selection for big data linear regression. Journal of the American Statistical Association, 1-13 (2018)


	
