# Batch mode active learning for efficient parameter estimation

Wei Zheng1, Ting Tian2, and Xueqin Wang3
1Department of Business Analytics and Statistics, University of Tennessee
2School of Mathematics, Sun Yat-sen University
3School of Management, University of Science and Technology of China

###### Abstract

For many tasks of data analysis, we may only have the information of the
explanatory variable and the evaluation of the response values are quite
expensive. While it is impractical or too costly to obtain the responses of
all units, a natural remedy is to judiciously select a good sample of units,
for which the responses are to be evaluated. In this paper, we adopt the
classical criteria in design of experiments to quantify the information of a
given sample regarding parameter estimation. Then, we provide a theoretical
justification for approximating the optimal sample problem by a continuous
problem, for which fast algorithms can be further developed with the guarantee
of global convergence. Our results have the following novelties: $(i)$ The
statistical efficiency of any candidate sample can be evaluated without
knowing the exact optimal sample; $(ii)$ It can be applied to a very wide
class of statistical models; $(iii)$ It can be integrated with a broad class
of information criteria; $(iv)$ It is much faster than existing algorithms.
$(v)$ A geometric interpretation is adopted to theoretically justify the
relaxation of the original combinatorial problem to continuous optimization
problem.

Keywords: Optimal design; active learning; sample; equivalence theorem;
precision medicine

##  1 Introduction.

For many tasks of data analysis, we may only have the information of the
explanatory variable while the response values are unavailable and expensive
to obtain. As a result, it is impractical or too costly to obtain the
responses of all units. With limited budget, we have to judiciously select a
good sample of the units to obtain the response and analyse this sample in
place of the full data. This problem has been extensively discussed in the
machine learning community by the name of active learning. See Settles (2009)
and Ren et al. (2021) for comprehensive surveys. Relatively, the study of
sampling methods targeted on statistical inferences is quite limited in
literature. Deng et al. (2009) studied the detection of money laundry activity
by using logistic regression model to select the informative bank customer in
a sequential procedure. In each step, one unit is selected from several pre-
screened candidate units according to the D-criterion and added to the sample.
The estimate of the parameters are correspondingly updated, based on which the
next point will be selected. This sequential process may not be practical when
the evaluation of the response is time consuming. Additionally, it fails to
take the information overlap between the multiple data points into account. In
such situations, we may need to select a big batch of points simultaneously
and the responses of the selected units are measured simultaneously,
especially when a parallel system is available. In Deng et al. (2009)’s
context, we shall confirm the money laundry activity for many bank customers
simultaneously. Another disadvantage of the sequential approach is that it
does not lead to the global optimization due to information overlap among
successively selected units.

The problem is formally stated as follows. Let ${\cal
X}=\\{x_{1},x_{2},...,x_{N}\\}$ be a set of all $N$ design points, and our
task is to select a sample $S\subset\\{1,2,...,N\\}$ of a given size $|S|=n<N$
to obtain the response $\\{y_{i}:i\in S\\}$. Then we shall use the data
$\\{(x_{i},y_{i}):i\in S\\}$ to train a parametric model $p(y|x,\theta)$,
where $\theta\in{\mathbb{R}}^{k}$ is the parameter of the distribution. The
information matrix for $\theta$ based on $x$ is
$M_{x}=-{\mathbb{E}}\partial^{2}log(p(y|x,\theta))/\partial\theta^{2}\in{\mathbb{R}}^{k\times
k}$. When the responses are independent and homoscedastic condition on the
explanatory variable, the information matrix of the sample $S$ is simply
$M(S)=\sum_{i\in S}M_{x_{i}}$. That means the asymptotic covariance matrix of
the maximum likelihood estimator (MLE), $\hat{\theta}$, is proportional to
$M(S)^{-1}$. For a fixed sample size, say $n$, we would like to choose the
optimal sample which maximizes $M(S)$ in some sense. Suppose we are interested
in $g(\theta)\in{\mathbb{R}}^{q}$, $1\leq q\leq k$, where $g$ is a
differential function, the asymptotic covariance matrix of its MLE,
$g(\hat{\theta})$, is given by

$\displaystyle\Sigma_{g}(S)$ $\displaystyle=$ $\displaystyle\frac{\partial
g(\theta)}{\partial\theta^{T}}M(S)^{-1}\left(\frac{\partial
g(\theta)}{\partial\theta^{T}}\right)^{T}.$

In the design context, the covariance matrix $\Sigma_{g}(S)$ can be minimized
under different criteria. For example, the criterion functions
$tr(\Sigma_{g}(S))$, $|\Sigma_{g}(S)|^{1/q}$ and
$\lambda_{max}(\Sigma_{g}(S))$ corresponds to the A, D and E criteria,
respectively. Here, $tr(\cdot)$ denotes the trace of a matrix; $|\cdot|$
denotes the determinant of a matrix; and $\lambda_{max}(\cdot)$ denotes the
largest eigenvalue of a matrix. In an effort to unify these criteria, Kiefer
(1974) defined the $\Phi_{p}$ function:

$\displaystyle\Phi_{p}(\Sigma_{g}(S))=\left(q^{-1}Tr\left[\Sigma_{g}(S)\right]^{p}\right)^{1/p},$
$\displaystyle 0\leq p<\infty.$ (1)

A sample is said to be $\Phi_{p}$-optimal if it solves the following problem

$\displaystyle\min_{S:|S|=n}\Phi_{p}(\Sigma_{g}(S)).$ (2)

By direct calculations, one can show that

$\displaystyle\lim_{p\rightarrow 0}\Phi_{p}(\Sigma_{g}(S))$ $\displaystyle=$
$\displaystyle|\Sigma_{g}(S)|^{1/q}.$

Hence, the $\Phi_{p}$ criterion reduces to D-criterion when $p=0$. Similarly,
it reduces to the A and E criteria when $p$ takes the values of $1$ and
$\infty$, respectively.

The combinatorial problem (2) is generally N-P hard, and the aim of this paper
is to develop a method of deriving a near-optimal solution in a fast manner.
The method shall also permit the evaluation of the lower bound of the
efficiency of any candidate solution. The rest of the paper is organized as
follows. Section 3 relaxes (2) to a continuous optimization problem and
provides a geometric interpretation of the connection between the two problems
under the D-criterion. Moreover, two forms of optimality conditions are
provided. Correspondingly, an algorithm is proposed for D-criterion. Section 4
generalizes the algorithm to the general $\Phi_{p}$-criterion and proves the
convergence of the algorithm. Section 5 provides some numerical examples to
illustrate that it runs a lot faster than existing algorithms.

##  2 Literature review

Given a large pool of unlabeled unit, active learning provides a way to
iteratively select the most informative unlabeled unites for the queries to
label. Most work on active learning have focused on selecting a single
unlabeled unit at a time. Among them, one main category is to measure the
informativeness of each unlabeled point. Zhang and Oles (2000) selects the
unit that provides the most information for the prediction of unlabelled data.
Houlsby et al. (2011) selects points that maximally decreases the expected
posterior entropy of the model parameters. They achieve this by adopting
Sebastiani and Wynn (2000)’s formulation on the equivalence of entropy changes
on the responses and parameters when a new unit is added to the data. This is
helpful to handle high dimensional parameters and heteroscedastic cases such
as the Guassian process. For other ideas, Lewis and Gale (1994) selects the
unit with conditional probability closest to $1/2$. Freund et al. (1997)
selects the unit on which a committee of classifiers disagree the most.
Campbell et al. (2000) and Tong and Koller (2001) suggest choosing the unit
closest to the classification boundary and the latter formulated it as a
version space reduction process. Scheffer et al. (2001), Wang and Hua (2011),
and Chen and Krause (2013) used the difference in prediction probability of
the most and second most probably class as the criterion.

The single unit selection demands for less computational resource, but also
has many drawbacks. It requires retraining with all units when each new unit
is added to the data. This is inefficient when the training is time consuming.
Since the next unit selection has to wait for the labeling of the previous
selected unit and this is a waste of resource and time when the parallel
labeling system is available. Also, there will be certain information overlap
among sequentially selected units since similar points that are ranked as most
informative or most plausible by other criteria will be selected. To address
this issue, some has proposed to select the subset of points simultaneously as
presentative examples for the whole population. One natural approach is based
on clustering (Li et al. (2012), Ienco et al. (2013)). The idea is that after
clustering is performed, the center of the clusters will be good
representative examples of the points in the cluster. Sener and Savarese
(2018) investigated the Core-set selection problem. The representative method
is more diversified than the information based method, but it tends to select
abundant points in high density area of the distribution.

Many have studied batch mode active learning (BMAL) to achieve the
informativeness and diversification goal jointly. Brinker (2003), Schohn and
Cohn (2000), and Xu et al. (2003) extended single instance selection
strategies that use support vector machines (SVM), and enforced diversity
among points within the margin of SVM. Guo and Schuurmans (2007) focused on
the semisupervised learning, where the objective function for the optimal
batch is the likelihood function of labeled data minus the entropy of the
unlabeled data expressed with the logistic regression as the working model.
The optimization is achieved by reformulated the objective function as an
integer programming problem. Hoi et al. (2006b, a) choose multiple instances
that efficiently reduce the Fisher information. Specifically, Hoi et al.
(2006a) extended Zhang and Oles (2000)’s work based on logistic regression to
the batch mode active learning and and took the criterion of the ratio of the
information matrices between all unlabeled samples and the selected unlabeled
sampled to be labeled. The ratio could be understood as the capability of
using the fitted model to predict all unlabeled samples once the selected
samples are labeled. They adopted the greedy algorithm based on the idea of
submodular function to reach the $1-1/e$ efficiency. The algorithm iteratively
optimize on each unit of the batch by fixing the remaining fixed. The
submodularity has been extensively applied in literature of active learning.
For example, Wei et al. (2015) connected it to likelihood functions of two
classifiers for the multistage BMAL. Ravi et al. (2016) studied BMAL for
lasso, where the objective function is the D-optimality plus penalty term
instead of the entropy. Kirsch et al. (2019) extended Houlsby et al. (2011) to
the batch model mutual information measurement. Deng et al. (2022) studied
semi-supervised clustering problem, where they select the most informative
pair to query based on the entropy criterion. They have all applied the
submodule function for finding the subset based on the criteria defined
therein. For deep learning models, Zhdanov (2019) used the K-means clustering
algorithm for the batch selection. Ghorbani et al. (2021) used Shapley value
to evaluate the importance of each data point. Specifically, a larger size
pool of data is preselected with high Shapley values under the KNN classifier
and a diversity algorithm is applied to further select the desired size of
batch from the pool.

Overall, many approaches use a variety of heuristics to guide the instance
selection process, where the selected batch should be informative about the
classification model while being diverse enough so that their information
overlap is minimized. For those formulating the BMAL as optimization problem
that aims to learn a good classifier directly, either an integer programing
method or the submodular algorithm with $1-1/e$ efficiency is employed to
solve the optimization problem. This is due to the NP-hard nature of subset
selection problem. Here, we adopt a new method which first relax the
combinatorial problem to a continuous optimization problem and a round
procedure is adopted to obtain the discrete solution. Mathematically, this
formulation is related to the constrained optimal design of experiments
historically. See Wynn (1977, 1982), Fedorov (1989), Müller and Pázman (1998),
Sahm and Schwabe (2001), Pronzato (2004, 2006), Uciński (2015) for example.
However, these methods are not applicable for the current BMAL problem. One
main reason is because exchange algorithm has been majorly adopted for the
optimization. This is not computational efficient in dealing with large
database. Meanwhile, models considered in literature are mostly of regression
type (e.g. logistic regression) so that the information matrix from the unit
is of the form $xx^{T}$, where $x$ is a column vector and T is the transpose
operator. Recently, the backward algorithm is proposed by Ouyang et al.
(2016a) when they search for the D-optimal sample for the logistic regression
model. In Ouyang et al. (2016b), they further applied the procedure to
precision medicine. Ouyang et al. (2017) extended the discussion to other
optimality criteria. Unfortunately, the backward algorithm is quite slow
especially when the population size increases. Besides the speed issue of
these two methods, the global statistical efficiency of the derived samples
remains unknown. Here, we have provided geometric interpretations of our
algorithm to justify the empirical evidenced close to 100% efficiency.

To summarize, our results have the following novelties. $(i)$ The statistical
efficiency of any candidate sample can be evaluated without knowing the exact
optimal sample; $(ii)$ It can be applied to a very wide class of statistical
models; $(iii)$ It can be integrated with a broad class of information
criteria; $(iv)$ It is much faster than existing algorithms. $(v)$ A geometric
interpretation is adopted to theoretically justify the relaxation of the
original combinatorial problem to continuous optimization problem.

##  3 The geometric point of view

In this section, we assume the information matrix at a single point $x$ is of
the form $M_{x}=xx^{T}$. Also, we focus on D-criterion and assume all
parameters are of interest, i.e. $p=0$ and $g(\theta)=\theta$. Note a sample,
say $S$, could be represented by a vector
$\delta=(\delta_{1},\delta_{2},...,\delta_{N})$, where $\delta_{i}=1$ if $i\in
S$ and $\delta_{i}=0$ otherwise. As a result, we have
$M(S)=\sum_{i=1}^{N}\delta_{i}x_{i}x_{i}^{T}$ and problem (2) reduces to

$\max_{\delta}\left|\sum_{i=1}^{N}\frac{\delta_{i}}{n}x_{i}x_{i}^{T}\right|,$
(3)

subject to $\delta\in\\{0,1\\}^{N}$ and $\sum_{i=1}^{N}\delta_{i}=n$. To
bypass the computational bottleneck of (3), it is natural to replace
$\delta_{i}/n$ by a continuous weight $w_{i}$ as follows.

$\max_{w\in\Omega}log|M(w)|.$ (P1)

where $M(w)=\sum^{N}_{i=1}w_{i}x_{i}x_{i}^{T}$ and
$\Omega=\\{w=\\{w_{1},w_{2},...,w_{N}\\}:w_{i}\geq 0,1\leq i\leq
N,\sum^{N}_{i=1}w_{i}=1\\}$. In approximate design theory, a vector
$w\in\Omega$ is called a measure, $M(w)$ is its information matrix, and the
solution of (P1) is called a D-optimal measure. However, (P1) is not a good
approximation of (3) since the number of points with positive weights (related
to $k$) is typically far less than $n$, i.e. the pre-specified number of
points to be selected. As shown in Figure 1, only two points are selected by
(P1) while $n=10$ may be needed.

Figure 1: The data consists of $100$ points generated from a two dimensional
normal distribution. The two solid points are the supporting points of the
solution of (P1) and the ten circled points are the supporting points of the
solution of (P2) and (P1-$\epsilon$). The dashed ellipsoid is the solution of
(D1) and the solid ellipsoid is the solution of (D2).

To address this issue, we shall first review on the geometric view of this
D-optimal design problem, and then explain our proposal in an intuitive way.
By Section 3 of Vandenberghe et al. (1998), the dual problem of (P1) is given
by

$\min_{W\in{\mathbb{R}}^{k\times k}}-log|W|,$ (D1)

subject to $W>0$ and $x_{i}^{T}Wx_{i}\leq k$, $1\leq i\leq N$. Note that the
set ${\cal E}(W)=\\{x|x^{T}Wx\leq k\\}$ features an ellipsoid centered at
origin and the region within it. Geometrically, (D1) tries to find the
smallest volume ellipsoid among all ellipsoids which contain all points in
${\cal X}$ and are centered at the origin. See the dashed ellipsoid in Figure
1 for ${\cal E}(W^{*})$, with $W^{*}$ being the solution of (D1).

Based on the duality between (P1) and (D1), we have $log|M(w)|\leq-log|W|$
with the equality holds if and only if $W^{-1}=M(w)$. As a result, we have the
identity $W^{*}=M(w^{*})^{-1}$, where $w^{*}=\\{w_{i}^{*},1\leq i\leq N\\}$ is
the solution of (P1). By Kiefer’s general equivalence theorem for
D-optimality, we have $w^{*}_{j}>0$ only for the points satisfying
$x_{j}^{T}W^{*}x_{j}=k$. Geometrically, it means only the points on the
boundary of ${\cal E}(W^{*})$ could be the support point of a D-optimal
measure. See the two solid circled points in the figure.

Note that the number of points on the boundary of ${\cal E}(W^{*})$ is
typically comparable to $k$ and hence far less than $n$. For the purpose of
selecting $n$ points from ${\cal X}$, the geometric representation of (D1)
inspires us to shrink the ellipsoid ${\cal E}(W^{*})$ until $n$ points emerge
on or outside the boundary of the new ellipsoid, as represented by the solid
ellipsoid in Figure 1. The collection of these points will be a reasonable,
though not necessarily optimal, solution for (3). See the $n=10$ circled
points outside the solid ellipsoid in the figure. To be more specific, we
define the soft-margin minimum volume ellipsoid problem as

$\min_{W\in{\mathbb{R}}^{p\times
p},\xi_{i}\in{\mathbb{R}}}-log|W|+\lambda\sum^{N}_{i=1}\xi_{i},$ (D2)

subject to $W>0$, $x_{i}^{T}Wx_{i}\leq k+\xi_{i}$ and $\xi_{i}\geq 0$, $1\leq
i\leq N$, with $\lambda$ being the penalty parameter. The introduction of the
slack variable $\xi_{i}$ allows some points to fall outside the ellipsoid.
Define the Lagrange
$L=-log|W|+\lambda\sum^{N}_{i=1}\xi_{i}-\sum^{N}_{i=1}\alpha_{i}(k+\xi_{i}-x_{i}^{T}Wx_{i})-\sum^{N}_{i=1}\mu_{i}\xi_{i}$,
from which we can derive the following Karush-Kuhn-Tucker (KKT) conditions.

$\displaystyle\frac{\partial L}{\partial
W}=-W^{-1}+\sum^{N}_{i=1}\alpha_{i}x_{i}x_{i}^{T}$ $\displaystyle=$
$\displaystyle 0,$ (4) $\displaystyle\frac{\partial
L}{\partial\xi_{i}}=\lambda-\alpha_{i}-\mu_{i}$ $\displaystyle=$
$\displaystyle 0,$ (5) $\displaystyle\alpha_{i},\mu_{i}$ $\displaystyle\geq$
$\displaystyle 0,1\leq i\leq N,$ (6)
$\displaystyle\mu_{i}\xi_{i}=\alpha_{i}(k+\xi_{i}-x_{i}^{T}Wx_{i})$
$\displaystyle=$ $\displaystyle 0,1\leq i\leq N.$ (7)

Plugging (4) and (5) back to the objective function of (D2), we derive its
dual problem as

$\max_{\alpha_{i}\in{\mathbb{R}},1\leq i\leq
N}log\left|\sum^{N}_{i=1}\alpha_{i}x_{i}x_{i}^{T}\right|+k(1-\sum\alpha_{i}),$
(P2)

subject to $0\leq\alpha_{i}\leq\lambda$. Here we utilize the equation

$\displaystyle\sum^{N}_{i=1}\alpha_{i}x_{i}^{T}\left(\sum^{N}_{i=1}\alpha_{i}x_{i}x_{i}^{T}\right)^{-1}x_{i}$
$\displaystyle=$ $\displaystyle k.$

As a result, the maximum of the objective function in (P2) shall be equal to
the minimum of the objective function in (D2). For the solutions of (P2) and
(D2), say $\alpha$ and $W$, we shall have the identity $W=M(\alpha)^{-1}$,
where $M(\alpha)=\sum^{N}_{i=1}\alpha_{i}x_{i}x_{i}^{T}$.

To this end, we can characterize the solution of (P2) in Theorem 1 below in
view of (6), (7) and the fact that $\xi_{i}>0$ if and only if
$x_{i}^{T}Wx_{i}>k$.

######  Theorem 1.

Suppose $\alpha=\\{\alpha_{1},...,\alpha_{N}\\}$ is the solution of (P2).
Then, for $1\leq i\leq N$, we have $(i)$ $\alpha_{i}=0$ if
$x_{i}^{T}M(\alpha)^{-1}x_{i}<k$. $(ii)$ $\alpha_{i}=\lambda$ if
$x_{i}^{T}M(\alpha)^{-1}x_{i}>k$. $(iii)$ $0\leq\alpha_{i}\leq\lambda$ if
$x_{i}^{T}M(\alpha)^{-1}x_{i}=k$.

It is worth noticing that the form of (P2) appear close to the objective
function in the literature of finite population sampling and constraint space
optimal design of experiment. See Wynn (1977, 1982), Fedorov (1989), Müller
and Pázman (1998), Sahm and Schwabe (2001), Pronzato (2004, 2006), Uciński
(2015) for relevant literature. One particular version of such problems can be
written as follows.

$\max_{w\in\Omega_{\epsilon}}log|M(w)|,$ (P1-$\epsilon$)

where $0<\epsilon<1$ and
$\Omega_{\epsilon}=\\{(w_{1},w_{2},...,w_{N})\in\Omega:0\leq
w_{i}\leq\epsilon,1\leq i\leq N\\}$. Here, we shall establish the connection
between (P2) and (P1-$\epsilon$), and then provide further insights,
improvements and extensions that are vital for our current context of BMAL.

For the connection, we shall prove that the solution path of (P2) is covered
by (P1-$\epsilon$) as illustrated as in Theorem 2. Relatively speaking, the
objective function in (P1-$\epsilon$) is simpler than that in (P2). By
choosing a proper value of $\epsilon$, e.g. $\epsilon=1/n$, we can guarantee
at least $n$ positive weights in the solution of (P1-$\epsilon$) and the
sample could be chosen according to the $n$ largest weights. Since our goal
here is to select the important points, it is not the absolute but the
relative values of $\alpha_{i}$ in (P2) that matters.

######  Theorem 2.

For given solution of (P2), say $\alpha=\\{\alpha_{1},...,\alpha_{N}\\}$,
under the penalty $\lambda$. Let $\epsilon=\lambda/\sum^{N}_{i=1}\alpha_{i}$,
the solution of (P1-$\epsilon$) will be given by
$w=\alpha/\sum^{N}_{i=1}\alpha_{i}$.

###### Proof.

Let $w_{i}=\alpha_{i}/\sum^{N}_{i=1}\alpha_{i}$, the target function of (P2)
becomes
$log|M(w)|+k(1-\sum^{N}_{i=1}\alpha_{i}+\log(\sum^{N}_{i=1}\alpha_{i}))$. By
fixing the second term at the solution of (P1-$\epsilon$), the maximization of
the first term can be achieved by solving (P1-$\epsilon$) under the constraint
of $0\leq
w_{i}=\alpha_{i}/\sum^{N}_{i=1}\alpha_{i}\leq\lambda/\sum^{N}_{i=1}\alpha_{i}=\epsilon$.
The constraint of $\sum^{N}_{i=1}w_{i}=1$ is inherently guaranteed by the
construction of $w_{i}$. ∎

To this end, we have justified approximating problem (3) by (P1-$\epsilon$).
While this latter problem has been studied in literature, there are several
questions and issues as listed below that are not resolved for the current
application. It is main purpose of this paper to address these issues.

  1. 1.

The upper bound $\epsilon$ serves to enforce at least $1/\epsilon$ points to
have positive weights, but it does not necessarily means the points such
selected is a good solution for the original problem (2). Do we have a
theoretical justification for (P1-$\epsilon$)?

  2. 2.

There seems to be no control for the maximum number of points with positive
weights. What if there are too many of them as compared to $n$?

  3. 3.

So far, the algorithms proposed in literature are statistically efficient, but
not computational efficient. Do we have a much faster algorithm that still
achieves the high statistical efficiency?

  4. 4.

So far, most work typically focus on a particular model and work on a
particular criterion, most frequently D-criterion. Can we propose an algorithm
that is readily applicable to any arbitrary model with interest in any
function ($g$) of the parameters under any $\Phi_{p}$-criterion?

The first issue has been largely addressed by the primal-dual argument coupled
with the geometric interpretation of the problem as illustrated by Figure 1.
Through the empirical results, we shall observe that our proposed algorithm
would achieve close 100% efficiency for all examples. The second issue is
addressed by Corollary 1, which indicates that the number of points with
positive weights will be very close to $n$ when we set $\epsilon=1/n$. The
solutions for the last two issues will be provided in Section 4.

######  Corollary 1.

Suppose $w=\\{w_{1},...,w_{N}\\}$ is the solution of (P1-$\epsilon$). Then
there exists a constant, say $c$, such that for $1\leq i\leq N$, $(i)$
$w_{i}=0$ if $x_{i}^{T}M(w)^{-1}x_{i}<c$. $(ii)$ $w_{i}=\epsilon$ if
$x_{i}^{T}M(w)^{-1}x_{i}>c$. $(iii)$ $0\leq w_{i}\leq\epsilon$ if
$x_{i}^{T}M(w)^{-1}x_{i}=c$.

This corollary can be obtained by investigating the linear proportionality
between the solutions of (P2) and (P1-$\epsilon$). It characterized the
solution of (P1-$\epsilon$) by the trichotomy structure as in Theorem 1, which
is the key for the development of efficient algorithms. Note that the proof of
Corollary 1 is, to some extend, inspired by the KKT condition arguments
adopted in supporter vector machine. As shown in Figure 2, only few points
will have the weight strictly between $0$ and $\epsilon$. In fact, these are
the points fall exactly on the small ellipsoid in Figure 1. Similarly, the
points with $w_{i}=\epsilon$ is exactly the points outside the small ellipsoid
in Figure 1. This addresses the second issue as mentioned earlier. The short
answer is that the total number of points with positive weight will be only
slightly larger than $n$, if not equal.

Figure 2: The plot of leverage score $x_{i}^{T}M(w)^{-1}x_{i}$ versus the
weight $w_{i}$ when $w$ is optimal for (P1-$\epsilon$).

To address the third issue, Theorem 3 provide a Kiefer’s type of equivalence
theorem for the optimality conditions of (P1-$\epsilon$), which sheds light on
the construction of Algorithm 1 below. Once the solution of (P1-$\epsilon$) is
derived, we shall project it to the solution space of (3) as a reasonable
solution for (3). As a result, our sample consists of the $n$ points with the
largest weights in the solution of (P1-$\epsilon$).

######  Theorem 3.

(i) A measure $w$ solves (P1-$\epsilon$) if and only if

$\displaystyle\max_{w^{\prime}\in\Omega_{\epsilon}}Tr(M(w)^{-1}M(w^{\prime}))$
$\displaystyle=$ $\displaystyle k$ (8)

with the maximum achieved by $w^{\prime}=w$.
(ii) If a measure $w$ is a solution of (P1-$\epsilon$), we have $w_{i}\leq
w_{j}$ whenever $x_{i}^{T}M(w)^{-1}x_{i}<x_{j}^{T}M(w)^{-1}x_{j}$.

###### Proof.

Let $L(w)=log|M(w)|$. By the concavity of $L(w)$, $w$ maximises $L(w)$ if and
only if

$\displaystyle 0\geq\left.\frac{\partial L((1-\alpha)w+\alpha
w^{\prime})}{\partial\alpha}\right|_{\alpha=0}=Tr(M(w)^{-1}M(w^{\prime}))-k$
(9)

for any measure $w^{\prime}\in\Omega_{\epsilon}$. That indicates

$\displaystyle\max_{w^{\prime}\in\Omega_{\epsilon}}Tr(M(w)^{-1}M(w^{\prime}))$
$\displaystyle\leq$ $\displaystyle k$ (10)

Since $Tr(M(w)^{-1}M(w))=k$ and

$\displaystyle
Tr(M(w)^{-1}M(w^{\prime}))=\sum^{N}_{i=1}w^{\prime}_{i}x_{i}^{T}M(w)^{-1}x_{i}$
(11)

is linear in $w^{\prime}$, hence we should also have

$\displaystyle\max_{w^{\prime}\in\Omega_{\epsilon}}Tr(M(w)^{-1}M(w^{\prime}))$
$\displaystyle\geq$ $\displaystyle k$ (12)

Then (8) is a direct result of (10) and (12). For part $(ii)$, suppose
$w_{i}>w_{j}$, (9) shall be violated by taking $w^{\prime}_{i}=w_{j}$,
$w^{\prime}_{j}=w_{i}$ and the rest components of $w^{\prime}$ being equal to
that of $w$. ∎

To introduce our algorithm for (P1-$\epsilon$), we define a measure, say
$\tilde{w}$, to be the steepest gradient (SG) measure of $w$ if it achieves
the maximum in (8). For a given measure, its SG measure represents the
direction of fastest improvement on the objective function. Particularly, a
measure is the SG measure of itself if it is the optimal solution of
(P1-$\epsilon$). Let $S_{n}(w)\subset\\{1,2,...,N\\}$ be the index of $n$
points with the largest value of $x_{i}^{T}M(w)^{-1}x_{i}$. By (11), we shall
be able to explicitly construct the SG measure of $w$ in the following way. If
$1/\epsilon$ is an integer, let $\tilde{w}_{i}=\epsilon$ for $i\in
S_{1/\epsilon}(w)$ and $\tilde{w}_{i}=0$ for $i\not\in S_{1/\epsilon}(w)$. In
general, set $\tilde{w}_{i}=\epsilon$ for $i\in S_{[1/\epsilon]}(w)$,
$\tilde{w}_{i}=1-\epsilon[1/\epsilon]$ for $i\in
S_{[1/\epsilon]+1}(w)-S_{[1/\epsilon]}(w)$ and $\tilde{w}_{i}=0$ for $i\not\in
S_{[1/\epsilon]+1}(w)$. Here, $[1/\epsilon]$ is the largest integer that is
smaller than $1/\epsilon$. The value $x_{i}^{T}M(w)^{-1}x_{i}$ is called the
leverage score in literature. It is the variance of predicted value of the
response at $x_{i}$ in a linear model. Hence, a SG measure tries to collect
the most unpredictable points from the dataset. Now we propose the following
algorithm for (P1-$\epsilon$).

Algorithm 1 Step 1: Initialization. Give the initial solution $w^{0}$ such
that $M(w^{0})>0$. Step 2: Update. At iteration $t\geq 1$, calculate the SG
measure of $w^{t-1}$, namely
$\tilde{w}^{t-1}=\\{\tilde{w}_{1}^{t-1},...,\tilde{w}_{N}^{t-1}\\}$. Let
$T^{t}_{1}=\\{i:w^{t-1}_{i}=\tilde{w}^{t-1}_{i}=\epsilon\\}$ and
$T^{t}_{2}=\\{i:w^{t-1}_{i}=\tilde{w}^{t-1}_{i}=0\\}$. Maximize $log|M(w)|$
subject to the constraints that $w_{i}=\epsilon$ for $i\in T^{t}_{1}$,
$w_{i}=0$ for $i\in T^{t}_{2}$ and $w\in\Omega_{\epsilon}$. Denote the
solution by $w^{t}$. Step 3: Check the optimality based on Theorem 3. If
$Tr(M(w^{t})^{-1}M(\tilde{w}^{t}))<k(1+v)$, with $v$ being a sufficiently
small positive constant, stop and return the value of $w^{t}$. Otherwise,
repeat Step 2 until the convergence criterion is met.

The computational complexity of deriving the SG measure is $O(k^{3}+Nk^{2})$
and the computational complexity of maximizing $log|M(w)|$ step 2 is
$O((m+n)k^{2}+k^{3})$, where $m$ is the number of elements neither in
$T^{t}_{1}$ nor in $T^{t}_{2}$. Hence it is important to choose a proper
initial solution so that $m$ is small in the beginning. This can be achieved
by selecting the set $S$ of $n$ points according to the largest leverage score
$x_{i}^{T}(\sum^{N}_{i=1}x_{i}x_{i}^{T})^{-1}x_{i}$ and let $w_{i}^{0}=1/n$
for $i\in S$. It is justified by the fact that the leverage score is
proportional to the derivative of $log|M(w)|$ with respective to $w_{i}$
evaluated at the uniform measure $w_{i}=1/N$, $1\leq i\leq N$. Hence, the
initial measure such selected is essentially the SG measure of the uniform
measure. As a result, we have $m\leq 2n$ and typically we have $m=o(n)$
especially in later iterations. To sum up, the computational complexity of
Algorithm 1 is $O(k^{3}+Nk^{2})$. For fixed dimension $k$, it is simply
$O(N)$. Another idea to expedite the algorithm is to update the measure by a
linear combination of the measure and its SG measure with the weight carefully
chosen. These ideas will be implemented in Algorithm 3 for the general
$\Phi_{p}$ optimality.

##  4 The $\Phi_{p}$ optimality

In this section, we try to address the last two issues as proposed in Section
3, that is to derive an efficient (both statistically and computationally)
algorithm to derive near-optimal solution of (2) in a unified framework for
any model, any function of parameters and any design criterion. Our framework
shall permit the evaluation of the statistical efficiency of the derived
samples as compared to the actual optimal one, even thought the latter may not
be known. The theoretical results regarding the global convergence is also
obtained. The speed of the algorithm will be illustrated by the examples in
the next section.

Inspired by the geometric justification of (P1-$\epsilon$) for D-criterion, we
propose to approximate (2) by the following continuous optimization problem.

$\displaystyle\min_{w\in\Omega_{\epsilon}}\Phi_{p}(\Sigma_{g}(w)),$ (13)

for $0\leq p<\infty$, where $\Phi_{p}$ is defined in (1) and

$\displaystyle\Sigma_{g}(w)$ $\displaystyle=$ $\displaystyle\frac{\partial
g(\theta)}{\partial\theta^{T}}M(w)^{-1}\left(\frac{\partial
g(\theta)}{\partial\theta^{T}}\right)^{T},$ $\displaystyle M(w)$
$\displaystyle=$ $\displaystyle\sum^{N}_{i=1}w_{i}M_{x_{i}}.$

Here, the information matrix $M_{x_{i}}$ from a single point $x_{i}$ does not
has to be of the particular form $x_{i}x_{i}^{T}$. When we take $g$ to be the
identity function and $p=0$, we have $\Phi_{p}(\Sigma_{g}(w))=|M(w)|^{-1/k}$.
Hence, (13) covers (P1-$\epsilon$) as a special case. For a sample $S$ of size
$n$, let $w(S)\in\Omega$ be the measure assigning weight $1/n$ to each element
of $S$. Obviously, we have
$\Phi_{p}(\Sigma_{g}(w(S)))=n\Phi_{p}(\Sigma_{g}(S))$. Note that
$w(S)\in\Omega_{\epsilon}$ for all $S\subset\\{1,2,...,N\\}$ if and only if
$\epsilon\geq 1/n$. In this case, the domain of (13) covers the domain of (1)
and hence we can evaluate the (lower bound of) efficiency of $S$

$\displaystyle{\cal E}_{p}(w)$ $\displaystyle:=$
$\displaystyle\frac{\Phi_{p}(\Sigma_{g}(w^{*}))}{\Phi_{p}(\Sigma_{g}(w(S)))},$
(14)

where $w^{*}$ is the solution of (13). As in Section 3, we shall sample the
$n$ points with the largest weights in $w^{*}$ as an approximate solution for
the minimization of (1). To make sure $w^{*}$ has at least $n$ positive
weights, we need the inequality $\epsilon\leq 1/n$. As a result, we would
suggest to simply set $\epsilon=1/n$ in general.

Now we shall generalize the results of Section 3 for the $\Phi_{p}$ criterion
and prove the global convergence of the generalized algorithm. The following
result plays a crucial role in the development of Algorithms 2 to be proposed.

######  Lemma 1.

Let $\Omega_{\epsilon}^{+}=\\{w\in\Omega_{\epsilon}:M(w)>0\\}$. For $0\leq
p<\infty$, $\Phi_{p}(\Sigma_{g}(w))$ is convex on the domain of
$\Omega_{\epsilon}^{+}$.

The proof follows by the same arguments as in Lemma 1 of Yang et al. (2013)
under the classical framework of approximate design theory, except that
$\xi_{0}$ therein shall be deleted. In Yang et al. (2013), the inclusion of
$\xi_{0}$ as the first stage design is necessary for the proof of the
convexity since they worked on the domain $\Omega$ for the second stage
design. Here, we bypass the need for the first stage design by shrinking the
domain from $\Omega_{\epsilon}$ to $\Omega_{\epsilon}^{+}$. Such adjustment is
crucial for the establishment of the theoretical convergence of our proposed
algorithms. By Lemma 1, a measure $w$ is a solution of (13) if and only if the
Fréchet direvitive of the objective function is nonnegative, namely

$\displaystyle\eta(w^{\prime},w):=\left.\frac{\partial\Phi_{p}(\Sigma_{g}((1-\alpha)w+\alpha
w^{\prime}))}{\partial\alpha}\right|_{\alpha=0}$ $\displaystyle\geq$
$\displaystyle 0.$ (15)

By direct calculation we have

$\displaystyle-\eta(w^{\prime},w)$ $\displaystyle=$
$\displaystyle\sum^{N}_{i=1}(w^{\prime}_{i}-w_{i})\phi(x_{i},w)$
$\displaystyle=$
$\displaystyle\sum^{N}_{i=1}w^{\prime}_{i}\phi_{p}(x_{i},w)-\Phi_{p}(\Sigma_{g}(w)),$
$\displaystyle\phi_{p}(x,w)$ $\displaystyle=$ $\displaystyle
q^{-1/p}\left(Tr\left[\Sigma_{g}(w)\right]^{p}\right)^{1/p-1}$
$\displaystyle\times Tr\left(\Sigma_{g}(w)^{p-1}\frac{\partial
g(\theta)}{\partial\theta^{T}}M(w)^{-1}M_{x}M(w)^{-1}\left(\frac{\partial
g(\theta)}{\partial\theta^{T}}\right)^{T}\right).$

The value of $\phi_{p}(x,w)$ can be considered as the $\Phi_{p}$ version of
the leverage. By using the same arguments as in Theorem 3 with (9) and (11)
therein replaced by (15) and (4), we have

######  Theorem 4.

(i) A measure $w$ solves (13) if and only if

$\displaystyle\max_{w^{\prime}\in\Omega_{\epsilon}}\sum^{N}_{i=1}w^{\prime}_{i}\phi_{p}(x_{i},w)$
$\displaystyle=$ $\displaystyle\Phi_{p}(\Sigma_{g}(w))$ (17)

with the maximum achieved by $w^{\prime}=w$.
(ii) If a measure $w$ is a solution of (P1-$\epsilon$), we have $w_{i}\geq
w_{j}$ whenever $\phi_{p}(x_{i},w)>\phi_{p}(x_{j},w)$.

Theorem 4 generalizes Theorem 3 in characterising the solution of (13).
Moreover, the trichotomy phenomena for (P1-$\epsilon$) also exists for (13)
with the details given by Corollary 2.

######  Corollary 2.

Suppose $w=\\{w_{1},...,w_{N}\\}$ is the solution of (13). Then there exists a
constant, say $c$, such that for $1\leq i\leq N$, $(i)$ $w_{i}=0$ if
$\phi_{p}(x_{i},w)<c$. $(ii)$ $w_{i}=\epsilon$ if $\phi_{p}(x_{i},w)>c$.
$(iii)$ $0\leq w_{i}\leq\epsilon$ if $\phi_{p}(x_{i},w)=c$.

To develop the algorithm for (13), we shall generalize the definition of SG
measure $\tilde{w}$ of $w$ in Section 3 to be the measure satisfying the
following equation.

$\displaystyle\eta(\tilde{w},w)$ $\displaystyle=$
$\displaystyle\min_{w^{\prime}\in\Omega_{\epsilon}}\eta(w^{\prime},w).$ (18)

The convergence of an algorithm typically relies on the convexity of the
objective function. Inspired by Lemma 1, we shall consider the measure which
achieves the minimum in (18) among $\Omega_{\epsilon}^{+}$ instead of
$\Omega_{\epsilon}$. Such a measure is said to be the positive steepest
gradient (PSG) measure of $w$. Since $n>>k$, the SG measure is almost always
the PSG measure in practice. For data generated by continuous distribution,
even $n=k$ will insure the positiveness of the SG measure with probability
$1$. We use the PSG measure in Algorithm 2, only for the need in the proof of
convergence in Theorem 5. In practice, we can directly adopt the SG measure
while its positiveness is checked during the process. The construction of SG
measure for the $\Phi_{p}$ criterion is similar to that for the D criterion,
only by replacing $x_{i}^{T}M(w)^{-1}x_{i}$ therein by $\phi_{p}(x_{i},w)$.
Particularly, when $1/\epsilon$ is an integer, the SG measure of $w$ is the
measure with its weight equal to $\epsilon$ for $1/\epsilon$ points with
largest values of $\phi_{p}(x_{i},w)$.

Algorithm 2 Step 1: Initialization. Give the initial solution $w^{0}$ such
that $M(w^{0})>0$. Step 2: Update. At iteration $t\geq 1$, calculate the PSG
measure of $w^{t-1}$, namely
$\tilde{w}^{t-1}=\\{\tilde{w}_{1}^{t-1},...,\tilde{w}_{N}^{t-1}\\}$. Minimize
$\Phi_{p}(\Sigma_{g}(w))$ subject to the constraints that $w_{i}=\epsilon$ for
$i\in T^{t}_{1}$, $w_{i}=0$ for $i\in T^{t}_{2}$ and $w\in\Omega_{\epsilon}$,
where $T^{t}_{1}$ and $T^{t}_{2}$ are defined as in Algorithm1. Denote the
solution by $w^{t}$. Step 3: Check the optimality based on Theorem 4. If
$\sum^{N}_{i=1}\tilde{w}_{i}^{t}\phi_{p}(x_{i},w)>(1-v)\Phi_{p}(\Sigma_{g}(w^{t}))$,
with $v$ being a sufficiently small positive constant, stop and return the
value of $w^{t}$. Otherwise, repeat Step 2 until the convergence criterion is
met.

The minimization in Step 2 can be carried out by any generic optimization
method, where the derivative of the objective function with respect to $w_{i}$
is $\phi_{p}(x_{i},w)$ and the domain is $\Omega_{\epsilon}$. Theorem 5 proves
the global convergence of Algorithm 2.

######  Theorem 5.

Let $w^{t}$ be the solution at iteration $t$ of Algorithm 2 and $\hat{w}$ be
the solution of (13), we have

$\displaystyle\lim_{t\rightarrow\infty}\Phi_{p}(\Sigma_{g}(w^{t}))$
$\displaystyle=$ $\displaystyle\Phi_{p}(\Sigma_{g}(\hat{w})).$ (19)

###### Proof.

Suppose (19) does not hold, we shall have a constant, say $\Delta>0$, such
that

$\displaystyle\Phi_{p}(\Sigma_{g}(w^{t}))-\Phi_{p}(\Sigma_{g}(\hat{w}))$
$\displaystyle\geq$ $\displaystyle\Delta,$ (20)

for all $t\geq 1$. By Lemma 1, we have

$\displaystyle\Phi_{p}(\Sigma_{g}((1-\alpha)w^{t}+\alpha\hat{w}))$
$\displaystyle\leq$
$\displaystyle(1-\alpha)\Phi_{p}(\Sigma_{g}(w^{t}))+\alpha\Phi_{p}(\Sigma_{g}(\hat{w})),$

for any $0\leq\alpha\leq 1$. By (15) and simple algebra and calculus, we have

$\displaystyle\eta(\hat{w},w^{t})$ $\displaystyle\leq$
$\displaystyle\Phi_{p}(\Sigma_{g}(\hat{w}))-\Phi_{p}(\Sigma_{g}(w^{t})),$ (21)

Let $\tilde{w}^{t}$ be the PSG measure of $w^{t}$. By (18), (20) and (21), we
have

$\displaystyle\eta(\tilde{w}^{t},w^{t})$ $\displaystyle\leq$
$\displaystyle\eta(\hat{w},w^{t})\leq-\Delta$ (22)

Denote

$\displaystyle\tau(\alpha,w^{\prime},w)$ $\displaystyle:=$
$\displaystyle\frac{\partial^{2}\Phi_{p}(\Sigma_{g}((1-\alpha)w+\alpha
w^{\prime}))}{\partial\alpha^{2}}$ (23)

By the same arguments for Theorem 3 of Yang et al. (2013), there exist a
constant $K$ such that $\sup\\{\tau(\alpha,\tilde{w}^{t},w^{t}):t\geq
1,0\leq\alpha\leq 1/2\\}=K<\infty$. Now we have the taylor’s expansion

$\displaystyle\Phi_{p}(\Sigma_{g}((1-\alpha)w^{t}+\alpha\tilde{w}^{t}))$
$\displaystyle=$
$\displaystyle\Phi_{p}(\Sigma_{g}(w^{t}))+\alpha\eta(\tilde{w}^{t},w^{t})+\alpha^{2}\tau(\alpha^{\prime},\tilde{w}^{t},w^{t})/2$
$\displaystyle\leq$
$\displaystyle\Phi_{p}(\Sigma_{g}(w^{t}))-\Delta\alpha+K\alpha^{2}/2$

where $0\leq\alpha^{\prime}\leq\alpha\leq 1/2$. By taking $\alpha=\Delta/K\leq
1/2$, we have

$\displaystyle\Phi_{p}(\Sigma_{g}((1-\alpha)w^{t}+\alpha\tilde{w}^{t}))$
$\displaystyle\leq$
$\displaystyle\Phi_{p}(\Sigma_{g}(w^{t}))-\Delta^{2}/(2K).$ (24)

However, the updating step of Algorithm 2 implies

$\displaystyle\Phi_{p}(\Sigma_{g}(w^{t+1}))\leq\Phi_{p}(\Sigma_{g}((1-\alpha)w^{t}+\alpha\tilde{w}^{t})).$
(25)

By (24) and (25), we have
$\Phi_{p}(\Sigma_{g}(w^{t}))=\Phi_{p}(\Sigma_{g}(w^{0}))-t\Delta^{2}/(2K)\rightarrow-\infty$
as $t$ goes to infinity, which contradict with the fact that
$\Phi_{p}(\Sigma_{g}(w))\geq 0$ for all $w\in\Omega$. ∎

Inspired by (24) and (25), we propose to speed up the updating step as
follows. For $\alpha$ in the neighbourhood of zero, we have the approximation

$\displaystyle\Phi_{p}(\Sigma_{g}((1-\alpha)w+\alpha\tilde{w}))$
$\displaystyle\approx$
$\displaystyle\Phi_{p}(\Sigma_{g}(w))+\alpha\eta(\tilde{w},w)$
$\displaystyle+\alpha^{2}\tau(0,\tilde{w},w)/2.$

Here, the function $\tau$ as defined in (23) is simply the second order
derivative. By the convexity of $\Phi_{p}$, we have $\tau(0,\tilde{w},w)\geq
0$ with equality achieved by the optimal measure. Hence, the right hand side
of (4) is minimized at $\alpha^{*}=-\eta(\tilde{w},w)/\tau(0,\tilde{w},w)$ and
we shall update $w$ by the new measure as
$(1-\alpha^{*})w+\alpha^{*}\tilde{w}$. To avoid the singularity issue, we can
add a very small amount to the denominator of $\alpha^{*}$. Since
$\eta(\tilde{w}^{t-1},w^{t-1})\leq 0$ with the equality achieved by the
optimal measure, we shall have $\alpha^{*}\geq 0$ in this case. To avoid the
negative case due to numerical error, we can truncate $\alpha^{*}$ by $0$ from
the left. Meanwhile, in order for the approximation to be valid, we should
truncate $\alpha^{*}$ in the right by a small value, say $0<r<1$. In practice,
$r=1/4$ works well. To summarize, we shall take

$\displaystyle\alpha^{*}$ $\displaystyle=$
$\displaystyle\min\left[r,\max\left(0,\frac{-\eta(\tilde{w},w)}{\tau(0,\tilde{w},w)+u}\right)\right]$
(27)

with $u>0$ being an arbitrarily small constant. Due to its simplicity, this
updating procedure $(1-\alpha^{*})w+\alpha^{*}\tilde{w}$ shall be much faster
than step 2 of Algorithm 2. Obviously, it can quickly boost up the speed of
updating process in the early stages. However, it will be converge slowly in
the later stages when precision requirement is high, that is when $v$ in step
3 of Algorithm 2 is small. The slow down phenomena could be explained by the
sparsity property of the optimal solution as characterized by Corollary 2.
Hence, we shall use the updating step of Algorithm 2 in the later stages. To
take advantages of both sides, we propose the following hybrid algorithm.

Algorithm 3 (Hybrid Algorithm) Step 1: Initialization. Take initial solution
$w^{0}$ to be the SG measure of the uniform measure. Step 2: Boost. At
iteration $t\geq 1$, calculate the SG measure of $w^{t-1}$, namely
$\tilde{w}^{t-1}$. Let
$w^{t}=(1-\alpha^{t})w^{t-1}+\alpha^{t}\tilde{w}^{t-1}$, where $\alpha^{t}$ is
calculated by (27). Step 3: Check. If
$\sum^{N}_{i=1}\tilde{w}_{i}^{t}\phi_{p}(x_{i},w)>(1-v_{0})\Phi_{p}(\Sigma_{g}(w^{t}))$,
with $v_{0}$ being a sufficiently small positive constant. If yes, go to step
4. Otherwise, repeat Step 2 until the criterion is met. Step 4: Optimization.
Run Algorithm 2 with the solution from Step 3 as the initial solution and the
precision parameter $v>0$.

In general, we take $v<v_{0}$ so that step 4 will further improve on the
measure provided by step 3. The flexibility of the algorithm also allows us to
take $v>v_{0}$ so that step 4 will be skipped. This can be useful when the
accuracy requirement is not super high. In practice, steps 1-3 of Algorithm 3
is very time effective for $v=0.001$ or larger. Since Algorithm 2 is embedded
into the hybrid algorithm, we shall still have the global convergence of the
latter.

######  Corollary 3.

Let $w^{t}$ be the solution at iteration $t$ of the hybrid algorithm and
$\hat{w}$ be the solution of (13), we have

$\displaystyle\lim_{t\rightarrow\infty}\Phi_{p}(\Sigma_{g}(w^{t}))$
$\displaystyle=$ $\displaystyle\Phi_{p}(\Sigma_{g}(\hat{w})).$ (28)

##  5 Examples

###  5.1 Logistic regression

Let $Y_{i}\in\\{0,1\\}$ and $z_{i}\in{\mathbb{R}}^{k}$, $1\leq i\leq N$, be
the binary response and the exploratory variable, respectively. In logistic
regression, the probability of success, $p_{i}=P(Y_{i}=1)$, is modelled by

$\displaystyle log(p_{i}/(1-p_{i}))=z_{i}^{T}\beta.$ (29)

The information matrix of a sample $S\subset\\{1,2,...,N\\}$ for $\beta$ is
$M(S)=\sum_{i\in S}x_{i}x_{i}^{T}$, where $x_{i}=\sqrt{p_{i}(1-p_{i})}z_{i}$.
Correspondingly, we have $M(w)=\sum_{i=1}^{N}w_{i}x_{i}x_{i}^{T}$ for a
measure $w\in\Omega_{\epsilon}$. This model coupled with the objective
function (3) is adopted by Ouyang et al. (2016a) when they tried to solve the
validation sampling problem for error-prone medical records. They proposed the
backward algorithm, that is to begin with $S^{0}=\\{1,2,...,N\\}$ and delete
one point each time until the sample size reduces to $n$. At iteration $t$,
they delete the point from $S^{t}$ with smallest value of
$x_{i}^{T}M(S^{t})x_{i}$, namely the point most predictable by the model
fitted to $S^{t}$. This is a type of greedy algorithm, which tries to minimize
the loss of information at each step of deletion. Besides the proposed
algorithms in this paper, we will also include the exchange algorithm in the
comparison. That is to iteratively delete the point with the smallest leverage
in the sample and add the point with the largest leverage outside the sample
in each updating step. We have argued that (3) can be approximated by (14)
with $p=0$ and $g$ being the identity function. The latter can be solved by
the hybrid algorithm proposed in this paper. We shall evaluate the (lower
bound of) D-efficiency of a candidate sample by

$\displaystyle{\cal E}_{0}(w(S))$ $\displaystyle=$
$\displaystyle\frac{|M(w(S))|^{1/k}}{|M(w^{*})|^{1/k}},$ (30)

where $w^{*}$ is the solution of (14).

An experiment to compare three algorithms in terms of both time and
statistical efficiency is carried out in the R environment on a desktop with
3.4 GHz cpu. The study is based on simulated data, where each point consists
of the constant term and a random vector from $10$ dimensional standard normal
distribution. The results are given in Table 1. While all the algorithms are
able to produce near-optimal samples, the time efficiency turns out to be
significantly different. The exchange algorithm only takes a fraction of the
time needed by the backward algorithm and the hybrid algorithm further takes a
fraction of what is needed by the exchange algorithm. It is possible to tune
the exchange algorithm to run faster, but it unlikely the improved version
could outperform the hybrid algorithm. We can also notice in the table that
the D-efficiency of the backward algorithm is slightly higher than the other
two by a tiny margin, that is because exchange and hybrid algorithms are set
to stop at $v=10^{-3}$ for $N=10,000$ and $v=10^{-2}$ for larger values of
$N$. If we set $v=10^{-6}$ for $N=10,000$, the hybrid algorithm produces a
sample with D-efficiency of $0.999995$ and it takes only $3.33$ seconds.

N | 10,000 | 100,000 | 1,000,000 | 10,000,000
---|---|---|---|---
Backward | 6.64 mins | 12.88 hours | NA | NA
0.9999939 | 0.9999895 | NA | NA
Exchange | 24.02 secs | 5.31 mins | 1.24 hours | 20.40 hours
0.9994489 | 0.9950084 | 0.9956839 | 0.9947914
Hybrid | 0.68 secs | 5.75 secs | 1.17 mins | 20.38 mins
0.9999827 | 0.999884 | 0.999649 | 0.9993768
Table 1: The time cost and D-efficiencies of the backward, exchange and hybrid
algorithms for selecting $n=1000$ points from populations of different sizes.
The dimension is $k=11$. Each point consists of the constant term and a random
vector from $10$ dimensional standard normal distribution.

While the backward algorithm works on the D-optimality and the case when all
the parameters of are interest, the hybrid algorithm on the other hand can
tackle with any $\Phi_{p}$-optimality regarding any smooth function of the
parameters. The time efficiency of the latter has been shown in Table 1. Now
we use it to investigate the influence of different design criteria on the
choice of samples. Here, we focus on the A and D criteria. For fixed
$N=10,000$ and various choices of $n$, we derive both the A-optimal and
D-optimal samples and evaluate their efficiencies against each other under the
opposite criteria. As shown by Table 2, an A-optimal sample does not
necessarily perform well under the D-criterion and vice versa. The discrepancy
between the two criteria is larger when the sample size to population size
ratio, i.e. $n/N$, is small. That is because less common points will be shared
by the two samples when there are relatively fewer points to be selected.

$n$ | 500 | 1000 | 3000 | 5000
---|---|---|---|---
A-efficiency | 0.6511 | 0.7064 | 0.8485 | 0.9439
D-effficiency | 0.8022 | 0.8287 | 0.8951 | 0.9484
Table 2: The A-efficiency of the D-optimal sample and the D-efficiency of the
A-optimal sample. The population is of fixed size 10,000. Each point consists
of the constant term and a random vector from $10$ dimensional standard normal
distribution.

Let $\lambda_{1},...\lambda_{k}$ be the eigenvalues of the inverse of the
information matrix $M^{-1}$. Recall that A-criterion aims to minimize
$\sum^{k}_{i=1}\lambda_{i}$ and D-criterion aims to minimize
$\Pi^{k}_{i=1}\lambda_{i}$. As pointed out by Ouyang et al. (2017), the
A-criterion treats the eigenvalues equally important while the D-criterion
allows certain eigenvalues to be large as long as some eigenvalues are small
enough. We can understand this by looking at the confidence region of the
parameters, which is in the form of
$\\{\beta:(\beta-\hat{\beta})^{T}M(\beta-\hat{\beta})\leq c\\}$ for some
constant $c$. It represents an ellipsoid, where the lengths of the axes are
proportional to the $\lambda_{i}$’s. In order for the volume of the ellipsoid
to be small, it is preferable to have a long but narrow shaped ellipsoid
rather than a round shaped ellipsoid. The former could be achieved by choosing
certain eigenvalues to be small enough while keeping others within some range.

###  5.2 Ordinal response: cumulative link model

Recently, Yang et al. (2017) studied the optimal designs for the cumulative
link model

$\displaystyle h(\gamma_{j})=\theta_{j}-z^{T}\beta,$ $\displaystyle
j=1,2,...,J-1,$ (31)

where $\gamma_{j}=P(Y\leq j)$ is the probability for the ordinal response
$Y\in\\{1,2,...,J\\}$ to be no larger than $j$ and $h$ is the link function.
They worked on the classical design problem where the weights, $w_{i}$’s, are
not bounded by $\epsilon$. When $h$ is taken to be the logit link function,
(31) is reduced to the proportional odds model. Perevozskaya et al. (2003)
studied optimal designs for this model when $z$ is univariate. For the general
model (31) with $z\in{\mathbb{R}}^{d}$, Yang et al. (2017) revealed an
interesting fact that the minimum number of experimental settings is $d+1$,
which is strictly less than the number of parameters $d+J-1$. This fact is of
more interest for a classical design problem than the sampling problem here,
since we are enforcing the constraint $w\leq\epsilon$ for the weights of the
designs points and a large number of supporting points is desired. By direct
calculations, the information matrix for the parameter
$\theta=(\beta^{T},\theta_{1},\theta_{2},...,\theta_{J-1})^{T}$ at the
covariate $z$ is

$\displaystyle M_{z}$ $\displaystyle=$
$\displaystyle\sum^{J}_{j=1}\frac{1}{\pi_{j}}x_{j}x_{j}^{T},$ $\displaystyle
x_{j}$ $\displaystyle=$
$\displaystyle\left(\begin{tabular}[]{c}$[f_{j}(z)-f_{j-1}(z)]z$\\\
$f_{j}(z)e_{j}-f_{j-1}(z)e_{j-1}$\end{tabular}\right),$

where $\pi_{j}=P(Y=j)$, $f_{j}(z)=(h^{-1})^{\prime}(\theta_{j}-z^{T}\beta)$
and $e_{j}$ is the $J-1$ dimensional vector with the $j$-th elements as $1$
and other elements as $0$. Here, we adopt the conventions
$f_{0}(z)=f_{J}(z)=0$ and $e_{0}=e_{J}=0$. For the proportional odds model, we
have $f_{j}(z)=exp(\theta_{j}-z^{T}\beta)/(1+exp(\theta_{j}-z^{T}\beta))^{2}$.
This representation of the information matrix is much simpler than that in
Yang, Tong and Mandal (2015) and is more suitable for the application of the
algorithms proposed in this paper. It also provides a mathematical explanation
for why the minimum number of designs points could be $d+1$ instead of
$d+J-1$. Note that the matrix $M_{z}$ based on a single point could be of rank
larger than $1$, while the rank can only be $1$ for traditional models such as
the example in Section 5.1. As a result, it is possible to collect $d+1$
points so that the summation of the information matrices based on them is of
full rank.

Now we apply the algorithm to the wine quality data set, which is publicly
available in the UCI Machine Learning Repository
(https://archive.ics.uci.edu/ml/datasets/Wine+Quality). The data were
collected to build a model regarding how the quality of the wine is affected
by various possible factors. Hence we shall be more interested in the
parameter $\beta$ in (31) rather than the $\theta_{j}$’s. The data includes
$11$ independent variables, which are all chemical components with their
levels easily measured by devices. They are all continuous variables, and we
shall standardize by their means and standard deviations so that the design
criteria could be validly adopted in the sampling process. The quality of each
wine will be rated by a wine expert by the integer scale from $0$ to $10$.
There are $N=4898$ wines in the data and Table 3 below provides the frequency
count of the quality ratings of wines. There are $7$ categories of response
observed in the full data and the categories of $3$ and $9$ are very rare. As
a result, there is a chance that some categories may not be observed in the
sample. This is not an issue since we will only be interested in $\beta$.

Quality | 3 | 4 | 5 | 6 | 7 | 8 | 9
---|---|---|---|---|---|---|---
Frequency | 20 | 163 | 1457 | 2198 | 880 | 175 | 5
Table 3: The frequency count of the quality ratings of the wines.

We use this data set only to illustrate our method instead of making any
scientific discoveries, hence we will pretend the response to be unobserved
before sampling. To implement the algorithm, we need to calculate the
information matrix of each single candidate point, which relies on the values
of the parameters. It is natural to propose a two stage procedure of sampling.
The first stage consists of a simple random sample based on which the
parameters will be estimated. The second stage tries to select the rest number
of points so that the combined sample is the most informative. Suppose MSE is
the primary target, we shall choose A-criterion as the measurement of
information. Another issue is the size allocation of the samples from the two
stages. If the percentage of the first stage sample, denoted by $r$, is too
large, the estimate of the parameters will be reliable, however, there is less
room to improve the quality of the sample on the second stage. On the other
hand, if $r$ is too small, the estimate of the parameters is subject to large
variability which induces further doubts on the determination of the second
stage sample. There is a trade off in the choice of $r$ for balancing these
two issues.

Since the data generation mechanism is unknown, we shall evaluate the
performance of the sampling methods through bootstrapping, where a big number,
say $B$, of datasets are generated by resampling from the original data with
replacement. A model fitted to the original data is treated as the true model
and the algorithm will be applied to to each resampled data to derive the
corresponding sample. The performance of a sampling method is evaluated by the
total MSE of all components of $\beta$ as defined by the quantity
$B^{-1}\sum^{B}_{b=1}||\hat{\beta}^{b}-\hat{\beta}||^{2}$, where $\hat{\beta}$
is the estimate of $\beta$ from the original data and $\hat{\beta}^{b}$ is the
estimate of $\beta$ based on the sample derived from the $b$-th bootstrapped
data. By a regular data analysis procedure, we find that model (31) is well
fitted to the original data with $11$ main effects and $17$ two way
interactions as the regressors. That means $\beta$ has $28$ components.
Further, we set the combined sample size to be $n=500$ and the number of
iterations to be $B=1000$.

Figure 3 gives the comparison between the derived sample and the random sample
in terms of the ratio of the total MSE as defined earlier when $r$ takes the
values of $0.2,0.3,...,0.9$. The smaller the ratio, the more advantage our
designed sample is as compared to the random sample. The value $r=0.1$ is not
included since the ratio is slightly larger than $1$ and hence excluded from
consideration. The smallest ratio as $0.4570$ is obtained when $r=0.4$. We
also calculate the ratio of MSE for each component of $\beta$ when $r=0.4$ and
draw a histogram of these values in Figure 4. It can been seen that only $3$
out of $28$ components have the ratio larger than one. The rest is smaller
than one with the majority below $0.6$.

Figure 3: The ratio of total MSE of all parameters between the two stage
designed sample and the random sample. Here the sample size is $500$ and the
data has $4898$ observations. The lowest ratio is $0.4570$ when the percentage
of the first stage sample is $40\%$. Figure 4: The histogram of ratios of MSE
for all parameters between the two stage designed sample and the random
sample. The sample size is $500$ and the first stage sample in the two stage
designed sample occupies $40\%$ of the total sample, namely of size $200$.

Similar to the analysis of wine quality data, we apply another data set from
the US National Health and Nutrition Examination Survey (NHANES) between 2007
and 2014 (Patel et al. 2016), a nationally representative, cross-sectional
biannual health survey, which included health questionnaire, laboratory (i.e.
urinary mercury) and clinical data. The gout as our target event, was
collected by self-reported questionnaire: "Has a doctor or other health
professional ever told you that you had gout?" Answering yes to this question
was identified positive for Gout. We remove the factors that differentiated
among individuals and those that have a majority (90%) of the observations
below a detection limit threshold as defined in the NHANES codebook. Also, we
discard the factors with missing values and the most important factors are
determined using random forest and known experience from clinical
professionals. As a result, there are a total of 8 factors with 16923
participants.

The factors are continuous, and standardized by their means and standard
deviations, which are uric acid (umol/L), red cell distribution width (%),
mean cell volume (fL), alanine aminotransferase ALT (IU/L), age in years of
the participant at the time of screening, waist circumference (cm), blood
pressure (mm Hg), and blood urea nitrogen (mg/dL). We follow the two-stage
process, where the first stage is a simple random sampling the original subset
data by $r$ between 0.3 and 0.9 to estimate the parameters in a logistic
model. The second stage is to select the remaining samples using hybrid
algorithm for D-criterion. The samples selected from two stages form to the
different sizes, in which we choose the sample sizes with 1500, 2250, 3000,
3750 and 4500 in total. After 1000 iteration bootstrapping, we obtain the
ratio of total MSE of all parameters between the two stage designed sample and
the random sample sizes from 1500 to 4500, respectively (Figure 5).

Figure 5: The ratio of total MSE of all parameters between the two stage
designed sample and the random sample. Here the majority of ratio is less than
1, except that $r$ takes values of 0.3 and 0.4 with total sample size 4500.

Figure 6 shows the histogram of the ratio of MSE between the two stage
designed sample and the random sample for each component in our logistic
model.

Figure 6: The ratio of total MSE of each parameters between the two stage
designed sample and the random sample. Here there is 8.57% component whose the
ratio of total MSE is larger than 1.

To distinguish the selection bias, we calculate the ratio of MSE for each
component including only variance. Figure 7 represents the histogram of the
ratio of MSE between the two stage designed sample and the random sample for
each component incorporating variance. There is only the intercept whose the
ratio of MSE for variance is larger than 1.

Figure 7: The ratio of MSE of each parameters between the two stage designed
sample and the random sample. Here there is 3.9% component whose the ratio of
MSE for variance is larger than 1.

## References

  * Brinker [2003] Klaus Brinker.  Incorporating diversity in active learning with support vector machines.  In _Proceedings of the 20th International Conference on Machine Learning (ICML-03)_ , pages 59–66, 2003.
  * Campbell et al. [2000] Colin Campbell, Nello Cristianini, Alex Smola, et al.  Query learning with large margin classifiers.  In _ICML_ , volume 20, 2000.
  * Chen and Krause [2013] Yuxin Chen and Andreas Krause.  Near-optimal batch mode active learning and adaptive submodular optimization.  In _International Conference on Machine Learning_ , pages 160–168. PMLR, 2013.
  * Deng et al. [2009] Xinwei Deng, V Roshan Joseph, Agus Sudjianto, and CF Jeff Wu.  Active learning through sequential design, with applications to detection of money laundering.  _Journal of the American Statistical Association_ , 104(487):969–981, 2009.
  * Deng et al. [2022] Yujia Deng, Yubai Yuan, Haoda Fu, and Annie Qu.  Query-augmented active metric learning.  _Journal of the American Statistical Association_ , pages 1–14, 2022\.
  * Fedorov [1989] Valerii V Fedorov.  Optimal design with bounded density: Optimization algorithms of the exchange type.  _Journal of Statistical Planning and Inference_ , 22(1):1–13, 1989.
  * Freund et al. [1997] Yoav Freund, H Sebastian Seung, Eli Shamir, and Naftali Tishby.  Selective sampling using the query by committee algorithm.  _Machine Learning_ , 28(2):133–168, 1997.
  * Ghorbani et al. [2021] Amirata Ghorbani, James Zou, and Andre Esteva.  Data shapley valuation for efficient batch active learning.  _arXiv preprint arXiv:2104.08312_ , 2021.
  * Guo and Schuurmans [2007] Yuhong Guo and Dale Schuurmans.  Discriminative batch mode active learning.  _Advances in Neural Information Processing Systems_ , 20, 2007.
  * Hoi et al. [2006a] Steven CH Hoi, Rong Jin, and Michael R Lyu.  Large-scale text categorization by batch mode active learning.  In _Proceedings of the 15th International Conference on World Wide Web_ , pages 633–642, 2006a.
  * Hoi et al. [2006b] Steven CH Hoi, Rong Jin, Jianke Zhu, and Michael R Lyu.  Batch mode active learning and its application to medical image classification.  In _Proceedings of the 23rd International Conference on Machine Learning_ , pages 417–424, 2006b.
  * Houlsby et al. [2011] Neil Houlsby, Ferenc Huszár, Zoubin Ghahramani, and Máté Lengyel.  Bayesian active learning for classification and preference learning.  _arXiv preprint arXiv:1112.5745_ , 2011.
  * Ienco et al. [2013] Dino Ienco, Albert Bifet, Indrė Žliobaitė, and Bernhard Pfahringer.  Clustering based active learning for evolving data streams.  In _International Conference on Discovery Science_ , pages 79–93. Springer, 2013.
  * Kiefer [1974] Jack Kiefer.  General equivalence theory for optimum designs (approximate theory).  _The Annals of Statistics_ , pages 849–879, 1974.
  * Kirsch et al. [2019] Andreas Kirsch, Joost Van Amersfoort, and Yarin Gal.  Batchbald: Efficient and diverse batch acquisition for deep bayesian active learning.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Lewis and Gale [1994] David D Lewis and William A Gale.  A sequential algorithm for training text classifiers.  In _SIGIR’94_ , pages 3–12. Springer, 1994.
  * Li et al. [2012] Xiao Li, Da Kuang, and Charles X Ling.  Active learning for hierarchical text classification.  In _Pacific-Asia Conference on Knowledge Discovery and Data Mining_ , pages 14–25. Springer, 2012.
  * Müller and Pázman [1998] Werner G Müller and Andrej Pázman.  Design measures and approximate information matrices for experiments without replications.  _Journal of Statistical Planning and Inference_ , 71(1-2):349–362, 1998.
  * Ouyang et al. [2016a] Liwen Ouyang, Daniel W Apley, and Sanjay Mehrotra.  A design of experiments approach to validation sampling for logistic regression modeling with error-prone medical records.  _Journal of the American Medical Informatics Association_ , 23(e1):e71–e78, 2016a.
  * Ouyang et al. [2016b] Liwen Ouyang, Daniel W Apley, and Sanjay Mehrotra.  Designed sampling from large databases for controlled trials.  _IIE Transactions_ , 48(12):1087–1097, 2016b.
  * Ouyang et al. [2017] Liwen Ouyang, Daniel W Apley, and Sanjay Mehrotra.  Batch sample design from databases for logistic regression.  _Quality and Reliability Engineering International_ , 33(1):87–101, 2017.
  * Patel et al. [2016] Chirag J Patel, Nam Pho, Michael McDuffie, Jeremy Easton-Marks, Cartik Kothari, Isaac S Kohane, and Paul Avillach.  A database of human exposomes and phenomes from the us national health and nutrition examination survey.  _Scientific Data_ , 3(1):1–10, 2016.
  * Perevozskaya et al. [2003] Inna Perevozskaya, William F Rosenberger, and Linda M Haines.  Optimal design for the proportional odds model.  _Canadian Journal of Statistics_ , 31(2):225–235, 2003.
  * Pronzato [2004] Luc Pronzato.  A minimax equivalence theorem for optimum bounded design measures.  _Statistics & Probability Letters_, 68(4):325–331, 2004.
  * Pronzato [2006] Luc Pronzato.  On the sequential construction of optimum bounded designs.  _Journal of Statistical Planning and Inference_ , 136(8):2783–2804, 2006.
  * Ravi et al. [2016] Sathya Narayanan Ravi, Vamsi Ithapu, Sterling Johnson, and Vikas Singh.  Experimental design on a budget for sparse linear models and applications.  In _International Conference on Machine Learning_ , pages 583–592. PMLR, 2016.
  * Ren et al. [2021] Pengzhen Ren, Yun Xiao, Xiaojun Chang, Po-Yao Huang, Zhihui Li, Brij B Gupta, Xiaojiang Chen, and Xin Wang.  A survey of deep active learning.  _ACM Computing Surveys (CSUR)_ , 54(9):1–40, 2021\.
  * Sahm and Schwabe [2001] Michael Sahm and Rainer Schwabe.  A note on optimal bounded designs.  In _Optimum Design 2000_ , pages 131–140. Springer, 2001.
  * Scheffer et al. [2001] Tobias Scheffer, Christian Decomain, and Stefan Wrobel.  Active hidden markov models for information extraction.  In _International Symposium on Intelligent Data Analysis_ , pages 309–318. Springer, 2001.
  * Schohn and Cohn [2000] Greg Schohn and David Cohn.  Less is more: Active learning with support vector machines.  In _ICML_ , volume 2, page 6. Citeseer, 2000.
  * Sebastiani and Wynn [2000] Paola Sebastiani and Henry P Wynn.  Maximum entropy sampling and optimal bayesian experimental design.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 62(1):145–157, 2000.
  * Sener and Savarese [2018] Ozan Sener and Silvio Savarese.  Active learning for convolutional neural networks: A core-set approach.  In _International Conference on Learning Representations_ , 2018.
  * Settles [2009] Burr Settles.  In _Active learning literature survey_. University of Wisconsin-Madison Department of Computer Sciences, 2009.
  * Tong and Koller [2001] Simon Tong and Daphne Koller.  Support vector machine active learning with applications to text classification.  _Journal of Machine Learning Research_ , 2:45–66, 2001\.
  * Uciński [2015] Dariusz Uciński.  An algorithm for construction of constrained d-optimum designs.  In _Stochastic Models, Statistics and Their Applications_ , pages 461–468. Springer, 2015.
  * Vandenberghe et al. [1998] Lieven Vandenberghe, Stephen Boyd, and Shao-Po Wu.  Determinant maximization with linear matrix inequality constraints.  _SIAM Journal on Matrix Analysis and Applications_ , 19(2):499–533, 1998.
  * Wang and Hua [2011] Meng Wang and Xian-Sheng Hua.  Active learning in multimedia annotation and retrieval: A survey.  _ACM Transactions on Intelligent Systems and Technology (TIST)_ , 2(2):1–21, 2011.
  * Wei et al. [2015] Kai Wei, Rishabh Iyer, and Jeff Bilmes.  Submodularity in data subset selection and active learning.  In _International Conference on Machine Learning_ , pages 1954–1963. PMLR, 2015.
  * Wynn [1977] H.P. Wynn.  Optimum designs for finite populations sampling.  In _Statistical Decision Theory and Related Topics II_ , pages 471–478. Academic Press, New York, 1977.
  * Wynn [1982] H.P. Wynn.  Optimum submeasures with application to finite population sampling.  In _Statistical Decision Theory and Related Topics III_ , pages 485–495. Elsevier, 1982.
  * Xu et al. [2003] Zhao Xu, Kai Yu, Volker Tresp, Xiaowei Xu, and Jizhi Wang.  Representative sampling for text classification using support vector machines.  In _European Conference on Information Retrieval_ , pages 393–407. Springer, 2003.
  * Yang et al. [2017] Jie Yang, Liping Tong, and Abhyuday Mandal.  D-optimal designs with ordered categorical data.  _Statistica Sinica_ , pages 1879–1902, 2017.
  * Yang et al. [2013] Min Yang, Stefanie Biedermann, and Elina Tang.  On optimal designs for nonlinear models: a general and efficient algorithm.  _Journal of the American Statistical Association_ , 108(504):1411–1420, 2013.
  * Zhang and Oles [2000] Tong Zhang and F Oles.  The value of unlabeled data for classification problems.  In _Proceedings of the 17th International Conference on Machine Learning,(Langley, P., ed.)_ , volume 20. Citeseer, 2000.
  * Zhdanov [2019] Fedor Zhdanov.  Diverse mini-batch active learning.  _arXiv preprint arXiv:1901.05954_ , 2019.