Adaptive Iterative Hessian Sketch via A-Optimal Subsampling

TL;DR

This paper introduces an improved deterministic iterative Hessian sketch method using A-optimal subsampling, enhancing efficiency and accuracy for large-scale least squares problems through novel initialization, preconditioning, and adaptive step sizing.

Contribution

The paper proposes a deterministic A-optimal subsampling approach to enhance the iterative Hessian sketch for large-scale least squares, including new initialization, preconditioning, and step size methods.

Findings

Outperforms existing accelerated IHS methods in experiments.

Provides a deterministic alternative to randomized sketching.

Demonstrates improved computational efficiency and accuracy.

Abstract

Iterative Hessian sketch (IHS) is an effective sketching method for modeling large-scale data. It was originally proposed by Pilanci and Wainwright (2016; JMLR) based on randomized sketching matrices. However, it is computationally intensive due to the iterative sketch process. In this paper, we analyze the IHS algorithm under the unconstrained least squares problem setting, then propose a deterministic approach for improving IHS via A-optimal subsampling. Our contributions are three-fold: (1) a good initial estimator based on the A-optimal design is suggested; (2) a novel ridged preconditioner is developed for repeated sketching; and (3) an exact line search method is proposed for determining the optimal step length adaptively. Extensive experimental results demonstrate that our proposed A-optimal IHS algorithm outperforms the existing accelerated IHS methods.

Figures

Tables

Table 1. Table 1: Averaged Δ Δ \Delta -values of the ridged preconditioner 𝑴 ​ ( λ , 𝜹 ) = n / m ​ ∑ i = 1 n δ i ​ 𝒙 i ​ 𝒙 i T + λ ​ 𝑰 d 𝑴 𝜆 𝜹 𝑛 𝑚 superscript subscript 𝑖 1 𝑛 subscript 𝛿 𝑖 subscript 𝒙 𝑖 superscript subscript 𝒙 𝑖 𝑇 𝜆 subscript 𝑰 𝑑 \bm{M}(\lambda,\bm{\delta})=n/m\sum_{i=1}^{n}\delta_{i}\bm{x}_{i}\bm{x}_{i}^{T}+\lambda\bm{I}_{d} and SRHT. The λ ~ ~ 𝜆 \tilde{\lambda} is chosen to be 0.1 ​ ∑ i = 1 n ‖ 𝒙 i ‖ 2 2 0.1 superscript subscript 𝑖 1 𝑛 subscript superscript norm subscript 𝒙 𝑖 2 2 0.1\sum_{i=1}^{n}\|\bm{x}_{i}\|^{2}_{2} for the normal distribution and 0.4 ​ ∑ i = 1 n ‖ 𝒙 i ‖ 2 2 0.4 superscript subscript 𝑖 1 𝑛 subscript superscript norm subscript 𝒙 𝑖 2 2 0.4\sum_{i=1}^{n}\|\bm{x}_{i}\|^{2}_{2} for other cases.

	Normal	Log-Normal	$t_2$	Mixture
	$d=50$
$\lambda =0$	-4.56	0.39	0.76	0.42
$\lambda =\tilde{\lambda }$	0.87	0.76	0.89	0.79
SRHT	0.55	0.52	0.75	0.70
	$d=100$
$\lambda =0$	-7.62	-0.41	0.63	-0.01
$\lambda =\tilde{\lambda }$	0.83	0.73	0.90	0.82
SRHT	-0.04	0.05	0.46	0.32

Table 2. Table 2: The averaged time for each approach achieving the precision 10 − 10 superscript 10 10 10^{-10} measured by ‖ 𝜷 ^ − 𝜷 ^ LS ‖ 2 subscript norm ^ 𝜷 superscript ^ 𝜷 LS 2 \|\hat{\bm{\beta}}-\hat{\bm{\beta}}^{\rm LS}\|_{2} , and the required numbers of iterations.

Method	Normal		Log-Normal		$t_2$		Mixture
	Time(s)	Iter	Time(s)	Iter	Time(s)	Iter	Time(s)	Iter
	$d=50$
IHS	15.28	18.30	14.83	18.25	15.63	18.43	15.21	18.38
acc-IHS	5.63	25.95	5.20	25.87	5.57	26.17	5.88	26.05
pwGradient	8.54	46.99	8.19	47.40	8.98	47.34	9.22	47.67
Aopt-IHS	3.01	10.27	3.90	14.97	3.75	12.65	4.74	17.39
	$d=100$
IHS	53.55	27.29	54.89	27.16	48.17	27.51	43.31	27.40
acc-IHS	19.59	40.64	19.59	40.51	17.11	41.01	15.42	40.92
pwGradient	—	diverge	—	diverge	—	diverge	—	diverge
Aopt-IHS	11.28	19.44	11.16	19.07	11.32	22.78	9.48	20.45

Equations

$${\hat{\mathbf{\beta }}}^{\mathrm{LS}}$$

$${\hat{\mathbf{\beta }}}^{\mathrm{CS}}$$

$${\hat{\mathbf{\beta }}}^{\mathrm{HS}}$$

$${\hat{\mathbf{\beta }}}_t^{\mathrm{IHS}}$$

$$\mathbf{S}=\sqrt{\frac{n}{m}}\mathbf{RHD},$$

$${\hat{\mathbf{\beta }}}_t$$

$$=$$

$${\hat{\mathbf{\beta }}}_t=\prod _{i=1}^t({\mathbf{I}}_d-{\mathbf{A}}_i){\hat{\mathbf{\beta }}}_0+\left[{\mathbf{I}}_d-\prod _{i=1}^t({\mathbf{I}}_d-{\mathbf{A}}_i)\right]{\hat{\mathbf{\beta }}}^{\mathrm{LS}}.$$

$$\Vert {\hat{\mathbf{\beta }}}_t-{\hat{\mathbf{\beta }}}^{\mathrm{LS}}{\Vert }_2\le {\left(\frac{\max \{{\epsilon }_1,{\epsilon }_2\}}{1-{\epsilon }_1}\right)}^t\Vert {\hat{\mathbf{\beta }}}_0-{\hat{\mathbf{\beta }}}^{\mathrm{LS}}{\Vert }_2,$$

$$1-{\epsilon }_1\le \frac{{\mathbf{a}}^T{\mathbf{M}}_t\mathbf{a}}{{\mathbf{a}}^T{\mathbf{X}}^T\mathbf{X}\mathbf{a}}\le 1+{\epsilon }_2,t=1,2,\dots$$

$${\hat{\mathbf{\beta }}}_t-{\hat{\mathbf{\beta }}}^{\mathrm{LS}}=\prod _{i=1}^t({\mathbf{I}}_d-{\mathbf{A}}_i)({\hat{\mathbf{\beta }}}_0-{\hat{\mathbf{\beta }}}^{\mathrm{LS}}).$$

$$\Vert {\hat{\mathbf{\beta }}}_t-{\hat{\mathbf{\beta }}}^{\mathrm{LS}}{\Vert }_2\le \prod _{i=1}^t\Vert {\mathbf{I}}_d-{\mathbf{A}}_i{\Vert }_2\Vert {\hat{\mathbf{\beta }}}_0-{\hat{\mathbf{\beta }}}^{\mathrm{LS}}{\Vert }_2.$$

$$\Vert {\mathbf{I}}_d-{\mathbf{A}}_i{\Vert }_2$$

$$\Vert {\mathbf{M}}_i^{-1}\mathbf{Q}{\Vert }_2=\Vert {\mathbf{Q}}^{1/2}{\mathbf{M}}_i^{-1}{\mathbf{Q}}^{1/2}{\Vert }_2,$$

$$\Vert {\mathbf{I}}_d-{\mathbf{Q}}^{-1}{\mathbf{M}}_i{\Vert }_2=\Vert {\mathbf{I}}_d-{\mathbf{Q}}^{-1/2}{\mathbf{M}}_i{\mathbf{Q}}^{-1/2}{\Vert }_2.$$

$$1-{\epsilon }_1\le \frac{{\mathbf{b}}^T{\mathbf{Q}}^{-1/2}{\mathbf{M}}_i{\mathbf{Q}}^{-1/2}\mathbf{b}}{{\mathbf{b}}^T\mathbf{b}}\le 1+{\epsilon }_2,$$

$$\Vert {\mathbf{M}}_i^{-1}\mathbf{Q}{\Vert }_2$$

$$\Vert {\mathbf{I}}_d-{\mathbf{Q}}^{-1}{\mathbf{M}}_i{\Vert }_2$$

$$\Vert {\mathbf{I}}_d-{\mathbf{A}}_i{\Vert }_2$$

$$\Vert {\hat{\mathbf{\beta }}}_t-{\hat{\mathbf{\beta }}}^{\mathrm{LS}}{\Vert }_2\le {\left(\frac{\max \{{\epsilon }_1,{\epsilon }_2\}}{1-{\epsilon }_1}\right)}^t\Vert {\hat{\mathbf{\beta }}}_0-{\hat{\mathbf{\beta }}}^{\mathrm{LS}}{\Vert }_2.$$

$${\hat{\mathbf{\eta }}}_t={\hat{\mathbf{\eta }}}_{t-1}+{\mathbf{M}}^{-1/2}{\mathbf{X}}^T(\mathbf{y}-\mathbf{X}{\mathbf{M}}^{-1/2}{\hat{\mathbf{\eta }}}_{t-1}),$$

$$\tilde{f}(\mathbf{\eta })$$

$${\hat{\mathbf{\beta }}}^{\mathrm{CS}}(\mathbf{\delta })$$

$$\text{mbox}cov({\hat{\mathbf{\beta }}}^{\mathrm{CS}}(\mathbf{\delta }))$$

$$\underset{\mathbf{\delta }\in \{0,1{\}}^n}{\min }\mathrm{Tr}[{\mathbf{M}}^{-1}(\mathbf{\delta })],\text{mbox}subjectto\sum _{i=1}^n{\delta }_i=m.$$

$$\mathrm{Tr}[{\mathbf{M}}^{-1}(\mathbf{\delta })]\le \frac{1}{{\lambda }_{\min }(\mathbf{Q})}\left[d+\frac{\kappa (\mathbf{Q})}{C}\sum _{i=1}^n(1-{\delta }_i)\Vert {\mathbf{x}}_i{\Vert }_2^2\right].$$

$$\mathrm{Tr}[{\mathbf{M}}^{-1}(\mathbf{\delta })]$$

$$({\mathbf{I}}_d-{\mathbf{U}}^T\bar{\mathbf{W}}\mathbf{U}{)}^{-1}={\mathbf{I}}_d+\sum _{k=1}^{\infty }({\mathbf{U}}^T\bar{\mathbf{W}}\mathbf{U}{)}^k.$$

$$\mathrm{Tr}\left[({\mathbf{I}}_d-{\mathbf{U}}^T\bar{\mathbf{W}}\mathbf{U}{)}^{-1}\right]$$

$$1-\Vert {\mathbf{U}}^T\bar{\mathbf{W}}\mathbf{U}{\Vert }_2$$

$${\lambda }_{\min }({\mathbf{Q}}^{-1/2}\mathbf{M}(\mathbf{\delta }){\mathbf{Q}}^{-1/2})$$

$$\mathrm{Tr}({\mathbf{U}}^T\bar{\mathbf{W}}\mathbf{U})$$

Full text

1 Introduction

Consider the linear model $\bm{y}=\bm{X}\bm{\beta}+\bm{\varepsilon}$ with the response vector $\bm{y}\in\mathbb{R}^{n}$, the design matrix $\bm{X}\in\mathbb{R}^{n\times d}$ and the noise term satisfying $\mathbb{E}(\bm{\varepsilon})=\bm{0}$ and $\mathbb{V}(\bm{\varepsilon})=\sigma^{2}\bm{I}_{n}$. The unknown parameter $\bm{\beta}\in\mathbb{R}^{d}$ can be efficiently estimated by the method of least squares:

[TABLE]

which is of computational complexity $O(nd^{2}+d^{3})$. For the massive data with $n\gg d$, the least squares method would become computationally expensive, and it may exceed the computing capacity.

Faster least squares approximation can be achieved by the randomized sketch proposed by Drineas et al., (2011). It relies on a proper random matrix $\bm{S}\in\mathbb{R}^{m\times n}$ with $m\ll n$. The ways to generate random matrix $\bm{S}$ are divided into three main categories: subsampling, random projection and their hybrid. A widely known subsampling method is based on statistical leverage (LEV) scores (Drineas et al.,, 2006; Mahoney et al.,, 2011; Drineas et al.,, 2012; Ma et al.,, 2015). The random projection approach includes the subsampled randomized Hadamard transformation (SRHT) (Drineas et al.,, 2011; Boutsidis and Gittens,, 2013) and the Clarkson–Woodruff sketch (Clarkson and Woodruff,, 2013). The hybrid approach is to combine subsampling and random projection, see e.g. McWilliams et al., (2014).

When a sketching matrix $\bm{S}$ is fixed, there exist different sketching schemes, including the classical sketch (CS), the Hessian sketch (HS) and the iterative Hessian sketch (IHS). The widely adopted CS uses the sketched data pair $(\bm{SX},\bm{Sy})$ for approximating $\hat{\bm{\beta}}^{\rm LS}$ by

[TABLE]

Pilanci and Wainwright, (2016) showed that $\hat{\bm{\beta}}^{\rm CS}$ is suboptimal in the sense that it has a substantially larger error than $\hat{\bm{\beta}}^{\rm LS}$ with respect to the ground truth $\bm{\beta}^{*}$. They introduced the HS estimator based on the partially sketched data $(\bm{SX},\bm{y})$,

[TABLE]

and furthermore the IHS estimator based on iterative sketched data $\{\bm{S}_{t}\bm{X},t=1,\ldots,N\}$,

[TABLE]

where $\bm{e}_{t}=\bm{y}-\bm{X}\hat{\bm{\beta}}^{\rm IHS}_{t}$, with the initial $\hat{\bm{\beta}}^{\rm IHS}_{0}$ being provided. Unlike the CS and HS estimators, the IHS estimator is guaranteed to converge to $\bm{\hat{\beta}}^{\rm LS}$ upon some good event conditions (Pilanci and Wainwright,, 2016).

The IHS can be interpreted as a first-order gradient descent method with a series of random preconditioners $\bm{M}_{t}=\bm{X}^{T}\bm{S}_{t}^{T}\bm{S}_{t}\bm{X}$ subject to the unit step length. The preconditioner is widely used to boost optimization algorithms (Knyazev and Lashuk,, 2007; Gonen et al.,, 2016). However, the IHS is computationally intensive since in every iteration $\bm{M}_{t}^{-1}$ has to be evaluated for a new random sketching matrix $\bm{S}_{t}$. To speed up the IHS, Wang and Xu, (2018) proposed the pwGradient method to improve the IHS by a fixed well-designed sketching matrix, in which case the IHS reduces to a first order method with a constant preconditioner. Meanwhile, Wang et al., (2017) proposed another accelerated version of IHS with the fixed sketching matrix, while adopting conjugate gradient descent. Note that all these sketch methods are based on randomized sketching matrices.

In this paper we propose to reformulate the IHS estimator $\hat{\bm{\beta}}^{\rm IHS}_{t}$ by a linear combination of the initial $\hat{\bm{\beta}}^{\rm IHS}_{0}$ and the full data estimator $\hat{\bm{\beta}}^{\rm LS}_{0}$. Such reformulation enables us to find a sufficient isometric condition on the sketching matrices so that $\hat{\bm{\beta}}^{\rm IHS}_{t}$ is guaranteed to converge to $\bm{\hat{\beta}}^{\rm LS}$ with geometric convergence. It then motivates us to propose a deterministic approach for improving the IHS based on $A$-optimal subsampling and adaptive step lengths, which modifies the original second-order IHS method to be an adaptive first-order method. In summary, we improve the IHS method with the following three contributions:

• •

Good initialization. A good initialization scheme can reduce the inner iteration rounds of IHS while still delivering the same precision. We suggest to initialize the IHS with the classical sketch based on the $A$-optimal deterministic subsampling matrix. To our best knowledge, this is the first attempt to take initialization into account for improving the IHS method.

• •

Improved preconditioner. It is critical to find a well-designed preconditioner so that it may be fixed during the iterative sketch process. We propose to construct the preconditioner from $A$-optimal subsample and refine it by adding a ridge term. Unlike complicated random projection-based methods, we obtain our preconditioner at a low cost by recycling the subsamples in initialization and make no assumption on the sample size.

• •

Adaptive step lengths. We modify the IHS to be an adaptive first-order method by using a fixed preconditioner subject to variable step lengths. The step lengths at each iteration are determined by the exact line search, which ensures the algorithm to enjoy the guaranteed convergence.

Through extensive experiments, our proposed method is shown to achieve the state-of-art performance in terms of both precision and speed when approximating both the ground truth $\bm{\beta}^{*}$ and the full data estimator $\bm{\hat{\beta}}^{\rm LS}$.

2 Reformulation of IHS

The original IHS algorithm is displayed in Algorithm 1. For simplicity, we omit the superscript from $\hat{\bm{\beta}}^{\rm IHS}$. As for the sketch matrix, we use SRHT as it is widely adopted in the literature. Following Tropp, (2011) and Lu et al., (2013), when $n=2^{k}$ where $k$ is a positive integer, the SRHT sketch matrix is given by

[TABLE]

where $\bm{R}$ is an $m\times n$ matrix with rows chosen uniformly without replacement from the standard bases of $\mathbb{R}^{n}$, $\bm{H}$ is a normalized Walsh–Hadamard matrix of size $n\times n$, and $\bm{D}$ is an $n\times n$ diagonal matrix with i.i.d. Rademacher entries. Note that when $n$ is not the power of two, $\bm{X}$ is padded with zeros until $n$ achieves the next greater power of two.

In what follows, we present a novel reformulation of IHS. At each step, the update formula in Algorithm 1 can be rewritten as

[TABLE]

By (4), we can derive the following lemma by mathematical induction.

Lemma 1

Given an initializer $\hat{\bm{\beta}}_{0}$ and a series of independent sketch matrices $\{\bm{S}_{i}\}_{i=1}^{t}$, for any positive integer $t$, we have

[TABLE]

where $\bm{A}_{t}=\bm{M}_{t}^{-1}\bm{X}^{T}\bm{X}$ and $\bm{M}_{t}=\bm{X}^{T}\bm{S}^{T}_{t}\bm{S}_{t}\bm{X}$.

Lemma 1 reveals that $\hat{\bm{\beta}}_{t}$ is a linear combination of the initial $\hat{\bm{\beta}}^{\rm IHS}_{0}$ and the full data estimator $\hat{\bm{\beta}}^{\rm LS}_{0}$, with the re-weighting matrices during the iterations. We can therefore study the convergence properties of $\hat{\bm{\beta}}_{t}$ by looking into $\bm{A}_{t}$ or $\bm{M}_{t}$. The following theorem provides a sufficient isometric condition for the convergence guarantee.

Theorem 1

Given an initial estimator $\hat{\bm{\beta}}_{0}$ and a set of sketch matrices $\{\bm{S}_{t}\}_{t=1}^{\infty}$, the IHS estimator $\hat{\bm{\beta}}_{t}$ converges to $\hat{\bm{\beta}}^{\rm LS}$ with geometric convergence,

[TABLE]

provided that for any $\varepsilon_{1}\in(0,1/2)$ and $\varepsilon_{2}\in(0,1-\varepsilon_{1})$ and for any \bm{a}\in\mathbb{R}^{d}$$\backslash\{\bm{0}\},

[TABLE]

Proof 1

From Lemma 1, we have

[TABLE]

Thus, it holds that

[TABLE]

Note that

[TABLE]

where $\bm{Q}=\bm{X}^{T}\bm{X}$. Since $\bm{M}_{i}^{-1}\bm{Q}$ shares the same eigenvalues with $\bm{Q}^{1/2}\bm{M}_{i}^{-1}\bm{Q}^{1/2}$ (See 6.54 in Seber, (2008)), we have

[TABLE]

[TABLE]

From the conditions (7), we know that

[TABLE]

for any \bm{b}=\bm{Q}^{1/2}\bm{a}\in\mathbb{R}^{d}$$\backslash\{\bm{0}\}. Thus, we have that

[TABLE]

So, (8) becomes

[TABLE]

It is clear that the rate $\in(0,1)$, so $\hat{\bm{\beta}}_{t}$ converges to $\hat{\bm{\beta}}^{\rm LS}$. $\square$

Theorem 1 is meaningful in different ways. Firstly, when $\varepsilon_{1}=\rho,\varepsilon_{2}=\rho/2$ and $\rho\in(0,1/2)$, it can be checked that Theorem 1 corresponds to the main result of Pilanci and Wainwright, (2016) under the so-called good event condition. It also mimics the Johnson-Lindenstrauss lemma (Johnson and Lindenstrauss,, 1984) when $\varepsilon_{1}=\varepsilon_{2}$. Moreover, our theorem makes no assumption about the randomness of $\bm{S}_{t}$ and it is applicable to all kinds of sketch matrices satisfying (7). In the case of random $\bm{S}_{t}$, let $\varepsilon_{1}=\varepsilon_{2}=\varepsilon$, $\bm{U}=\bm{X}\bm{Q}^{-1/2}$ be the orthonormal basis of $\bm{X}$’s column space. It suffices to obtain (9) and hence (7) once $\bm{S}_{t}$ satisfies $\|\bm{I}_{d}-\bm{U}^{T}\bm{S}_{t}^{T}\bm{S}_{t}\bm{U}\|_{2}\leq\varepsilon$. Theorem 2.4 in Woodruff et al., (2014) showed one way to construct such random $\bm{S}_{t}$ with high probability.

The original IHS method requires calculating $\bm{M}_{t}$ with time complexity $O(nd\log(d))$ repeatedly for each $t$; see Section 2.5 of Pilanci and Wainwright, (2016). Theorem 1 indicates that a fixed sketch matrix such $\bm{S}_{t}=\bm{S}_{1},\forall t\geq 2$ can also ensure the convergence if the condition (7) is satisfied. This result provides enables us to reduce the computational cost for the original IHS method. For example, Wang and Xu, (2018) proposed the pwGradient method to improve the IHS by constructing a well-designed sketch algorithm, which can be viewed as an application of Theorem 1 by letting $\varepsilon_{1}=2\theta-\theta^{2}$, $\varepsilon_{2}=2\theta+\theta^{2}$ and $\theta\in(0,1/4)$.

In the meanwhile, we can also interpret the IHS approach based on a transformed space. Let the preconditioner $\bm{M}_{t}=\bm{M}$, multiply $\bm{M}^{1/2}$ to the both sides of (2), and denote $\bm{B}=\bm{M}^{-1/2}\bm{X}^{T}\bm{X}\bm{M}^{-1/2}$ and $\hat{\bm{\eta}}_{t}=\bm{M}^{1/2}\hat{\bm{\beta}}_{t}$ for $t=1,2,\ldots$. Then, we have

[TABLE]

which corresponds to the gradient descent update when minimizing the following least squares objective

[TABLE]

The one-to-one mapping between $\hat{\bm{\beta}}_{t}$ and $\hat{\bm{\eta}}_{t}$ ensures that one can show the convergence of $\hat{\bm{\beta}}_{t}$ and obtain its optima via equivalently analyzing $\hat{\bm{\eta}}_{t}$. Based on this observation, Wang et al., (2017) proposed the acc-IHS method by fixing the preconditioner and replacing the gradient descent with the conjugate counterpart in the transformed parameter space. It is worth mentioning that both the pwGradient and acc-IHS methods belong to the randomized approach based on random projections, while the SRHT sketching needs to operate on all the entries of $\bm{X}$. We hence seek to improve the IHS method in a more efficient and deterministic way.

3 Adaptive IHS with $A$-Optimal Subsampling

We extend the concept of sketch matrix from randomized settings to deterministic settings, by introducing $\bm{\delta}\in\mathbb{R}^{n}$ to indicate where an observation is selected, i.e., $\delta_{i}=1$ if sample $(\bm{x}_{i},y_{i})$ or $\bm{x}_{i}$ is included, $\delta_{i}=0$ otherwise. It is assumed that $\sum_{i=1}^{n}\delta_{i}=m$, which implies that the corresponding sketch matrix satisfies $\bm{S}^{T}\bm{S}=m^{-1}\mbox{diag}(\bm{\delta})$. It is our objective to find a good $\bm{\delta}$ subject to certain optimality criterion.

3.1 $A$-Optimal Classical Sketch

The convergence of $\hat{\bm{\beta}}_{t}$ also depends on $\|\hat{\bm{\beta}}_{0}-\hat{\bm{\beta}}^{\rm LS}\|_{2}$ as shown in (6), which motivates us to find a good initializer. We achieve this goal by proposing an $A$-optimal estimator under the classical sketch scheme. Suppose that we select a subset of $m$ observations. The least squares estimator based on the subdata and the corresponding covariance matrix are given by

[TABLE]

Let $\bm{M}(\bm{\delta})=\sum_{i=1}^{n}\delta_{i}\bm{x}_{i}\bm{x}_{i}^{T}$. Following the $A$-optimality criterion in experimental design (Pukelsheim,, 1993), we seek the subdata as indicated by $\bm{\delta}$ that minimizes the averaged variance of $\hat{\bm{\beta}}(\bm{\delta})$, which is proportional to the trace of $\bm{M}^{-1}(\bm{\delta})$. Formally, our goal can be formulated as the following discrete optimization problem,

[TABLE]

However, it is NP-hard to solve it exactly, so we turn to derive a upper bound of $\mbox{Tr}[\bm{M}^{-1}(\bm{\delta})]$, which leads to Algorithm 2.

Theorem 2

Let $\bm{Q}=\bm{X}^{T}\bm{X}$, $\kappa(\cdot)$ denote the condition number, $\mathcal{S}=\{\bm{\delta}|\bm{M}(\bm{\delta})>0,\sum_{i=1}^{n}\delta_{i}=m\}$ be the collection of all feasible $\bm{\delta}$ solutions and $\lambda_{\min}$ correspond to the smallest eigenvalue. Assuming that $C=\inf_{\bm{\delta}\in\mathcal{S}}\lambda_{\min}(\bm{M}(\bm{\delta}))>0$, for any given $\bm{\delta}\in\mathcal{S}$,

[TABLE]

Proof 2

Let $\bm{W}={\rm diag}(\bm{\delta})$, $\bar{\bm{W}}={\rm diag}(\bm{1}-\bm{\delta})$, $\lambda_{\max}$ be the largest eigenvalue, $\bm{X}=\bm{U}\bm{\Sigma}\bm{V}^{T}$ denote the SVD of $\bm{X}$ where $\bm{U}\in\mathbb{R}^{n\times d}$, $\bm{\Sigma}\in\mathbb{R}^{d\times d}$ and $\bm{V}\in\mathbb{R}^{d\times d}$. Therefore, we have

[TABLE]

It is easy to check that $\|\bm{U}^{T}\bar{\bm{W}}\bm{U}\|_{2}\leq\|\bm{U}^{T}\bm{U}\|_{2}\|\bar{\bm{W}}\|_{2}\leq 1$ where the first inequality follows $\bm{U}^{T}\bar{\bm{W}}\bm{U}\leq\bm{U}^{T}\bm{U}\|\bar{\bm{W}}\|_{2}$ and 10.47 in Seber, (2008). We can further conclude $\|\bm{U}^{T}\bar{\bm{W}}\bm{U}\|_{2}<1$ as $\bm{\delta}\in\mathcal{S}$. According to the extension of Corollary 5.6.16. in Horn and Johnson, (2012), we have

[TABLE]

So the trace in (13) can be further bounded as follows,

[TABLE]

Note that the denominator

[TABLE]

where the second and third equations follow the definition of eigenvalue and $\bm{X}$’s SVD respectively. With 6.76 in Seber, (2008), the denominator is lower bounded by

[TABLE]

Similarly, for the nominator, we have

[TABLE]

The inequality follows by plugging above two bounds into (14). $\square$

By Theorem 2, we can seek an approximately $A$-optimal design by shrinking its upper bound, i.e., maximizing $\sum_{i=1}^{n}\delta_{i}\|\bm{x}_{i}\|_{2}^{2}$. Our subsampling approach is described as follows.

Remark 1

The time for norm calculation is $O(nd)$, while sorting requires on average $O(n)$ operations (Martınez,, 2004). For step 3, it costs $O(md^{2}+d^{3})$ to obtain $\hat{\bm{\beta}}^{\rm CS}(\bm{\delta})$. In total, the time complexity of Algorithm 2 is $O(nd+md^{2})$. The time can be further reduced to $O(nd)$ when it comes to $n>md$, which is a common scenario for massive data. Furthermore, this constraint can also be a reference for the sketch size determination as it implies that $m$ should not exceed $n/d$ for an efficient complexity. From the time complexity perspective, our algorithm is as efficient as the $D$-optimality based method in Wang et al., (2018), while ours is more suitable for parallel computation.

Remark 2

In practice, when the data matrix $\bm{X}$ is centered and scaled initially, our algorithm tends to choose the extreme samples with the farthest distances to the center, which is consistent with the conclusion in Wang et al., (2018).

As shown in Figure 1, the performance of the $A$-optimal method uniformly dominates that of randomized sketch matrices and classical sketching scheme. Furthermore, the MSE of our approach decreases with an increase of $n$. In this regard, our $A$-optimal estimator serves as a good initialization for further enhancements.

3.2 Improved Preconditioner

Note that $\hat{\bm{\beta}}$ is mainly determined by a set of matrices $\{\bm{A}_{i}=\bm{M}_{i}^{-1}\bm{X}^{T}\bm{X}\}_{i=1}^{N}$. Rather than specifying $\bm{M}_{i}$’s by repeatedly sketching in IHS or fixing $\bm{M}_{i}=\bm{M}$ in pwGradient and acc-IHS, we make a compromise by defining $\bm{M}_{i}=\alpha_{i}^{-1}\bm{M}$, that is,

[TABLE]

where $\bm{A}=\bm{M}^{-1}\bm{X}^{T}\bm{X}$. We consider the $A$-optimal Hessian sketch in order to specify a deterministic preconditioner $\bm{M}$. Similarly, we can obtain the estimator of Hessian sketch and its covariance matrix based on an optimality criterion,

[TABLE]

where $\bm{M}(\bm{\delta})=n/m\sum_{i=1}^{n}\delta_{i}\bm{x}_{i}\bm{x}_{i}^{T}$. The next theorem provides a upper bound for the trace of the covariance matrix.

Theorem 3

Under the same conditions of Theorem 2, for any $\bm{\delta}\in\mathcal{S}$, we have

[TABLE]

Proof 3

Note that

[TABLE]

Therefore the above inequality can be easily obtained by applying Theorem 2. $\square$

Theorem 3 shows that the $A$-optimal design of Hessian sketched estimator can also be approximated via finding the samples with largest $\ell_{2}$ norm. The samples selected for initialization can be recycled for constructing the preconditioner. We find that adding a ridge term can further improve our preconditioner

[TABLE]

The rationale is demonstrated by Figure 2. Note that the effectiveness of the precondtioner $\bm{M}$ is measured by $\kappa(\bm{M}^{-1}\bm{X}^{T}\bm{X})=\kappa(\bm{M}^{-1/2}\bm{X}^{T}\bm{X}\bm{M}^{-1/2})=\kappa(\bm{B})$ (Benzi,, 2002), which should be close to 1. And $\kappa(\bm{B})$ can be visualized by the contour plot of $\tilde{f}(\bm{\eta})$ in the two-dimensional cases. Specifically, the flat degree of the contour plot indicates the size of the condition number. The larger the condition number, the more circular the contour plot. The optimal transformation is $\bm{M}=\bm{X}^{T}\bm{X}$ as shown by Figure 2(b). Observing Figures 2(a) and 2(c), $\kappa(\bm{B})$ is still close to the $\kappa(\bm{X}^{T}\bm{X})$, i.e., the preconditioner fails. After ridging, our preconditioner shrinks the major radius in Figure 2(c) and render the transformation to reach optimality as shown by Figure 2(d).

3.3 Exact Line Search

Using the improved preconditioner of the form (16), the adaptive first-order IHS estimator (15) can be rewritten as

[TABLE]

for $t=1,2,\ldots,N$. To determine $\alpha_{t}$ adaptively, we note that it actually corresponds to the learning rate of the first order method in the transformed space. Recall that the learning rate in (10) is fixed as unit 1, which may be suboptimal and does not lead to the largest descent for every iteration. We hence fulfill the potential of IHS by exact line search. Specifically, let

[TABLE]

be the update direction at the $t$th iteration, the optimal step lengths are determined by adaptively minimizing the univariate function $\psi(\alpha_{t})=\tilde{f}(\hat{\bm{\eta}}_{t-1}+\alpha_{t}\bm{d}_{t})$. It is easy to show that

[TABLE]

where $f(\bm{\beta})=1/2\|\bm{X}\bm{\beta}-\bm{y}\|_{2}^{2}$. Moreover, with the above $\{\alpha_{t}\}_{t=1}^{N}$, the first-order method becomes steepest descent in the transformed space with the guaranteed convergence; see Theorem 3.3 in Nocedal and Wright, (2006). Thus it also converges in the original space by the one-to-one mapping between $\hat{\bm{\eta}}_{t}$ and $\hat{\bm{\beta}}_{t}$.

In summary, we present the adaptive $A$-optimal IHS in Algorithm 3. The algorithm uses a pre-specified loop number $N$, and it can be modified by early stopping strategy, e.g. when the difference between $\hat{\bm{\beta}}_{t}$ and $\hat{\bm{\beta}}_{t-1}$ is below a tolerance threshold.

4 Numerical Experiments

4.1 Simulated Data Analysis

Data Generation.

We generate data by following the experimental setups of Wang et al., (2018). All data are generated from the linear model with the true parameter $\bm{\beta}^{*}$ being $d$ i.i.d $N(0,1)$ variates and $\sigma^{2}=9$. Let $\bm{\Sigma}$ be the covariance matrix where its $(i,j)$ entry is $\Sigma_{ij}=0.5^{I(i\neq j)}$ for $i,j\in[d]$, and $I(\cdot)$ is the indicator function. We consider four multivariate distributions for covariates $\{\bm{x}_{i}\}_{i=1}^{n}$.

1. (1)

A multivariate normal distribution $N(\bm{0},\bm{\Sigma})$. 2. (2)

A multivariate log-normal distribution $\mbox{\rm LN}(\bm{0},\bm{\Sigma})$ which is generated by taking the exponential transformation of $N(\bm{0},\bm{\Sigma})$. 3. (3)

A multivariate $t$ distribution with $2$ degrees of freedom $t_{2}(\bm{0},\bm{\Sigma})$. 4. (4)

A mixture distribution composed of five different distributions $N(\bm{1},\bm{\Sigma})$, $t_{2}(\bm{0},\bm{\Sigma})$, $t_{3}(\bm{0},\bm{\Sigma})$, $\mbox{\rm Unif}(\bm{0},\bm{2})$, \mbox{\rm LN}({\color[rgb]{0,0,0}{}\bm{0}},\bm{\Sigma}) with equal proportions, where $\mbox{\rm Unif}(\bm{0},\bm{2})$ represents $d$ elements from independent uniform distributions.

To remove the effect of the intercept, we center the data. Unless particularly stated, we set $n=2^{17},m=1000$ and present the results averaged over $R=1000$ replications.

Choice of Preconditioner.

To assess the quality of the preconditioner $\bm{M}$ that is a positive definite matrix, we first define the $\Delta(\cdot)$ measure as the ratio of improvement with respect to the condition number,

[TABLE]

For a good preconditioner, its $\Delta$ measure should be positive and close to 1. In contrast, a preconditioner with negative $\Delta$ value would worsen the condition number $\kappa(\cdot)$.

To choose the ridge parameter $\lambda$, we calculate the averaged $\Delta$ values on different distributions versus $\lambda$ as a proportion of $\sum_{i=1}^{n}\|\bm{x}_{i}\|^{2}_{2}$. Specifically, let $d\in\{50,100\}$ and the proportion take values from $0.1$ to $1$ with step size $0.1$. The results are plotted in Figure 3, which shows that the curves are similar on heavy-tailed distributions but not so consistent for the normal case. Through extensive experiments, as a rule of thumb, we suggest that $\lambda$ takes value $0.1\sum_{i=1}^{n}\|\bm{x}_{i}\|^{2}_{2}$ for relatively concentrated data and $0.4\sum_{i=1}^{n}\|\bm{x}_{i}\|^{2}_{2}$ for heavy-tailed distributions.

Table 1 presents the averaged $\Delta$ values of the proposed preconditioner (16) for different distributions with dimensionality $d=50$ and $d=100$, respectively. We use $\tilde{\lambda}$ to denote $0.1\sum_{i=1}^{n}\|\bm{x}_{i}\|^{2}_{2}$ for the normal distribution and $0.4\sum_{i=1}^{n}\|\bm{x}_{i}\|^{2}_{2}$ for other cases. The $\Delta$ results are compared with $\lambda=0$ and the SRHT scheme. It is clear that the preconditioner with parameter $\tilde{\lambda}$ uniformly dominates in all cases. Such superiority is especially significant for the cases with $d=100$, i.e., when $d/m$ is relatively large. Furthermore, note that when $\bm{M}=\lambda\bm{I}_{d}$, $\kappa(\bm{M}^{-1/2}\bm{X}^{T}\bm{X}\bm{M}^{-1/2})=\kappa(\bm{X}^{T}\bm{X}/\lambda)=\kappa(\bm{X}^{T}\bm{X})$ and $\Delta=0$. Therefore the ridged combination indeed brings significant improvement on the preconditioner.

Comparative Study.

We compare the proposed deterministic $A$-optimal IHS to the original randomized IHS Pilanci and Wainwright, (2016), as well as two other improved IHS methods: acc-IHS Wang et al., (2017) and pwGradient Wang and Xu, (2018). The initial estimators of these benchmark methods are set to be the vector zero as suggested in the original papers. Two empirical MSE criteria are used to evaluate the algorithm:

[TABLE]

where $\hat{\bm{\beta}}_{ti}$ denotes the estimator of the $t$th iteration in the $i$th round, and the trimmed mean with fraction 0.025 is applied to handle the effects of abnormal situations. We record MSEs from the [math]th iteration (i.e. initialization stage) and set the sketch dimension to be $m$ at every iteration. Note that the randomized IHS requires independent sketch matrix for each iteration, while the other methods reuse the same sketch matrix for each iteration.

Figures 4 and 5 illustrate the advantages of our proposed Aopt-IHS method. Firstly, the superiority of initialization can be thoroughly reflected by the flat curve in Figure 4, as well as the lower starting values in Figure 5. During the iterations, the outperformance of the proposed method can be reflected in Figure 5. Specifically, given the convergent property of $\hat{\bm{\beta}}_{t}$, the curve slope reflects the convergence rate while the intercept or the starting point measures the effect of the initialization. From the plots, the proposed Aopt-IHS method not only outperforms in the initialization but also demonstrates a decent convergence rate that is competitive to pwGradient and it is even sharper than IHS for the normal case.

To compare the computational time, we present in Table 2 the averaged run time for each method for achieving the precision $10^{-10}$ in $\|\hat{\bm{\beta}}-\hat{\bm{\beta}}^{\rm LS}\|_{2}$, as well as the required numbers of iterations. Note that the time for computing $\bm{X}^{T}\bm{y}$ in the update formula is not counted for all methods. Table 2 shows that our method can attain the same precision with fewer iterations and shorter time.

4.2 Real Data Analysis

Finally, we apply the proposed method to a real food intake dataset111Data can be found in https://www.icpsr.umich.edu/icpsrweb/ICPSR/studies/21960#.. This data, consisting of three sub-files, records the information on one-day dietary intakes of male residents in the United States with age from 19 to 50 in 1985. We consider the largest food intake file that includes how much and when the food items are taken. This sub-file has $n=10177$ samples and we are interested in the linear regression model between the food energy ($y$) and the four features: total protein ($x_{1}$), fat ($x_{2}$), carbohydrate ($x_{3}$), and alcohol ($x_{4}$). Both the feature and target variables are centered from we feed them to the linear model.

Using the IHS methods with sketch size as $m=500$ and $\lambda=0.1\sum_{i=1}^{n}\|\bm{x}_{i}\|^{2}_{2}$, we evaluate their performances through the Euclidean distance between the estimated parameter and the full data least square estimator $\|\hat{\bm{\beta}}-\hat{\bm{\beta}}^{\rm LS}\|_{2}$. Same as the simulation study above, we take the trimmed average of Euclidean errors over $R=1000$ repetitions for the randomized methods. The comparative results are plotted in Figure 6 and and it is shown that our proposed Aopt-IHS has the best performance during the beginning iterations. It outperforms both acc-IHS and pwGradient for iterations up to 20. In this real data analysis, the randomized IHS has the best convergence rate, but computationally it is the least efficient method.

5 Conclusion

We reformulate the IHS method as an adaptive first-order optimization method, by using the idea of optimal design of experiments for subdata selection. To our best knowledge, this is the first attempt to improve the IHS method in a deterministic manner while maintaining a decent speed and precision. The numerical experiments confirm the superiority of the proposed approach.

There are several open problems worth of further investigation, including the theoretical properties of the ridged preconditioner according to the conditions derived in Theorem 1. It is of our future interest to investigate how the ridge term affects the convergence rate of Algorithm 3. Other than using the $A$-optimal design, it is also interesting to investigate the $D$-optimal or other types of optimal designs for the purpose of subdata selection, in particular when the data are heterogeneously distributed.

Appendix A: Extra Results with Different Ridged Preconditioners

We further compare the ridged preconditioner $\bm{M}=\frac{n}{m}\sum_{i=1}^{n}\delta_{i}\bm{x}_{i}\bm{x}_{i}^{T}+\lambda\bm{I}_{d}$ with its two components, the non-ridged term and the scaled identity matrix. Three preconditioners are evaluated through ${\rm MSE}_{2}$ under our proposed algorithm framework. We only consider the identity matrix $\bm{M}=\bm{I}$ since any scaling multiplier $\lambda$ in $\bm{M}=\lambda\bm{I}$ can be canceled out during the update of $\hat{\bm{\beta}}_{t}$. The results strengthen that the ridging operation may enhance the preconditioner performance.

Appendix B: Extra Results with the Same Proposed Initial Estimator

In this section, we perform some additional experiments where all the methods are initialized by our proposed $A$-optimal estimator. The subsample size is fixed as $m=1000$. These experiments further justify that the proposed Aopt-IHS method generally enjoys the better convergent rates than the benchmark methods.