On greedy heuristics for computing D-efficient saturated subsets

TL;DR

This paper introduces two new greedy algorithms for selecting saturated subsets of vectors that maximize the volume of the parallelotope, improving $D$-optimal experimental design and ellipsoid algorithms.

Contribution

The paper proposes two novel greedy heuristics for constructing $D$-optimal saturated subsets, with geometric insights and comparisons to existing methods.

Findings

The heuristics outperform random and regularized greedy methods.

The methods are effective for large sets and alternative criteria.

Analytic and geometric analysis supports the heuristics' effectiveness.

Abstract

Let $\mathcal{F}$ be a set consisting of $n$ real vectors of dimension $m \leq n$. For any saturated, i.e., $m$-element, subset $\mathcal{S}$ of $\mathcal{F}$, let $\mathrm{vol}(\mathcal{S})$ be the volume of the parallelotope formed by the vectors of $\mathcal{S}$. A set $\mathcal{S}^*$ is called a $D$-optimal saturated subset of $\mathcal{F}$, if it maximizes $\mathrm{vol}(\mathcal{S})$ among all saturated subsets of $\mathcal{F}$. In this paper, we propose two greedy heuristics for the construction of saturated subsets performing well with respect to the criterion of $D$-optimality: an improvement of the method suggested by Galil and Kiefer for the initiation of $D$-optimal experimental design algorithms, and a modification of the Kumar-Yildirim method, the original version of which was proposed for the initiation of the minimum-volume enclosing ellipsoid algorithms. We provide geometric and analytic insights into the two methods, and compare them to the commonly used random and regularized greedy heuristics. We also suggest variants of the greedy methods for a large set $\mathcal{F}$, for the construction of $D$-efficient non-saturated subsets, and for alternative optimality criteria.

Equations

$$\varphi :2^{\mathcal{F}}\to [0,\infty );\varphi (\mathcal{S})={\det }^{1/m}(\mathbf{M}(\mathcal{S})).$$

$${\xi }^{\ast }\in {\mathrm{argmax}}_{\xi \in \Xi }{\det }^{1/m}\left({\int }_{\mathcal{F}}\mathbf{f}{\mathbf{f}}^{\prime }d\xi (\mathbf{f})\right),$$

$${\varphi }_s^{\ast }:=s{\det }^{1/m}\left({\int }_{\mathcal{F}}\mathbf{f}{\mathbf{f}}^{\prime }d{\xi }^{\ast }(\mathbf{f})\right).$$

$${\mathbf{f}}_{(j)}\in {\mathrm{argmax}}_{\mathbf{f}\in \mathcal{F}\setminus {\mathcal{S}}_j}\varphi (\mathbf{M}({\mathcal{S}}_j)+\mathbf{f}{\mathbf{f}}^{\prime }+\delta {\mathbf{I}}_m).$$

$${\mathbf{f}}_{(j)}\in {\mathrm{argmax}}_{\mathbf{f}\in \mathcal{F}\setminus {\mathcal{S}}_j}{\mathbf{f}}^{\prime }(\mathbf{M}({\mathcal{S}}_j)+\delta {\mathbf{I}}_m{)}^{-1}\mathbf{f}.$$

$${\mathbf{f}}_{(j)}\in {\mathrm{argmax}}_{\mathbf{f}\in \mathcal{F}\setminus {\mathcal{S}}_j}\prod _{i=1}^{j+1}{\lambda }_i(\mathbf{M}({\mathcal{S}}_j)+\mathbf{f}{\mathbf{f}}^{\prime }),$$

$$\mathbf{f}_{(j)}\in\mathrm{argmax}_{\mathbf{f}\in\mathcal{F}\setminus\mathcal{S}_{j}}\det\big{(}\mathbf{F}(\mathcal{S}_{j}\cup\mathbf{f})\mathbf{F}^{\prime}(\mathcal{S}_{j}\cup\mathbf{f})\big{)},$$

$${\mathbf{f}}_{(j)}\in {\mathrm{argmax}}_{\mathbf{f}\in \mathcal{F}\setminus {\mathcal{S}}_j}{\mathbf{f}}^{\prime }({\mathbf{I}}_m-\mathbf{P}({\mathcal{S}}_j))\mathbf{f},$$

$$\prod _{i=1}^{j+1}{\lambda }_i(\mathbf{M}({\mathcal{S}}_j)+\mathbf{f}{\mathbf{f}}^{\prime })=\det \left[\begin{array}{cc}{\mathbf{V}}^{\prime }\mathbf{V}& {\mathbf{V}}^{\prime }\mathbf{f}\\ {\mathbf{f}}^{\prime }\mathbf{V}& {\mathbf{f}}^{\prime }\mathbf{f}\end{array}\right],$$

$${\mathbf{f}}^{\prime }({\mathbf{I}}_m-\mathbf{P}({\mathcal{S}}_j))\mathbf{f}=\Vert ({\mathbf{I}}_m-\mathbf{P}({\mathcal{S}}_j))\mathbf{f}{\Vert }^2,$$

$${\mathrm{argmax}}_{\mathbf{f}\in \mathcal{F}\setminus {\mathcal{S}}_j}{\mathbf{f}}^{\prime }({\mathbf{I}}_m-{\mathbf{F}}^{\prime }(\mathbf{F}{\mathbf{F}}^{\prime }+\delta {\mathbf{I}}_j{)}^{-1}\mathbf{F})\mathbf{f}.$$

$${\delta }^{m-j-1}\det (\left[\begin{array}{cc}\mathbf{F}{\mathbf{F}}^{\prime }& \mathbf{Ff}\\ {\mathbf{f}}^{\prime }{\mathbf{F}}^{\prime }& {\mathbf{f}}^{\prime }\mathbf{f}\end{array}\right]+\delta {\mathbf{I}}_{j+1}).$$

$${\mathrm{argmax}}_{\mathbf{f}\in \mathcal{G}}\varphi (\mathbf{M}+\delta {\mathbf{I}}_m+\mathbf{f}{\mathbf{f}}^{\prime })\subseteq {\mathrm{argmax}}_{\mathbf{f}\in \mathcal{G}}{\mathbf{f}}^{\prime }({\mathbf{I}}_m-\mathbf{P})\mathbf{f}$$

$${\mathrm{argmax}}_{\mathbf{f}\in \mathcal{G}}{\mathbf{f}}^{\prime }({\mathbf{I}}_m-\mathbf{P})\mathbf{f}=\{\mathbf{f}\in \mathcal{G}:{\mathbf{f}}^{\prime }({\mathbf{I}}_m-\mathbf{P})\mathbf{f}>h-ϵ\}.$$

$$\underset{\delta \to 0^{+}}{\lim }{\mathbf{f}}^{\prime }({\mathbf{I}}_m-{\mathbf{F}}^{\prime }(\mathbf{F}{\mathbf{F}}^{\prime }+\delta {\mathbf{I}}_r{)}^{-1}\mathbf{F})\mathbf{f}={\mathbf{f}}^{\prime }({\mathbf{I}}_m-\mathbf{P})\mathbf{f}$$

$$|{\mathbf{f}}^{\prime }({\mathbf{I}}_m-{\mathbf{F}}^{\prime }(\mathbf{F}{\mathbf{F}}^{\prime }+\delta {\mathbf{I}}_r{)}^{-1}\mathbf{F})\mathbf{f}-{\mathbf{f}}^{\prime }({\mathbf{I}}_m-\mathbf{P})\mathbf{f}|<ϵ/2$$

$${\tilde{\mathbf{f}}}_j$$

$$=(\mathbf{I}-{\tilde{\mathbf{P}}}_{j-1}-\dots -{\tilde{\mathbf{P}}}_1){\mathbf{f}}_j,$$

$$=(\mathbf{I}-{\tilde{\mathbf{P}}}_{j-1})\cdots (\mathbf{I}-{\tilde{\mathbf{P}}}_1){\mathbf{f}}_j.$$

$$\frac{\varphi ({\mathcal{S}}_{GK})}{\varphi ({\mathcal{S}}^{\ast })}\ge \frac{\varphi (\tilde{\mathcal{S}})}{\varphi ({\mathcal{S}}_B)}=\frac{{\prod }_i\Vert {\tilde{\mathbf{f}}}_i{\Vert }^{2/m}}{m{\prod }_i\Vert {\tilde{\mathbf{f}}}_i{\Vert }^{2/m}}=\frac{1}{m}.$$

$${\mathbf{f}}_{(j)}\in {\mathrm{argmax}}_{\mathbf{f}\in \mathcal{F}\setminus {\mathcal{S}}_j}|{\mathbf{f}}^{\prime }\mathbf{b}|$$

$$\mathrm{eff}({\mathcal{S}}_{KY})\ge \frac{\pi }{4m(\Gamma (1+m/2){)}^{2/m}}\ge \frac{\pi }{2m^2}.$$

$$\mathrm{vol}(K)\le |\det (2\mathbf{F}({\mathcal{S}}_{KY}))|=2^m\det (\mathbf{M}({\mathcal{S}}_{KY}){)}^{1/2},$$

$$\frac{1}{m^{m/2}}\mathrm{vol}({\mathrm{MVEE}}^{\ast })$$

$$\mathrm{vol}({\mathrm{MVEE}}^{\ast })=\eta (\det \mathbf{M}({\mathcal{S}}^{\ast }){)}^{1/2},$$

$$\frac{\eta (\det \mathbf{M}({\mathcal{S}}^{\ast }){)}^{1/2}}{m^{m/2}}\le 2^m(\det (\mathbf{M}({\mathcal{S}}_{KY}){)}^{1/2}$$

$$\frac{{\Phi }_D(\mathbf{M}({\mathcal{S}}_{KY}))}{{\Phi }_D(\mathbf{M}({\mathcal{S}}^{\ast }))}\ge \frac{{\eta }^{2/m}}{4m}=\frac{\pi }{4m(\Gamma (1+m/2){)}^{2/m}}.$$

$$\mathrm{eff}({\mathcal{S}}_{KY})\ge \frac{\pi }{4m(\Gamma (1+m/2){)}^{2/m}}\approx \frac{e\pi }{2m^2}.$$

$$P[\mathcal{R}\text{ is non-singular}]=P[{\mathbf{g}}_1,\dots ,{\mathbf{g}}_{km}\text{ span }{\mathbb{R}}^m]\ge$$

$$P[{\mathbf{h}}_1,\dots ,{\mathbf{h}}_{km}\text{ span }{\mathbb{R}}^m]\ge 1-P[{\cap }_t{A}_t]\ge 1-p^k,$$

$${\mathrm{argmax}}_{\mathbf{f}\in \mathcal{F}\setminus \mathcal{S}}{\left(\sum _{i=1}^{r+1}{\lambda }_i^{-1}(\mathbf{M}(\mathcal{S})+\mathbf{f}{\mathbf{f}}^{\prime })\right)}^{-1}=$$

$${\mathrm{argmax}}_{\mathbf{f}\in \mathcal{F}\setminus \mathcal{S}}\frac{{\mathbf{f}}^{\prime }({\mathbf{I}}_m-\mathbf{P}(\mathcal{S}))\mathbf{f}}{1+{\mathbf{f}}^{\prime }{\mathbf{M}}^{+}(\mathcal{S})\mathbf{f}},$$

Full text

Keywords: $D$-optimality; Optimal experimental design; Minimum-volume enclosing ellipsoid; Subsampling; Greedy heuristic

1 Introduction

Consider a set $\mathcal{F}$ consisting of vectors $\mathbf{f}_{1},\ldots,\mathbf{f}_{n}\in\mathbb{R}^{m}$ which we will call “regressors”, and let $2\leq m\leq n$. Our aim is to select an $s$-element subset of $\mathcal{F}$, $s\leq n$, focusing on the “saturated” subsets, that is, $s=m$. The saturated subsets are natural initial solutions of various algorithms applied in statistics and computational geometry; note that in some of the applications the set $\mathcal{F}$ is the result of a random process, in others it is a deterministic set obtained by an algebraic construction. Typically $n\gg m$.

To measure the quality of the subsets of $\mathcal{F}$, saturated or non-saturated, we consider an “information-based” $D$-optimality criterion, which is widely used in the theory of optimal experimental design (e.g., [5], [20], [22], and [2]).

Define the information matrix of $\mathcal{S}\subseteq\mathcal{F}$ as $\mathbf{M}(\mathcal{S})=\sum_{\mathbf{f}\in\mathcal{S}}\mathbf{f}\mathbf{f}^{\prime}$. If $\mathcal{S}=\emptyset$, set $\mathbf{M}(\mathcal{S})=\mathbf{0}_{m\times m}$. Then, the $D$-criterion is

[TABLE]

The $D$-optimal $s$-element subset $\mathcal{S}^{*}$ is any subset of $\mathcal{F}$ that maximizes $\phi(\mathcal{S})$ in the class of all $s$-element subsets of $\mathcal{F}$.

To make the problem meaningful, we assume that $\mathcal{F}$ is “non-singular” in the sense that $\mathbf{M}(\mathcal{F})$ is a non-singular matrix, i.e., the $n\times m$ matrix $\mathbf{F}(\mathcal{F})$ with rows $\mathbf{f}^{\prime}_{i}$, $\mathbf{f}_{i}\in\mathcal{F}$, has full column rank. This assumption is clearly equivalent to the existence of a saturated non-singular subset of $\mathcal{F}$.

Note that for the case of saturated subsets, we can equivalently define $D$-optimality as follows: A saturated subset $\mathcal{S}^{*}$ is $D$-optimal if and only if it maximizes $|\det(\mathbf{F}(\mathcal{S}))|$. Geometrically, a $D$-optimal saturated subset corresponds to the maximum-volume parallelotope formed by an $m$-tuple of vectors from $\mathcal{F}$. For a general $s\geq m$, a $D$-optimal $s$-element subset of $\mathcal{F}$ maximizes the mean squared volume of the $\binom{s}{m}$ parallelotopes formed by the vectors of $\mathcal{S}$. This is a direct consequence of the Cauchy-Binet formula.

In the theory of optimal design of experiments, a $D$-optimal $s$-element subset is called a replication-free $D$-optimal exact experimental design of size $s$ (e.g., [23], cf. [6])111For the saturated case, the problem of $D$-optimal exact design with replication and the problem of $D$-optimal exact designs without replications are equivalent.. Computing such sets, or designs, is generally a difficult problem of discrete optimization; for a review of optimal design algorithms, see, e.g., Chapter 12 of [2], [9] and [18].

Let $\Xi$ be the set of all probability measures on $\mathcal{F}$. Clearly, if we find some

[TABLE]

then we obtain the following upper bound on $\phi(\mathcal{S}^{*})$:

[TABLE]

The probability measure $\xi^{*}$ from (2) is called the $D$-optimal approximate design for the linear regression model with regressors $\mathcal{F}$. The $D$-optimal approximate designs and, consequently, the bounds $\phi^{*}$, can be computed via a number of specialized algorithms of convex optimization (e.g., [34], [33], [11], cf. Chapter 9 of [21]).

The $D$-efficiency of a set $\mathcal{S}_{1}$ relative to a non-singular set $\mathcal{S}_{2}$ is defined as $\phi(\mathcal{S}_{1})/\phi(\mathcal{S}_{2})$. If $\mathcal{S}^{*}$ is a $D$-optimal $s$-element subset, then the relative efficiency of an $s$-element subset $\mathcal{S}$ with respect to $\mathcal{S}^{*}$ will be denoted by $\mathrm{eff}(\mathcal{S})$. Note that $\phi(\mathcal{S})/\phi^{*}_{s}\leq\mathrm{eff}(\mathcal{S})\leq 1$.

The algorithms for $D$-optimal subsets can also be used for various purposes other than constructing maximum-volume (configurations of) parallelepippeds or designs of statistical experiments. For instance, we may be interested in computing the ellipsoid $\mathrm{MVEE}_{0}(\mathcal{X})$222MVEE stands for “Minimum-Volume Enclosing Ellipsoid”. which has the minimum volume among all ellipsoids covering the set $\mathcal{X}=\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\}\subset\mathbb{R}^{d}$ and centred in $\mathbf{0}_{d}$, or the ellipsoid $\mathrm{MVEE}(\mathcal{X})$ with the minimum volume among completely all ellipsoids covering $\mathcal{X}$. It turns out that the computation the MVEEs is equivalent to the computation of a $D$-optimal approximate design for the statistical model with regressors $\mathbf{f}_{i}=\mathbf{x}_{i}$ in the case of $\mathrm{MVEE}_{0}(\mathcal{X})$, and for the statistical model with regressors $\mathbf{f}_{i}=(\mathbf{x}_{i}^{\prime},1)^{\prime}$ in the case of $\mathrm{MVEE}(\mathcal{X})$. For details, see [26], [28], [29], [1], [30] or the introduction of [11].

Another application of $D$-efficient subsets is in the information-based data subselection. Consider a large sample consisting of the points $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\in\mathbb{R}^{d}$ (e.g., the vectors of covariates) and let $\mathbf{f}_{i}=f(\mathbf{x}_{i})$, $i=1,\ldots,n$, be the corresponding regressors, where $f$ is an appropriate “regression function”. Then, one can select an informative subsample of size $s$ by choosing $s$ regressors $\mathbf{f}_{i}$ (and the corresponding points $\mathbf{x}_{i}$) that make up a $D$-optimal or $D$-efficient subset (e.g., [4] and [32]).

To compute an information-based subsample, [32] provides a fast specialized deterministic algorithm, primarily for sizes $s=2kd$, $k\in\mathbb{N}$. Note, however, that the algorithm, although fast and potentially very useful for big data, produces subsets that are not saturated, and generally not particularly $D$-efficient, in some cases even singular. The authors of [4] also construct a subsample motivated by criteria used in experimental design, by solving two optimization problems: the first one on a possibly continuous superset of $\mathcal{F}$, and then second one on $\mathcal{F}$.

As shown above, there is a large number of algorithms for constructing $D$-optimal or $D$-efficient subsets of a general size. However, these algorithms often require a small and non-singular initiation subset, preferably as efficient as possible. In particular, the saturated subsets can be used to initiate algorithms for both optimal approximate designs and optimal exact designs. In many such algorithms, such as [34], [33], and [2], the the initial solution suggested by the authors is either random or based on a regularized greedy heuristic (cf. Subsections 2.1 and 2.2), which may not be the ideal method. It was the main motivation of this paper to analyse and compare alternative methods to generate such initial subsets.

2 Greedy heuristics for saturated subsets

A general form of the greedy procedure for the construction of saturated subsets is the following:

**Scheme of greedy heuristics for saturated subsets

**1. Set $j\leftarrow 0$ and $\mathcal{S}_{0}\leftarrow\emptyset$

2. Select $\mathbf{f}_{(j)}$ from $\mathcal{F}\setminus\mathcal{S}_{j}$   **(S)

**3. Set $\mathcal{S}_{j+1}\leftarrow\mathcal{S}_{j}\cup\{\mathbf{f}_{(j)}\}$

5. If $j+1<m$ set $j\leftarrow j+1$ and return to 2

6. Output $\mathcal{S}_{m}$

The main difficulty for the greedy construction of saturated subsets is that during the process of “saturation” we work with subsets with fewer than $m$ elements, and such subsets $\mathcal{S}$ are necessarily singular. That is, $\phi(\mathcal{S})=0$ and the quality of such subsets cannot be directly compared via $\phi$. In particular, the general greedy heuristic, which in each step adds the regressor $\mathbf{f}_{(j)}$ that provides the highest increase in the objective function (e.g., [19]) cannot be used. In [24], the author circumvented this deficiency by considering optimality criteria that are not constantly zero for singular subsets, but we will not use this approach.

The methods described next solve the saturated subset problem in different ways, leading to various specifications of the critical step 2 which we call a saturation step and denote it by (S).

The heuristic outlined above is sometimes called a “forward” method. However, there also exist “drop” or “backward” methods (e.g., [23] and [2], Section 12.4), where we begin with the full set $\mathcal{F}$ and delete the least informative regressors one by one. These, however, turn out to be slow, because their asymptotic complexity is $O(n^{2}-nm)$ and in applications $n$ is usually a large number.

2.1 Random subsampling

The simplest method for obtaining a saturated subset is to select $m$ regressors uniformly at random. Formally, in step (S) we select a random $\mathbf{f}_{(j)}$ from $\mathcal{F}\setminus\mathcal{S}_{j}$. This method is extremely fast, but the random saturated subsets tend to have low efficiency. It is even possible that such a saturated $\mathcal{S}$ is singular. The case of a singular random subset is not only possible in theory, but also in applications, as shown in Example 1.

Example 1**.**

Consider the classical factorial linear regression model with $m$ factors, each with levels $-1,1$, with main factors and no constant term. In this situation the set of regressors is $\mathcal{F}=\{-1,1\}^{m}$. Suppose that we intend to construct a saturated experimental design at random. Then, the probability of obtaining a singular subset (i.e., a singular design) can be undesirably high. For example, for $m=3,4,5$, and $6$ this probability is $\frac{3}{7}\approx 0.429,\frac{223}{455}\approx 0.490,\frac{3285}{6293}\approx 0.522$, and $\frac{175795}{334707}\approx 0.525$.

To obtain more efficient subsets, methods for nonuniform random subsampling were proposed. In particular, the leverage-based random subsampling (e.g., [17] and [16]) selects $\mathbf{f}_{(j)}=\mathbf{f}_{i}$ from $\mathcal{F}\setminus\mathcal{S}_{j}$ with probability $\pi_{i}$ that is proportional to its so-called leverage score $\mathbf{f}_{i}^{\prime}\mathbf{M}^{-1}(\mathcal{F})\mathbf{f}_{i}$.333Leverages based on subsets of $\mathcal{F}$, called sensitivities in optimal design of experiments (e.g., [5]) naturally occur in most $D$-optimum design algorithms because they have an important analytic and statistical interpretation. Note that a non-uniform subsampling can be much more computationally demanding than the uniform sampling. Moreover, the leveraging methods are also not guaranteed to produce a non-singular subset, especially if one seeks a saturated subset: Example 2 demonstrates this.

Example 2**.**

In Example 1, the information matrix $\mathbf{M}(\mathcal{F})$ is proportional to $\mathbf{I}_{m}$ and $\mathbf{f}_{i}^{\prime}\mathbf{f}_{i}=m$ for all $i$. Therefore, all regressors have the same leverage score which means that the leveraging method is equivalent to the uniform random subset selection, including the probabilities of generating a singular subset as listed in Example 1.

A non-zero probability of generating a singular subset for a given $\mathcal{F}$ is not necessarily a damaging feature444Unless, of course, the probability of generating a singular subset is $1$ or very close to $1$.; we can simply repeatedly use the method until we hit a non-singular subset. See Subsection 3.3 for further discussion on this “multi-start” approach.

2.2 The regularized greedy heuristic (RGH)

Except for the simple or the weighted random choice of a saturated set $\mathcal{S}\subseteq\mathcal{F}$, the most common is the procedure based on a regularization of the information matrix of a sub-saturated set (see, e.g., Section 11.2 of [2]; cf. [31]). The regularized greedy heuristic (RGH) specifies step (S) as

[TABLE]

Here, $\delta$ is a small positive constant; in the context of experimental design, [2] suggest to take $\delta$ between $10^{-4}$ and $10^{-6}$. The matrix determinant lemma directly implies that (4) is equivalent to

[TABLE]

A short R-implementation of the RGH based on an efficient computation of (5) is given in Appendix. The asymptotic complexity of this algorithm is $O(nm^{3})$.

The numerical experience shows that the RGH usually finds a reasonably $D$-efficient saturated subset $\mathcal{S}$. However, a little known fact is that, similarly to the random methods, the RGH is not guaranteed to provide a non-singular saturated subset, even if one exists:

Example 3**.**

Consider the RGH with $\delta=10^{-4}$ for $\mathbf{f}_{a}=(1,0,0)^{\prime}$, $\mathbf{f}_{b}=(0,1,0)^{\prime}$, $\mathbf{f}_{c}=(1,1,0)^{\prime}$ and $\mathbf{f}_{d}=(0,0,10^{-5})^{\prime}$. Then, RGH first chooses $\mathbf{f}_{c}$, then $\mathbf{f}_{a}$ (or $\mathbf{f}_{b}$) and then $\mathbf{f}_{b}$ (or $\mathbf{f}_{a}$). The resulting subset $\mathcal{S}=\{\mathbf{f}_{a},\mathbf{f}_{b},\mathbf{f}_{c}\}$ is singular, although a non-singular subset $\{\mathbf{f}_{a},\mathbf{f}_{b},\mathbf{f}_{d}\}$ of size $3$ does exist.

That is, we cannot fully rely on the non-singularity of the design produced by RGH. Moreover, unlike the random methods, the RGH is deterministic (cf. Subsection 3.3), therefore it has no meaning to run it multiple times in search for an improvement. On the other hand, RGH can be also applied to criteria other than $D$-optimality in a straightforward manner. We will show in the following sections that to construct a $D$-efficient saturated subset, it may be preferable to use modifications of the original Galil-Kiefer and the Kumar-Yildirim methods.

2.3 The Galil-Kiefer method (GKM)

In [7] (see also [2], Section 12.4) Galil and Kiefer proposed to specify step (S) as

[TABLE]

where $\lambda_{1}(\mathbf{M})\geq\cdots\geq\lambda_{m}(\mathbf{M})$ denote the eigenvalues of a non-negative definite $m\times m$ matrix $\mathbf{M}$. That is, instead of maximizing the product of all eigenvalues $\phi(\mathcal{S})$, which is zero, we maximize the product of the $j+1$ largest eigenvalues of $\mathbf{M}(\mathcal{S})$, which can be non-zero. Let $\mathcal{S}\subseteq\mathcal{F}$ be an $r$-element set. Using the well known fact (e.g. 6.54 of [25]) that the $r$ largest eigenvalues of $\mathbf{M}(\mathcal{S})=\mathbf{F}^{\prime}(\mathcal{S})\mathbf{F}(\mathcal{S})$ are equal to the $r$ eigenvalues of $\mathbf{F}(\mathcal{S})\mathbf{F}^{\prime}(\mathcal{S})$, we see that (6) is equivalent to

[TABLE]

which was also noted in [7]. Formula (7) replaces the computation of eigenvalues by the simpler computation of determinants.

Although the approach of Galil and Kiefer appears to be reasonable, it seems that it has not been much analysed in the literature and used in applications. This is probably caused by the fact that the determinant form (7) is still computationally expensive. The paper [7] does not provide a more efficient solution for the update of $\mathcal{S}_{j}$; the authors simply state that the standard formulas can be “implemented with appropriate changes”.

In the sequel, we show a more efficient and geometrically interpretable form of the optimization step. The basic insight is that the determinant form of the GKM has a simple geometric interpretation - we choose $\mathbf{f}_{(j)}$ to maximize the $j$-dimensional volume of the parallelotope spanned in $\mathbb{R}^{m}$ by the vectors of $\mathcal{S}_{j}\cup\mathbf{f}$, $\mathbf{f}\in\mathcal{F}\setminus\mathcal{S}_{j}$. Alternatively, we can formulate it as follows.

Theorem 1**.**

The GKM always produces a non-singular saturated subset of $\mathcal{F}$ if such a subset exists. Moreover, the optimization (6) (or (7)) is equivalent to

[TABLE]

where $\mathbf{P}(\mathcal{S}_{j})$ is the projector on the column space of $\mathbf{M}(\mathcal{S}_{j})$.

Proof.

If $\mathbf{f}\in\mathcal{C}(\mathbf{F}^{\prime}(\mathcal{S}_{j}))$ then $\mathbf{F}(\mathcal{S}_{j}\cup\mathbf{f})\mathbf{F}^{\prime}(\mathcal{S}_{j}\cup\mathbf{f})$ is not of full rank and its determinant is zero; here $\mathcal{C}$ denotes the column space. It follows that the GKM always chooses $\mathbf{f}\not\in\mathcal{C}(\mathbf{F}^{\prime}(\mathcal{S}_{j}))$ and by induction, we obtain that the rank of $\mathbf{M}(\mathcal{S}_{j})$ is $j$ for $j=1,\ldots,m$. Therefore, $\mathbf{M}(\mathcal{S}_{j})$ can be expressed as $\mathbf{V}\mathbf{V}^{\prime}$, where $\mathbf{V}$ is an $m\times j$ matrix of rank $j$. The positive eigenvalues of $\mathbf{M}(\mathcal{S}_{j})+\mathbf{f}\mathbf{f}^{\prime}=(\mathbf{V},\mathbf{f})(\mathbf{V},\mathbf{f})^{\prime}$ are equal to the positive eigenvalues of $(\mathbf{V},\mathbf{f})^{\prime}(\mathbf{V},\mathbf{f})$. That is,

[TABLE]

which is equal to $(\mathbf{f}^{\prime}\mathbf{f}-\mathbf{f}^{\prime}\mathbf{V}(\mathbf{V}^{\prime}\mathbf{V})^{-1}\mathbf{V}^{\prime}\mathbf{f})\det(\mathbf{V}^{\prime}\mathbf{V})$ by the block determinant formula (e.g. [12], Theorem 13.3.8). Moreover, $\det(\mathbf{V}^{\prime}\mathbf{V})$ does not depend on $\mathbf{f}$ and $\mathbf{V}(\mathbf{V}^{\prime}\mathbf{V})^{-1}\mathbf{V}^{\prime}=\mathbf{P}(\mathcal{S}_{j})$. ∎

As a side result, we obtained that the GKM is guaranteed to produce a non-singular subset.

Clearly,

[TABLE]

therefore Theorem 1 implies that the algebraic approach based on the positive eigenvalues of the information matrix is geometrically a successive projection method: it projects the regressors to the orthogonal complement of $\mathcal{C}(\mathbf{M}(\mathcal{S}_{j}))$ and then finds the projected regressor of the largest norm. In an intuitive sense, such a regressor provides the most information not already captured by $\mathcal{S}_{j}$.

Another observation is that the GKM is, roughly stated, a limit of the RGH as $\delta$ tends to zero. First, we need a technical lemma.

Lemma 1**.**

Let $j\geq 2$, $\mathbf{F}=\mathbf{F}(\mathcal{S}_{j})$, $\delta>0$. The set on the right hand side of (4) is equal to

[TABLE]

Proof.

Let $\mathbf{M}=\mathbf{M}(\mathcal{S}_{j})$. Observe that $\lambda_{i}(\mathbf{M}+\mathbf{f}\mathbf{f}^{\prime}+\delta\mathbf{I}_{m})=\delta+\lambda_{i}(\mathbf{F}^{\prime}\mathbf{F}+\mathbf{f}\mathbf{f}^{\prime})$. The non-zero eigenvalues of $\mathbf{F}^{\prime}\mathbf{F}+\mathbf{f}\mathbf{f}^{\prime}$ coincide with the non-zero eigenvalues of $(\mathbf{F}^{\prime},\mathbf{f})^{\prime}(\mathbf{F}^{\prime},\mathbf{f})$. It follows that $\det(\mathbf{M}+\mathbf{f}\mathbf{f}^{\prime}+\delta\mathbf{I}_{m})$ is equal to $\delta^{m-j-1}\prod_{i=1}^{j+1}(\delta+\lambda_{i}((\mathbf{F}^{\prime},\mathbf{f})^{\prime}(\mathbf{F}^{\prime},\mathbf{f})))$, which can be expressed as

[TABLE]

The term $\delta^{m-j-1}$ can be ignored in the maximization, and using the formula for determinant of a block matrix, (10) is equivalent to $\mathrm{argmax}_{\mathbf{f}\in\mathcal{F}\setminus\mathcal{S}_{j}}\det(\mathbf{F}\mathbf{F}^{\prime}+\delta\mathbf{I}_{j})(\mathbf{f}^{\prime}\mathbf{f}+\delta-\mathbf{f}^{\prime}\mathbf{F}^{\prime}(\mathbf{F}\mathbf{F}^{\prime}+\delta\mathbf{I}_{j})^{-1}\mathbf{F}\mathbf{f})$. Since the first term of the expression does not depend on $\mathbf{f}$, and the additive term $+\delta$ does not affect the maximization, we obtain the desired result. ∎

Theorem 2**.**

The GKM is the limit of the RGH in the sense that for any finite set $\mathcal{G}\subset\mathbb{R}^{m}$, an $r\times m$ matrix $\mathbf{F}$, an $m\times m$ matrix $\mathbf{M}=\mathbf{F}^{\prime}\mathbf{F}$, and the orthogonal projector $\mathbf{P}$ on $\mathcal{C}(\mathbf{M})$ there is some $\delta_{*}>0$ such that

[TABLE]

for all positive $\delta<\delta_{*}$.

Proof.

Let $h:=\max_{\mathbf{f}\in\mathcal{G}}\mathbf{f}^{\prime}(\mathbf{I}_{m}-\mathbf{P})\mathbf{f}$ and let $\epsilon>0$ be such that

[TABLE]

For any matrix $\mathbf{A}$, the pseudoinverse satisfies $\mathbf{A}^{+}=\lim_{\delta\to 0^{+}}(\mathbf{A}^{\prime}\mathbf{A}+\delta\mathbf{I})^{-1}\mathbf{A}^{\prime}$ (e.g. [12], Theorem 20.7.1) and $\mathbf{F}^{\prime}(\mathbf{F}^{\prime})^{+}=\mathbf{P}$, therefore

[TABLE]

for any $\mathbf{f}$. Since $\mathcal{G}$ is finite, there is some $\delta_{*}$ such that

[TABLE]

for all $\delta<\delta_{*}$ and all $\mathbf{f}\in\mathcal{G}$. Let $\delta<\delta_{*}$ and let $\mathbf{f}_{*}\in\mathcal{G}$ be such that $\mathbf{f}_{*}\notin\mathrm{argmax}_{\mathbf{f}\in\mathcal{G}}\mathbf{f}^{\prime}(\mathbf{I}_{m}-\mathbf{P})\mathbf{f}$. Then $\mathbf{f}_{*}^{\prime}(\mathbf{I}_{m}-\mathbf{P})\mathbf{f}_{*}\leq h-\epsilon$, therefore $\mathbf{f}_{*}^{\prime}(\mathbf{I}_{m}-\mathbf{F}^{\prime}(\mathbf{F}\mathbf{F}^{\prime}+\delta\mathbf{I}_{r})^{-1}\mathbf{F})\mathbf{f}_{*}<h-\epsilon/2$. It follows that $\mathbf{f}_{*}\notin\mathrm{argmax}_{\mathbf{f}\in\mathcal{G}}\mathbf{f}^{\prime}(\mathbf{I}_{m}-\mathbf{F}^{\prime}(\mathbf{F}\mathbf{F}^{\prime}+\delta\mathbf{I}_{r})^{-1}\mathbf{F})\mathbf{f}_{*}$, i.e., $\mathbf{f}_{*}\notin\mathrm{argmax}_{\mathbf{f}\in\mathcal{G}}\phi(\mathbf{M}+\delta\mathbf{I}_{m}+\mathbf{f}\mathbf{f}^{\prime})$ by the previous lemma. ∎

If the regressors $\mathbf{f}$ of $\mathcal{G}$ are the result of iid sampling from a continuous $m$-dimensional distribution, then the set on the RHS of (11) is a singleton with probability $1$, i.e., both sides of (11) coincide. However, for a symmetric $\mathcal{G}$, the right hand side can be a large set, and the solutions of the RGH can form a strict subset of those of GKM for any $\delta>0$; see Section 4, specifically panels (a), (b), (d) in Figure 2.

The GKM is not only guaranteed to produce a non-singular saturated subset, it is guaranteed to produce a subset with efficiency at least $1/m$, as we show next. First, a technical lemma.

Lemma 2**.**

Let $\mathbf{f}_{1},\ldots,\mathbf{f}_{m}\in\mathbb{R}^{m}$. Then, $\det(\sum_{j}\mathbf{f}_{j}\mathbf{f}_{j}^{\prime})=\det(\sum_{j}\tilde{\mathbf{f}}_{j}\tilde{\mathbf{f}}_{j}^{\prime})$, where $\tilde{\mathbf{f}}_{1}=\mathbf{f}_{1}$ and $\tilde{\mathbf{f}}_{j}$ for $j>1$ can be defined by any of the three equivalent formulas:

[TABLE]

In the expressions above $\mathbf{P}_{1:(j-1)}$ is the projection matrix to $\mathcal{C}(\mathbf{f}_{1},\ldots,\mathbf{f}_{j-1})$ and $\tilde{\mathbf{P}}_{i}$ is the projection matrix to $\mathcal{C}(\tilde{\mathbf{f}}_{i})$.

Proof.

From formula (13), it is clear that the transformation from $\mathbf{F}=(\mathbf{f}_{1},\ldots,\mathbf{f}_{m})$ to $\tilde{\mathbf{F}}=(\tilde{\mathbf{f}}_{1},\ldots,\tilde{\mathbf{f}}_{m})$ is in fact a Gram-Schmidt orthogonalization, where the obtained orthogonal vectors $\tilde{\mathbf{f}}_{j}$ are not normalized. That is, $\mathbf{F}=\tilde{\mathbf{F}}\mathbf{R}$ for some upper triangular matrix $\mathbf{R}$ with ones on the diagonal. It follows that $\det(\mathbf{F})=\det(\tilde{\mathbf{F}})$, which yields the desired result $\det(\mathbf{F}\mathbf{F}^{\prime})=\det(\tilde{\mathbf{F}}\tilde{\mathbf{F}}^{\prime})$.

Since $\tilde{\mathbf{f}}_{1},\ldots,\tilde{\mathbf{f}}_{m}$ are orthogonal, we have $\tilde{\mathbf{P}}_{i}\tilde{\mathbf{P}}_{k}=\mathbf{0}$ for any $i\neq k$, and hence (14) is obtained. Formula (12) holds because $\mathcal{C}(\mathbf{f}_{1},\ldots,\mathbf{f}_{j-1})=\mathcal{C}(\tilde{\mathbf{f}}_{1},\ldots,\tilde{\mathbf{f}}_{j-1})$. ∎

Theorem 3**.**

Let $\mathcal{S}_{GK}$ be the saturated subset of $\mathcal{F}$ produced by the GKM. Then $\mathrm{eff}(\mathcal{S}_{GK})\geq\frac{1}{m}$.

Proof.

Let $\mathcal{S}_{GK}=\{\mathbf{f}_{1},\ldots,\mathbf{f}_{m}\}$ and define $\tilde{\mathcal{S}}=\{\tilde{\mathbf{f}}_{1},\ldots,\tilde{\mathbf{f}}_{m}\}$, where the $\tilde{\mathbf{f}}_{j}$’s are given by Lemma 2. Then $\phi(\mathcal{S}_{GK})=\phi(\tilde{\mathcal{S}})$. Let us also work with “normalized” information matrix $\mathbf{N}(\mathcal{S})=\sum_{\mathbf{f}\in\mathcal{S}}\mathbf{f}\mathbf{f}^{\prime}/s$, where $s$ is the size of $\mathcal{S}$. Without loss of generality, let us rotate the coordinate axes of $\mathbb{R}^{m}$, so that $\tilde{\mathbf{f}}_{j}=\|\tilde{\mathbf{f}}_{j}\|\mathbf{e}_{j}$, which is possible because $\tilde{\mathbf{f}}_{1},\ldots,\tilde{\mathbf{f}}_{m}$ are orthogonal. Then, $\det(\sum_{i}\tilde{\mathbf{f}}_{i}\tilde{\mathbf{f}}_{i}^{\prime})=\prod_{i}\|\tilde{\mathbf{f}}_{i}\|^{2}$. Moreover, the set $\mathcal{S}_{B}$ that consists of all points of the form $(\pm\|\tilde{\mathbf{f}}_{1}\|,\ldots,\pm\|\tilde{\mathbf{f}}_{m}\|)^{\prime}$ satisfies $\det(\mathbf{N}(\mathcal{S}_{B}))=\prod_{i}\|\tilde{\mathbf{f}}_{i}\|^{2}$ and $\mathbf{N}(\mathcal{S}_{B})\succeq\mathbf{N}(\mathcal{S})$ for any $\mathcal{S}\subseteq\mathcal{F}$ because all the elements of $\mathcal{F}$ lie in $\mathrm{conv}(\mathcal{S}_{B})$. In particular, $\mathbf{N}(\mathcal{S}_{B})\succeq\mathbf{N}(\mathcal{S}^{*})$, where $\mathcal{S}^{*}$ is the $D$-optimal subset of size $m$. It follows that $\det(\mathbf{M}(\mathcal{S}^{*}))=m^{m}\det(\mathbf{N}(\mathcal{S}^{*}))\leq m^{m}\det(\mathbf{N}(\mathcal{S}_{B}))=m^{m}\prod_{i}\|\tilde{\mathbf{f}}_{i}\|^{2}$. Therefore,

[TABLE]

∎

As a side result, we obtain an explicit formula for the $D$-optimality value of $\mathcal{S}_{GK}$: $\phi(\mathcal{S}_{GK})=\prod_{i}\|\tilde{\mathbf{f}}_{i}\|^{2/m}$, where the projected regressors $\tilde{\mathbf{f}}_{i}$ are constructed during the GKM process.

A naïve implementation of the GKM using (8) is relatively time consuming as it requires the calculations of projection matrices $\mathbf{P}(\mathcal{S}_{j})$ and the projections of the $n$ regressors. However, (8) depends only on the norms of the projected regressors $\tilde{\mathbf{f}}$ and formula (14) shows that the projected regressors can be calculated successively by applying projections to the one-dimensional subspaces orthogonal to $\tilde{\mathbf{f}}_{j}$. Therefore, the algorithm can “forget” the original regressors, which allows for the following efficient implementation of the GKM:

**Efficient GKM

**1. Set $\textit{ind}\leftarrow\emptyset$

2. $j\leftarrow\arg\max_{i=1,\ldots,n}\lVert\mathbf{f}_{i}\rVert$

3. Set $\textit{ind}\leftarrow\textit{ind}\cup\{j\}$

4. If $j<m$, update $\mathbf{f}_{i}\leftarrow\mathbf{f}_{i}-[\mathbf{f}_{j}^{\prime}\mathbf{f}_{i}/(\mathbf{f}_{j}^{\prime}\mathbf{f}_{j})]\mathbf{f}_{j}$, set $j\leftarrow j+1$ and return to 2

5. Output ind

The efficient formulation of the GKM has asymptotic complexity $O(nm^{2})$. Its implementation in the programming language of R is in Appendix.

The GKM is non-random which is, in some situations, a drawback.555For instance, if we use the saturated subsets for initiation of deterministic heuristics, a randomized initiation method is crucial for the approach of multiple restarts. Galil and Kiefer were aware of this problem and they proposed a randomisation based on first selecting some points at random and only then starting to apply the GKM. However, such an approach may lead to singular $m$-point subsets, if the first random points were chosen “unluckily”.

Nevertheless, the GKM can be also randomised as follows: in step (S), instead of selecting regressor $\mathbf{f}_{(j)}$ that maximizes $v(\mathbf{f}):=\mathbf{f}^{\prime}(\mathbf{I}_{m}-\mathbf{P}(\mathbf{M}(\mathcal{S}_{j}))\mathbf{f}$, we may select $\mathbf{f}_{(j)}$ at random, with probabilities of the particular regressors $\mathbf{f}$ proportional to $v^{\alpha}(\mathbf{f})$, where $\alpha>0$; the higher $\alpha$ is, the less random is the procedure.

Because the randomised GKM never chooses regressors that lie in $\mathrm{Span}(\mathbf{f}_{1},\ldots,\mathbf{f}_{j-1})$, it also arrives at a non-singular subset if such a subset exists by the same argument as in the non-random case.

2.4 The Kumar-Yildirim method (KYM)

Originally, the Kumar-Yildirim method was proposed in [15] to find initial points for computing the $\mathrm{MVEE}$, which need not be centred at the origin. In each iteration, the method chooses $2$ points $\mathbf{f}_{j_{1}},\mathbf{f}_{j_{2}}\in\mathcal{F}$ maximizing and minimizing the scalar product $\mathbf{f}^{\prime}\mathbf{b}$ for a randomly chosen $\mathbf{b}$. Thus, the output is a set of $2m$ points.

But the $D$-optimal approximate design problem is equivalent to computing the $\mathrm{MVEE}_{0}$ centred at the origin and containing $\mathbf{f}_{1},\ldots,\mathbf{f}_{n}$. For such a problem, one can apply the original Kumar-Yildirim method with formally added points $-\mathbf{f}_{1},\ldots,-\mathbf{f}_{n}$. Then, the method selects in each step both $\mathbf{f}_{j}$ and $-\mathbf{f}_{j}$ for some $j$, which is redundant. It follows that it is enough to choose one of the two points.

Hence, the modified KYM can be stated as a special case of the general greedy heuristic where step (S) is

[TABLE]

and $\mathbf{b}$ is a random direction from $\mathcal{C}^{\perp}(\mathcal{S}_{j})$.

It is not difficult to see that the output of the KYM is a non-singular saturated subset. However, we can formulate a stronger statement as follows.

Theorem 4**.**

Let $\mathcal{S}_{KY}$ be the saturated subset of $\mathcal{F}$ produced by the KYM. Then

[TABLE]

Proof.

The authors of [15] obtain the efficiency bound for the output of their algorithm by applying the results of the paper [3]. There, it is shown that $|\det(\mathbf{z}_{1},\ldots,\mathbf{z}_{m})|/m!\leq\mathrm{vol}(K)\leq|\det(\mathbf{z}_{1},\ldots,\mathbf{z}_{m})|$, where each $\mathbf{z}_{j}$ is equal to the difference $\mathbf{f}_{j_{1}}-\mathbf{f}_{j_{2}}$ as in the KYM, but for a general convex body $K$. By considering the convex body $K=\mathrm{conv}(\pm\mathbf{f}_{1},\ldots,\pm\mathbf{f}_{n})$666$K$ is called the Elfving set in the optimal experimental design literature, we obtain for $\mathcal{S}_{KY}$ that

[TABLE]

where $2\mathbf{F}(\mathcal{S}_{KY})$ appears in the inequality because in each iteration of the algorithm in [3], one chooses $\mathbf{z}_{j}=\mathbf{f}_{(j)}-(-\mathbf{f}_{(j)})=2\mathbf{f}_{(j)}$ for the $\mathbf{f}_{(j)}$ maximizing $|\mathbf{f}_{i}^{\prime}\mathbf{b}|$.

Let $\mathrm{MVEE}^{*}=\mathrm{MVEE}(K)$. In [13], John showed that $K\supseteq\mathrm{MVEE}(K)/\sqrt{m}$ if $K$ is centred at the origin; hence $\mathrm{vol}(K)\geq\mathrm{vol}(\mathrm{MVEE}^{*})/m^{m/2}$ in our settings. Combining the previous results, we obtain that

[TABLE]

From the results mentioned, e.g., in [15] it follows that

[TABLE]

where $\mathcal{S}^{*}$ is the $D$-optimal saturated subset and $\eta=\pi^{m/2}/\Gamma(1+m/2)$ is the volume of the unit ball. Then,

[TABLE]

and

[TABLE]

∎

For a large $m$, the Stirling approximation of $\Gamma$ yields

[TABLE]

Note that by a direct application of the efficiency limit provided by Lemma 3.1 of [15], one would obtain a weaker bound $\mathrm{eff}(\mathcal{S}_{KY})\geq m^{-3}$ than that in Theorem 4. The bound in Theorem 4 is tighter because Kumar and Yildirim deal with general sets, and the centrally symmetric sets considered in this paper allow for stronger results.

An efficient implementation of the KYM in the programming language of R is given in Appendix. This implementation has asymptotic complexity $O(nm^{2})$, i.e., the same as the GKM implementation, although non-asymptotically the KYM tends to be faster than GKM (unless $m\gtrapprox n/2$).

Interestingly, the main idea of the KYM is similar to the projective interpretation of the GKM: both procedures seek to find a “good” new regressor in the space $\mathcal{C}(\mathbf{M}(\mathcal{S}))^{\perp}$ not yet spanned by the chosen regressors. However, the respective approaches for measuring what is a good regressor differ.

Note that both GKM and KYM are invariant under rotation, and, unlike the RGH, they are also invariant with respect to a centred dilation $\mathbf{f}\to c\mathbf{f}$, $c>0$, i.e., invariant under the common change of the units of measurement, if the regressors correspond to explanatory variables.

3 Variants and extensions

3.1 Subsets of sizes other than $m$

The greedy heuristics for saturated subsets can be easily modified to generate efficient non-saturated subsets, i.e., subsets of sizes $s$ other than $m$.

The problem of obtaining efficient supersaturated experimental designs (e.g., [8]) is very similar to problem of constructing $D$-efficient subsets with $s<m$. In fact, the main idea of the Galil-Kiefer initiation for $D$-optimality, i.e. that $\det(\mathbf{F}(\mathcal{S})\mathbf{F}^{\prime}(\mathcal{S}))$ can be maximized instead of $\det(\mathbf{F}^{\prime}(\mathcal{S})\mathbf{F}(\mathcal{S}))$, is the same as the approach to $D$-optimal supersaturated designs in [14]. Note, however, that the motivation of Jones and Majumdar in [14] for working with $\mathbf{F}(\mathcal{S})\mathbf{F}^{\prime}(\mathcal{S})$ arises from a minimum bias estimator of the model parameters, and that they provide analytical results on optimal supersaturated designs in main-effects model, not design algorithms.

Based on a greedy heuristic for saturated subsets, we can formulate a general scheme for a greedy heuristic for arbitrary size subsets $s$. Simply stated, if $s<m$ we can pre-maturely stop the greedy heuristics and if $s>m$, we can repeatedly perform the heuristic for a subset of size $m$ (or less than $m$ in the last run) and after each run of the method remove the $m$ obtained regressors from $\mathcal{F}$. We can generalize this idea such that at each selection step, we select not a single regressor but a batch of $b$ regressors, which can be computationally advantageous. More precisely:

**Scheme of greedy heuristics for subsets of size $s$

**1. Set $j\leftarrow 0$ and $\mathcal{S}_{0}\leftarrow\emptyset$

2. Set $a\leftarrow\min(s-j,b)$

3. Select an $a$-element subset $\mathcal{F}_{j}$ of $\mathcal{F}\setminus\mathcal{S}_{j}$   **(S)

**4. Set $\mathcal{S}_{j+a}\leftarrow\mathcal{S}_{j}\cup\mathcal{F}_{j}$

5. If $j<s$ set $j\leftarrow j+a$ and return to 2.

6. Output $\mathcal{S}_{s}$.

Of course here we will also need to keep and properly update some parameters which regulate step (S). Special case $b=1$, $s=m$ provides the scheme of the heuristics for saturated subsets, but note that the IBOSS mentioned in the introduction and the original KYM are also special cases of this general scheme. The multitude of various specifications of the above scheme and the comparison of the resulting algorithms is left for further research.

3.2 A pre-selection strategy for a big $n$

Let $\mathcal{R}:=\{\mathbf{g}_{1},\ldots,\mathbf{g}_{km}\}$ be a uniformly randomly selected $km$-element subset of $\mathcal{F}$, where $1\leq k\leq n/m$. Let us estimate the probability that $\mathcal{R}$ is non-singular. Let $\mathbf{h}_{1},\ldots,\mathbf{h}_{km}$ be uniformly randomly chosen from $\mathcal{F}$ with replacement. For $t=0,\ldots,k-1$ let $A_{t}$ be the event that $[\mathbf{h}_{tm+1},\ldots,\mathbf{h}_{(t+1)m}]$ is a singular $m\times m$ matrix. Clearly, $A_{0},\ldots,A_{k-1}$ are independent, therefore

[TABLE]

where $p$ is the probability of each of $A_{t}$. That is, the probability that a $km$-element uniformly selected random subset of $\mathcal{F}$ is singular vanishes at least exponentially with respect to $k$.

At the same time, the probability $p$ is itself not too close to $1$ for most of sets $\mathcal{F}$ used in applications. Even for such a symmetric set as $\mathcal{F}=\{-1,1\}^{m}$ from Example 1 it is possible to show that $p\leq(3/4+O(1))^{m}$; cf. [27].

This motivates the following strategy to obtain an efficient saturated subset of $\mathcal{F}$:

Randomly pre-select $\tilde{\mathcal{F}}$ of size $km$ from $\mathcal{F}$. 2. 2.

Use either GKM or KYM to construct a saturated subset of $\tilde{\mathcal{F}}$.

Note that the size of $\tilde{\mathcal{F}}$ does not depend on $n$, i.e., it can be small also for an enormous set $\mathcal{F}$. Even if $k$ is only moderately large, the probability that the resulting saturated subset is singular can be negligibly small. For instance, if $k=50$ and $p=3/4$, the probability of a singular resulting subset is less than $10^{-6}$.

We will show in Section 4 that the pre-selection strategy is a good compromise between the speed and efficiency.

3.3 Multi-run approach

Let H be a randomized heuristic. Since H can generate subsets of various efficiencies, a natural idea is to run H multiple times and then select the best subset obtained. Note that in this respect the randomness inherent in the result can be a positive, not a negative feature of an algorithm.

Denote $Y_{i}$, $i=1,2,\ldots$, the iid random variables representing the efficiencies of the subsets obtained by independent runs of H. Let $0<\mathrm{eff}\lesssim 1$ be such that $P[Y_{1}\geq\mathrm{eff}]>0$. Then it can be shown that we need roughly $\ln(2)/P[Y_{1}\geq\mathrm{eff}]$ runs of H for the median of $\max_{1\leq i\leq N}Y_{i}$ to be $\mathrm{eff}$. In this sense, two heuristics are similarly good, if they have similar ratios $t^{-1}P[Y_{1}\geq\mathrm{eff}]$, where $t$ is the time of a single run of H plus the time needed to compute the quality of the resulting subset.

The performance of the multi-run approach is therefore given by the speed of H as well as the tails of the distribution of $Y_{1}$, and both quantities may be difficult to assess. However, numerical simulations for the randomized heuristics that we consider in this paper suggest that the variability of $Y_{1}$ tends to be rather small, therefore $P[Y_{1}\geq\mathrm{eff}]$ is mostly determined by the median of $Y_{1}$. Hence, in the assessment of the performance of a heuristic, the median efficiency of the subsets resulting from a single run is the more important factor than the speed. See Section 4 for numerical results.

3.4 Eigenvalue-based criteria

We say that a criterion $\phi:2^{\mathcal{F}}\to\mathbb{R}$ is eigenvalue-based if it depends only on the eigenvalues of $\mathbf{M}(\mathcal{S})$. For instance, $D$-optimality can be expressed as $\phi(\mathcal{S})=\left(\prod_{i}\lambda_{i}(\mathbf{M}(\mathcal{S}))\right)^{1/m}$. Besides $D$-optimality, the second most important eigenvalue-based criterion is $A$-optimality defined as $\phi_{A}(\mathcal{S})=\left(\sum_{i}\lambda_{i}^{-1}(\mathbf{M}(\mathcal{S}))\right)^{-1}$ for a non-singular $\mathbf{M}(\mathcal{S})$ and $\phi_{A}(\mathcal{S})=0$ for a singular $\mathbf{M}(\mathcal{S})$ (e.g., [22], Chapter 6).

Clearly, the basic idea of the Galil-Kiefer method can be used for any eigenvalue-based criterion – we can simply specify (S) as selecting $\mathbf{f}_{(j)}$ that maximizes an appropriate function of the non-zero eigenvalues of $\mathbf{M}(\mathcal{S}_{j})+\mathbf{f}\mathbf{f}^{\prime}$. In some cases, like $D$-optimality, such maximization can be performed much more easily than via numerical evaluation of the eigenvalues. For an example, we can utilize the following theorem for $A$-optimality:

Theorem 5**.**

Let $\mathbf{M}(\mathcal{S})$ be of rank $r$. Then,

[TABLE]

Proof.

The positive eigenvalues of $\mathbf{M}(\mathcal{S})+\mathbf{f}\mathbf{f}^{\prime}=(\mathbf{V},\mathbf{f})(\mathbf{V},\mathbf{f})^{\prime}$ are equal to the positive eigenvalues of $(\mathbf{V},\mathbf{f})^{\prime}(\mathbf{V},\mathbf{f})$ with the notation as in the proof of Theorem 1. Hence

[TABLE]

which is equal to $\mathrm{tr}(\mathbf{V}^{\prime}\mathbf{V})^{-1}+s^{-1}\mathbf{f}^{\prime}\mathbf{V}(\mathbf{V}^{\prime}\mathbf{V})^{-2}\mathbf{V}^{\prime}\mathbf{f}+s^{-1}$, where $s=\mathbf{f}^{\prime}(\mathbf{I}_{m}-\mathbf{V}(\mathbf{V}^{\prime}\mathbf{V})^{-1}\mathbf{V}^{\prime})\mathbf{f}=\mathbf{f}^{\prime}(\mathbf{I}_{m}-\mathbf{P}(\mathbf{M}(\mathcal{S})))\mathbf{f}$, based on the formula for the inverse of a block matrix ([12], Theorem 8.5.11). Observing that $\mathbf{V}(\mathbf{V}^{\prime}\mathbf{V})^{-2}\mathbf{V}^{\prime}=\mathbf{M}^{+}$ yields that the studied maximization problem is equivalent to the minimization of

[TABLE]

Because $\mathrm{tr}(\mathbf{V}^{\prime}\mathbf{V})^{-1}$ does not depend on $\mathbf{f}$, the minimization of (16) is equivalent to the maximization of (5). ∎

That is, the GKM for $A$-optimality is very similar to the GKM for $D$-optimality, but $A$-optimality requires an additional “correction term” to be added to (8).

It is unclear how can the KYM be tailored to $A$-optimality and other eigenvalue-based criteria. However, the efficiency of the $D$-optimal approximate design with respect to a very general class of eigenvalue-based criteria is bounded from below by $1/m$ ([10]) which suggests that the KYM (as well as the GKM for $D$-optimality) will produce reasonably efficient subsets even in the case of other criteria than $D$-optimality. This is also supported by our numerical experience. We leave the development of the heuristics for general eigenvalue-based criteria to further research.

4 Numerical comparisons

To provide a brief but comprehensive comparison of the performances of the greedy methods we generate the sets $\mathcal{F}$ at random, as well as deterministically:

We first generate a covariance matrix $\Sigma$ from the Wishart distribution $\mathcal{W}_{d}(\mathbf{I}_{d},d)$ and then $\mathbf{x}_{i}$ independently from $\mathcal{N}_{d}(\mathbf{0}_{d},\Sigma)$. We set $\mathcal{F}=\{(\mathbf{x}^{\prime}_{1},1)^{\prime},\ldots,(\mathbf{x}^{\prime}_{n},1)^{\prime}\}\subset\mathbb{R}^{m}$, $m=d+1$. 2. 2.

We set $\mathcal{F}=\{-1,1\}^{m}$.

The efficiencies of the saturated subsets $\mathcal{S}$ are calculated as $\phi(\mathcal{S})/\phi^{*}_{m}$, where $\phi^{*}_{m}$ is computed via the REX algorithm suggested in [11]. We used a computer with a 64-bit Windows 10 system running an Intel Core i7-5500U CPU processor at $2.40$ GHz with $8$ GB of RAM. For the RGH, we used the regularization parameter $\delta=10^{-4}$.

The results are exhibited in Figures 1 and 2. The numerical observations can be summarized as follows:

• •

Both RGH and GKM tend to provide highly efficient subsets. For a random $\mathcal{F}$ the results of GKM and RGH usually coincide and for $\mathcal{F}=\{-1,1\}^{m}$, RGH sometimes gives slightly more efficient subsets than GKM. On the other hand, GKM tends to be as fast as or faster than RGH.

• •

The KYM is somewhat faster but less efficient than GKM.

• •

As one might expect, the random initiation (RND) is the fastest and the least efficient method. The leverage-weighted random sampling (RNDl), performs badly in all analysed cases. For $\mathcal{F}=\{-1,1\}^{m}$ both RND and RNDl occasionally fail to the extent that they provide singular subsets.

• •

GKM and KYM applied to a random subsample of size $50m$ (cf. Subsection 3.2), denoted by GKMf and KYMf in the figures, are both very fast and efficient. In some cases, the speed of GKMf and KYMf is only slightly smaller than the speed of the random sub-sampling methods, yet, with a few exceptions, the efficiency is the same or almost the same as that of GKM and KYM applied to the full set $\mathcal{F}$. Moreover, the computation times and the resulting efficiencies are relatively stable under changes of $n$ and $m$.

To compare the degree of the efficiency improvement using the multi-run approach described in Subsection 3.3, we computed the “time profiles” of the increase in efficiency of the best subset found by a repeated independent application of the individual heuristics. Evidently, for all procedures the efficiency increases very slowly (notice that the horizontal axis denotes the logarithm of the time), i.e., for the performance of the multi-run method the efficiency of a single-run is crucial. Naturally, the numerical differences between the studied methods depend on a multitude of factors, such as the actual implementation used and the set $\mathcal{F}$ itself. Nevertheless, the exhibited figures are representative of the general behaviour of the heuristics that we observed.

5 Discussion

We showed that a good approach to the construction of $D$-efficient saturated subsets is to use the efficient version of the Galil-Kiefer method or the modified Kumar-Yildirim method. Nevertheless, the regularized greedy method and the random sampling can also be reasonable choices for some applications (the regularized greedy method if we are more concerned about the efficiency than the speed and the random sampling in the opposite situation). These heuristics can be directly generalized for the construction of non-saturated subsets. For a large size of the initial set, the pre-sampling strategy can be recommended. To minimize the risk of a singular subset, we can use the multi-run approach. Finally, the regularized greedy method and the closely related Galil-Kiefer method can be modified for the use with a different criterion than $D$-optimality.

Appendix - R Scripts

Regularized greedy heuristic

function(Fx, del=1e-4) {

Fx is the nxm matrix of regressors

m <- ncol(Fx); one <- rep(1, m) M <- del * diag(m) j <- which.max((Fx^2) %% one) ind <- rep(0, m); ind[1] <- j for(i in 2:m) { M <- M + tcrossprod(Fx[j, ]) j <- which.max(((Fx %% t(chol(solve(M))))^2) %*% one) ind[i] <- j } return(ind) }

Galil-Kiefer method

function(Fx) {

Fx is the nxm matrix of regressors

m <- ncol(Fx) v2 <- (Fx^2) %% rep(1, m) j <- which.max(v2) ind <- rep(0, m); ind[1] <- j for(i in 2:m) { scp <- Fx %% Fx[j, ] Fx <- Fx - scp %*% t(Fx[j, ]) / v2[j] v2 <- v2 - scp^2 / v2[j] j <- which.max(v2); ind[i] <- j } return(ind) }

Kumar-Yildirim method

function(Fx) {

Fx is the nxm matrix of regressors

m <- ncol(Fx); P <- diag(m) j <- which.max(abs(Fx %% rnorm(m))) ind <- rep(0, m); ind[1] <- j for(i in 2:m) { fx <- P %% Fx[j,] P <- P - tcrossprod(fx) / sum(fx^2) j <- which.max(abs(Fx %% (P %% rnorm(m)))) ind[i] <- j } return(ind) }