Optimal designs of experiments are developed to reduce the variance of estimated model parameters [24] by optimizing a function of the information matrix usually under a (generalized) linear model assumption. Examples include the determinant, the trace, and the minimum eigenvalue of the information matrix, which result in D-, A-, and E-optimal designs, respectively. In many areas, such as clinical trials or user behavior studies in IT companies, the experiments are conducted to investigate the treatment effects over different users by designing controlled experiments, i.e., divide the users into different groups and assigning each group with a treatment or control. For example, optimal design of experiments is applied to design clinical trials in the era of precision medicine, where the patient response may depend on some biomarkers (covariates). In fact, there is a significant growing evidence that the response of the patients may depend on their covariates in a wide range of applications. For instance, [26] showed that the top ten highest-grossing drugs in the US are only effective in 4% to 25% of the patient population, a shocking result for the current clinical practice. [6] applied a modified D-optimal design to early stage clinical trials when the covariates are uncontrollable. [19] applied a weighted L-optimal design for late stage Phase III clinical trials with observed covariates. As another emerging application, Tech Giants such as Google, Facebook, and LinkedIn frequently use A/B testing for marketing, web design, and data-driven decision making [29]. In fact, applying controlled experiments such as A/B testing experiments in business has improved the performance of companies [18]. Also, there are some methodological contributions for optimal design of experiments applied for A/B testing. For example, [5] studied an offline and online standard A/B testing with the goal of maximizing the precision of least square estimations.

The literature of optimal design of experiments for linear models, which we also use in this study, is vast. From the optimal design perspective, [2] investigated the problem with two treatments, and developed the D-optimal design for both covariates and treatment allocation under the assumption of a linear model with interactions between covariates and treatments. [30] generalized the work of [2] to experiments with multiple treatments. Following the same model setting, [25] developed the A- and E-optimal designs for treatments and covariates.

In practice, the user covariates may not be controllable, due to the restriction in the user recruiting process. Thus, an alternative problem is to allocate the treatments to a given set of users with observed covariates. From the optimal design perspective, the aim is to obtain an exact design by optimizing a design objective with the treatment allocation as the design variable given fixed covariates, which often results in an optimization problem with binary decision variables. For two treatment allocation, under the linear model assumption without interaction between covariates and treatments, [16] proposed computationally tractable algorithms to search D- and A- optimal designs. Also, with a particular interest of minimizing the variance of the estimated treatment effect (i.e., a D s optimal design criterion), [5] demonstrated that the problem can be solved efficiently based on a generalization of the MAX-CUT semi-definite programming (SDP) relaxation of [14]. Given the linear model with interaction between treatments and covariates, [31] proposed an optimal design objective for personalized decision making with two treatments and a corresponding approximation solution approach to obtain near optimal allocations efficiently.

In this paper, we generalize the work in [31] to multiple treatment allocation for personalized decision making with observational covariates. Specifically, we consider the situation that there are a finite number of treatment options and subjects with given covariates. We assume that the response of a subject assigned to a treatment follows a linear model which includes the interaction between covariates and treatments. Accordingly, we can define the best treatment option for each subject groups, which can be estimated with data. This setup facilitates precision decision making. However, by including more than two treatments, the resulting design objective can not be simplified by the approximation approach developed for the two treatment cases in [31]. Since the design objective matches the E-optimal design criterion, we reformulate the optimization problem as a SDP with binary decision variables. We use a YALMIP and MOSEK based optimization solver for SDP in a branch-and-bound scheme to provide the optimal design. In our numerical study, we assess the quality of the optimal designs provided by the solver.

A decision-making problem in practice can often be simplified to choose a treatment from a finite candidate set { 1 , … , K } that has the best performance measure μ k , i.e., (2.1) k ∗ ∈ argmax k ∈ { 1 , … , K } μ k . The values of μ k ’s are unknown, which can be estimated based on collected user responses to their treatment allocations. Particularly, given n users, the experimenter allocates one of the K treatments to each of the users. The response of the i-th user is denoted by y i . Assume that y i is a continuous scalar that follows a model given by (2.2) y i = ∑ k = 1 K x i k μ k + ε i , for i = 1 , … , n , where ε i is an additive error term, and x i k ∈ { 0 , 1 } with x i k = 1 indicating that the k-th treatment is allocated to the i-th user. We require that ∑ k = 1 K x i k = 1 to guarantee that each user is only exposed to a single treatment. Under the model assumption in (2.2), the performance measure μ k represents the average treatment performance over the target user population, which can be estimated by μ ˆ k given responses y i ’s and treatment allocation x i k ’s. Then a data-driven solution to (2.1) is given by (2.3) k ˆ ∗ ∈ argmax k ∈ { 1 , … , K } μ ˆ k . The accuracy of this decision is associated with the accuracy of the estimates μ ˆ k ’s. To reduce the uncertainty of the decision, we can reduce the variances of the estimates μ ˆ k ’s.

Under the context of personalized decision making, a user is associated with observed covariates z = ( z 1 , … , z p ) ⊤ ∈ Z ⊂ R p , where Z is the covariates space of the target population. Examples of the covariates include the demographic information, social behavior and network connections. The treatment performance can often be related to the covariates values. Let μ k ( z ) be a function of the users’ covariates z representing the personalized treatment performance measure. In this paper, we further assume that μ k ( z ) is a linear function over z. Therefore, the linear model in (2.2) is represented by (2.4) y i = ∑ k = 1 K x i k μ k ( z i ) + ε i = ∑ k = 1 K x i k z i ⊤ β k + ε i , for i = 1 , … , n , where β k = ( β 1 k , … , β p k ) ⊤ is a vector of linear coefficients. The first feature in each z i is fixed to be one as the linear intercept. Similar to [30], this model incorporates the interaction between treatment and covariates, which enables the estimation of heterogeneity of treatment performance over covariates z.

Under this model, the problem of personalized decision making becomes (2.5) k ∗ ( z ) = argmax k ∈ { 1 , … , K } μ k ( z ) = argmax k ∈ { 1 , … , K } z ⊤ β k for any z ∈ Z . With estimated parameters, a data-driven solution to (2.5) can be obtained for any z ∈ Z, i.e., (2.6) k ˆ ∗ ( z ) = argmax k = 1 , 2 , … , K z ⊤ β k ˆ Similar to (2.3), to improve the accuracy of personalized decisions, it is desired to reduce the variance of estimated μ k ( z ) = z ⊤ β k for each k ∈ { 1 , … , K } and each z ∈ Z. Following the optimal design perspective in [31], we can reduce the variances by searching an optimal allocation of the treatments x i k ’s in the presence of fixed covariates of recruited users as illustrated in Figure 1. Next, we investigate the design criterion that characterizes the variance of estimated μ k ( z ) = z ⊤ β k for each k ∈ { 1 , … , K } and any z ∈ Z.

Let y = ( y 1 , … , y n ) ⊤ be the responses of n users. Correspondingly, let Z be an n × p matrix, whose rows are the users’ covariates z 1 ⊤ , … , z n ⊤ (the first entry is loaded by one as the intercept). The treatment allocation is recorded in an n × K matrix X, where the i-th row of X is denoted by x i ⊤ = ( x i 1 , … , x i K ) ⊤ representing the treatment allocation of the i-th user. We assume that the responses are generated under the model in (2.4). The model covariates matrix is then given by x 1 ⊤ ⊗ z 1 ⊤ ⋮ x n ⊤ ⊗ z n ⊤ which we assume has a full column rank. The notation ⊗ denotes the Kronecker product. By stacking β k ’s in a vector b = ( β 1 ⊤ , … , β K ⊤ ) ⊤ of size K p, we can express the least squares estimator of b by (3.1) b ˆ = β ˆ 1 ⊤ ⋯ β ˆ K ⊤ = x 1 ⊤ ⊗ z 1 ⊤ ⋮ x n ⊤ ⊗ z n ⊤ ⊤ x 1 ⊤ ⊗ z 1 ⊤ ⋮ x n ⊤ ⊗ z n ⊤ − 1 · x 1 ⊤ ⊗ z 1 ⊤ ⋮ x n ⊤ ⊗ z n ⊤ ⊤ y .

We assume that the additive error term ε i ’s in (2.4) are independent and identically distributed random variables with mean zero and variance σ 2 . The variance-covariance matrix of the estimated b in (3.1) can be expressed by var ( b ˆ ) = σ 2 x 1 ⊤ ⊗ z 1 ⊤ ⋮ x n ⊤ ⊗ z n ⊤ ⊤ x 1 ⊤ ⊗ z 1 ⊤ ⋮ x n ⊤ ⊗ z n ⊤ − 1 . Notice that x 1 ⊤ ⊗ z 1 ⊤ ⋮ x n ⊤ ⊗ z n ⊤ ⊤ x 1 ⊤ ⊗ z 1 ⊤ ⋮ x n ⊤ ⊗ z n ⊤ = ∑ i = 1 n ( x i ⊗ z i ) ( x i ⊤ ⊗ z i ⊤ ) = ∑ i = 1 n ( x i x i ⊤ ) ⊗ ( z i z i ⊤ ) = ∑ i = 1 n z i z i ⊤ x i 1 0 0 ⋱ 0 0 ∑ i = 1 n z i z i ⊤ x i K , which is the result of that x i x i ⊤ is a K × K diagonal matrix with diagonal entries x i 1 , … , x i K . Then the variance-covariance matrix of b ˆ can be simplified to a block diagonal matrix: (3.2) var ( b ˆ ) = σ 2 diag ∑ i = 1 n z i z i ⊤ x i 1 − 1 … ∑ i = 1 n z i z i ⊤ x i K − 1 .

For any k ∈ { 1 , … , K } and z ∈ Z, we have that var z ⊤ β ˆ k = z ⊤ var β ˆ k z = σ 2 z ⊤ ∑ i = 1 n z i z i ⊤ x i k − 1 z . As noted earlier, we assume that the experiments are conducted with n existing users. Thus, the personalized information matrix Z is observed, and the experimental design problem is to determine the treatment allocation matrix X. A necessary condition for positive definiteness of the variance-covariance matrix requires that ∑ i = 1 n x i k ≥ p for any k ∈ { 1 , … , K } and n ≥ K p.

Our goal is to minimize the variance of z ⊤ β ˆ k for any k and z to improve the accuracy of the data-driven decision in (2.6). Instead of minimizing all the variances simultaneously, we formulate the optimization problem by minimizing the worst case, i.e., the maximum variance of any k ∈ { 1 , … , K } and z ∈ Z. Therefore, the design criterion is given by (3.3) max z ∈ Z max k = 1 , 2 , … , K var z ⊤ β ˆ k = max z ∈ Z max k = 1 , 2 , … , K σ 2 z ⊤ ∑ i = 1 n z i z i ⊤ x i k − 1 z . The optimal design is then given by solving the problem (3.4a) min X max z ∈ Z max k = 1 , 2 , … , K z ⊤ ∑ i = 1 n z i z i ⊤ x i k − 1 z (3.4b) s.t. ∑ i = 1 n x i k ≥ p , k = 1 , … , K , (3.4c) ∑ k = 1 K x i k = 1 , i = 1 , 2 , … , n , (3.4d) X ∈ { 0 , 1 } n × K , where ∑ i = 1 n x i k ≥ p is required for the positive definiteness of the matrix ∑ i = 1 n z i z i ⊤ x i k , and the constraint ∑ k = 1 K x i k = 1 ensures that each user is exactly assigned to one treatment.

Although the optimal design criterion has a closed-form expression, the corresponding optimization problem is still a challenge to solve directly as a mixed-integer minimax problem. We assume that the covariates space Z is given by a unit ball, i.e., Z = { z ∈ R p ∣ ‖ z ‖ ≤ 1 }. For any X that gives a positive definite variance-covariance matrix, we have that max ‖ z ‖ ≤ 1 max k = 1 , 2 , … , K z ⊤ ∑ i = 1 n z i z i ⊤ x i k − 1 z = max k = 1 , 2 , … , K max ‖ z ‖ ≤ 1 z ⊤ ∑ i = 1 n z i z i ⊤ x i k − 1 z = max k = 1 , 2 , … , K max ‖ z ‖ = 1 z ⊤ ∑ i = 1 n z i z i ⊤ x i k − 1 z = max k = 1 , 2 , … , K λ max ∑ i = 1 n z i z i ⊤ x i k − 1 = max k = 1 , 2 , … , K 1 / λ min ∑ i = 1 n z i z i ⊤ x i k = max k = 1 , 2 , … , K 1 / min ‖ z ‖ = 1 z ⊤ ∑ i = 1 n z i z i ⊤ x i k z , where λ max and λ min produce the maximum and minimum eigenvalues of a matrix, respectively. Accordingly, the optimization problem in (3.4) is equivalent to (3.5a) max X min k = 1 , 2 , … , K λ min ∑ i = 1 n z i z i ⊤ x i k (3.5b) s.t. ∑ i = 1 n x i k ≥ p , k = 1 , … , K , (3.5c) ∑ k = 1 K x i k = 1 , i = 1 , 2 , … , n , (3.5d) X ∈ { 0 , 1 } n × K . This design objective is known as the E-optimal design in the literature (e.g., [24] and [9]), i.e., min k = 1 , 2 , … , K λ min ∑ i = 1 n z i z i ⊤ x i k = λ min ∑ i = 1 n z i z i ⊤ x i 1 0 0 ⋱ 0 0 ∑ i = 1 n z i z i ⊤ x i K .

Next, we provide a SDP-based solution approach to solve (3.5) with observed covariates Z and a binary decision X.

Following the reformulation technique for E-optimal design in [8], we cast problem (3.5) as a SDP with binary decision variables. The reformulation is based on the following observation. Let A be a real p × p symmetric matrix. It is well known in matrix analysis that for any δ ∈ R, λ min ( A − δ I ) = λ min ( A ) − δ, where I is the p × p identity matrix. Therefore, λ min ( A ) can be expressed as the largest value of δ such that λ min ( A − δ I ) ≥ 0, i.e., λ min ( A ) = max { δ ∣ A − δ I ⪰ 0 } , where A − δ I ⪰ 0 means that A − δ I is a positive semi-definite matrix.

We apply the above observation to problem (3.5). For k = 1 , … , K, λ min ∑ i = 1 n z i z i ⊤ x i k = max δ k | ∑ i = 1 n z i z i ⊤ x i k − δ k I ⪰ 0 .

As a result, min k = 1 , … , K λ min ∑ i = 1 n z i z i ⊤ x i k = min k = 1 , … , K max δ k | ∑ i = 1 n z i z i ⊤ x i k − δ k I ⪰ 0 = max λ | λ ≤ δ k , ∑ i = 1 n z i z i ⊤ x i k − δ k I ⪰ 0 , k = 1 , … , K = max λ | ∑ i = 1 n z i z i ⊤ x i k − λ I ⪰ 0 , k = 1 , … , K . Consequently, problem (3.5) can be equivalently written as (4.1a) max X , λ λ (4.1b) s.t. ∑ i = 1 n z i z i ⊤ x i k − λ I ⪰ 0 , k = 1 , 2 , … , K , (4.1c) ∑ i = 1 n x i k ≥ p , k = 1 , … , K , (4.1d) ∑ k = 1 K x i k = 1 , i = 1 , 2 , … , n , (4.1e) X ∈ { 0 , 1 } n × K .

Note that, the reformulation technique is the same as used in the SDP reformulation proposed for the approximate E-optimal design in Section 7.5 of [8]. However, the discrete feature in problem (3.5) cannot be averted by simply considering a continuous relaxation as in [8], even when n is large. Thus, the implementation of a SDP problem with continuous decisions cannot solve our problem. As a mixed binary SDP, the reformulation (4.1) can be solved by a branch-and-bound (BnB) algorithm, where each relaxation problem in a BnB node is solved by an SDP solver [13]. A built-in solver in YALMIP [20] provides such a BnB implementation for MATLAB [21], and MOSEK [1] can be used as the embedded SDP solver.

In practice, the experimenter often requires a (roughly) balanced allocation across the K treatments. If K can be divided by n, we can include balanced constraints (4.2) ∑ i = 1 n x i k = n / K , k = 1 , … , K , which also potentially accelerates the computation time of solving this problem by reducing the number of feasible solutions. If K can not be divided by n, a set of roughly balanced constraints can be added similarly. For convenience, we set n as a multiple of K and include the balanced constraints in our numerical study. We lastly remark that the original problem in (3.4) is equivalent to this reformulation under the constraint that Z = { z ∈ R p ∣ ‖ z ‖ ≤ 1 }. However, for a given set of covariates, we can rescale the covariates with the largest possible norm of the covariates in Z to make sure that this constraint is satisfied.

In this section, we illustrate the numerical performance of the proposed solution approach based on synthetic datasets. Our aim is to demonstrate the quality of the resulting optimal design in terms of the optimality objectives and the prediction performance.

To improve the accuracy of personalized decision, this paper investigates the optimal design to minimize the maximum variance of estimated personalized treatment effects over different treatments and different covariates values. The resulting design objective matches the E-optimal design criterion. To provide the optimal design of multiple treatment allocation in the presence of observed covariates, a SDP solution approach is applied to solve the optimization problem. The proposed solution approach is able to provide optimal designs efficiently for n from 100–300 as shown in our numerical results. We point out potential future directions. First, it is desired to develop more efficient reformulations or optimization algorithms to handle the proposed design objective with larger sizes. Second, in practice, it is useful to investigate online optimal allocation assuming that the users are recruited in an online dynamic fashion under the proposed design objective.