Iterative Hessian sketch: Fast and accurate solution approximation for constrained least-squares

Introduction

There are various ways in which the quality of the approximate solution $\widetilde{x}$ can be assessed. One standard way is in terms of the minimizing value of the quadratic cost function $f$ defining the original problem (1), which we refer to as cost approximation. In terms of $f$-cost, the approximate solution $\widetilde{x}$ is said to be $\varepsilon$-optimal if

We refer to these measures as solution approximation.

Now of course, a cost approximation bound (3) can be used to derive guarantees on the solution approximation error. However, it is natural to wonder whether or not, for a reasonable sketch size, the resulting guarantees are “good”. For instance, using arguments from Drineas et al. , for the problem of unconstrained least-squares, it can be shown that the same conditions ensuring a $\varepsilon$-accurate cost approximation also ensure that

is required. This scaling is undesirable in the regime $n\gg d$, where the whole point of sketching is to have the sketch dimension $m$ much lower than $n$.

Now the alert reader will have observed that the preceding argument was only rough and heuristic. However, the first result of this paper (Theorem 1) provides a rigorous confirmation of the conclusion: whenever $m\ll n$, the classical least-squares sketch (2) is sub-optimal as a method for solution approximation. Figure 1 provides an empirical demonstration of the poor behavior of the classical least-squares sketch for an unconstrained problem.

This sub-optimality holds not only for unconstrained least-squares but also more generally for a broad class of constrained problems. Actually, Theorem 1 is a more general claim: any estimator based only on the pair $(SA,Sy)$—an infinite family of methods including the standard sketching algorithm as a particular case—is sub-optimal relative to the original least-squares estimator in the regime $m\ll n$. We are thus led to a natural question: can this sub-optimality be avoided by a different type of sketch that is nonetheless computationally efficient? Motivated by this question, our second main result (Theorem 2) is to propose an alternative method—known as the iterative Hessian sketch—and prove that it yields optimal approximations to the least-squares solution using a projection size that scales with the intrinsic dimension of the underlying problem, along with a logarithmic number of iterations. The main idea underlying iterative Hessian sketch is to obtain multiple sketches of the data $(S^{1}A,...,S^{N}A)$ and iteratively refine the solution where $N$ can be chosen logarithmic in $n$.

The remainder of this paper is organized as follows. In Section 2, we begin by introducing some background on classes of random sketching matrices, before turning to the statement of our lower bound (Theorem 1) on the classical least-squares sketch (2). We then introduce the Hessian sketch, and show that an iterative version of it can be used to compute $\varepsilon$-accurate solution approximations using $\log(1/\varepsilon)$-steps (Theorem 2). In Section 3, we illustrate the consequences of this general theorem for various specific classes of least-squares problems, and we conclude with a discussion in Section 4. The majority of our proofs are deferred to the appendices.

For the convenience of the reader, we summarize some standard notation used in this paper. For sequences $\{a_{t}\}_{t=0}^{\infty}$ and $\{b_{t}\}_{t=0}^{\infty}$, we use the notation $a_{t}\preceq b_{t}$ to mean that there is a constant (independent of $t$) such that $a_{t}\leq C\,b_{t}$ for all $t$. Equivalently, we write $b_{t}\succeq a_{t}$. We write $a_{t}\asymp b_{t}$ if $a_{t}\preceq b_{t}$ and $b_{t}\preceq a_{t}$.

Main results

In this section, we begin with background on different classes of randomized sketches, including those based on random matrices with sub-Gaussian entries, as well as those based on randomized orthonormal systems and random sampling. In Section 2.2, we prove a general lower bound on the solution approximation accuracy of any method that attempts to approximate the least-squares problem based on observing only the pair $(SA,Sy)$. This negative result motivates the investigation of alternative sketching methods, and we begin this investigation by introducing the Hessian sketch in Section 2.3. It serves as the basic building block of the iterative Hessian sketch (IHS), which can be used to construct an iterative method that is optimal up to logarithmic factors.

The second type of randomized sketch we consider is randomized orthonormal system (ROS), for which matrix multiplication can be performed much more efficiently.

Given a probability distribution $\{p_{j}\}_{j=1}^{n}$ over $[n]=\{1,\ldots,n\}$, another choice of sketch is to randomly sample the rows of the extended data matrix $\begin{bmatrix}A&y\end{bmatrix}$ a total of $m$ times with replacement from the given probability distribution. Thus, the rows of $S$ are independent and take on the values

for some constant $\alpha$ independent of $n$.

In the following section, we present a lower bound that applies to all the three kinds of sketching matrices described above.

2 Sub-optimality of classical least-squares sketch

We begin by proving a lower bound on any estimator that is a function of the pair $(SA,Sy)$. In order to do so, we consider an ensemble of least-squares problems, namely those generated by a noisy observation model of the form

With this set-up, we have the following result:

where $M_{1/2}$ is the $1/2$-packing number of $\mathcal{C}_{0}$ in the semi-norm $\|\cdot\|_{A}$.

The proof, given in Appendix A, is based on a reduction from statistical minimax theory combined with information-theoretic bounds. The lower bound is best understood by considering some concrete examples:

Consequently, the sketch dimension $m$ must grow proportionally to $n$ in order for the sketched solution to have a mean-squared error comparable to the original least-squares estimate. This is highly undesirable for least-squares problems in which $n\gg d$, since it should be possible to sketch down to a dimension proportional to $\operatorname{rank}(A)=d$. Thus, Theorem 1 this reveals a surprising gap between the classical least-squares sketch (2) and the accuracy of the original least-squares estimate.

In contrast, the sketching method of this paper, known as iterative Hessian sketching (IHS), matches the optimal mean-squared error using a sketch of size $d+\log(n)$ in each round, and a total of $\log(n)$ rounds; see Corollary 2 for a precise statement. The red curves in Figure 1 show that the mean-squared errors ($\|\widehat{x}-x^{*}\|_{2}^{2}$ in panel (a), and $\|\widehat{x}-x^{*}\|_{A}^{2}$ in panel (b)) of the IHS method using this sketch dimension closely track the associated errors of the full least-squares solution (blue curves). Consistent with our previous discussion, both curves drop off at the $n^{-1}$ rate.

Since the IHS method with $\log(n)$ rounds uses a total of T=\log(n)\big{\{}d+\log(n)\} sketches, a fair comparison is to implement the classical method with $T$ sketches in total. The black curves show the MSE of the resulting sketch: as predicted by our theory, these curves are relatively flat as a function of sample size $n$. Indeed, in this particular case, the lower bound (10)

showing we can expect (at best) an inverse logarithmic drop-off. $\diamondsuit$

This sub-optimality can be extended to other forms of constrained least-squares estimates as well, such as those involving sparsity constraints.

On the other hand, the $\frac{1}{2}$-packing number $M$ of the set $\mathcal{C}_{0}$ can be lower bounded as \log M\succsim s\log\big{(}\frac{ed}{s}\big{)}; see Appendix D.2 for the details of this calculation. Consequently, in application to this particular problem, Theorem 1 implies that any estimator $x^{\dagger}$ based on the pair $(SA,Sy)$ has mean-squared error lower bounded as

Again, we see that the projection dimension $m$ must be of the order of $n$ in order to match the mean-squared error of the constrained least-squares estimate $x^{\mbox{\tiny{LS}}}$ up to constant factors. By contrast, in this special case, the sketching method developed in this paper matches the error $\|x^{\mbox{\tiny{LS}}}-x^{*}\|_{2}$ using a sketch dimension that scales only as s\log\big{(}\frac{ed}{s}\big{)}+\log(n); see Corollary 3 for the details of a more general result. $\diamondsuit$

where $\gamma_{j}(X)$ denotes the $j^{th}$ singular value of $X$. This observation motivates a standard relaxation of the rank constraint using the nuclear norm $|\!|\!|X|\!|\!|_{{\mbox{\tiny{nuc}}}}:\,=\sum_{j=1}^{\min\{d_{1},d_{2}\}}\gamma_{j}(X)$.

Accordingly, let us consider the constrained least-squares problem

where $|\!|\!|\cdot|\!|\!|_{{\mbox{\tiny{fro}}}}$ denotes the Frobenius norm on matrices, or equivalently the Euclidean norm on its vectorized version. Let $\mathcal{C}_{0}$ denote the set of matrices with rank $r<\frac{1}{2}\min\{d_{1},d_{2}\}$, and Frobenius norm at most one. In this case, we show in Appendix D that the constrained least-squares solution $X^{\mbox{\tiny{LS}}}$ satisfies the bound

As with the previous examples, we see the sub-optimality of the sketched approach in the regime $m<n$. In contrast, for this class of problems, our sketching method matches the error $\|X^{\mbox{\tiny{LS}}}-X^{*}\|_{A}$ using a sketch dimension that scales only as $\{r(d_{1}+d_{2})+\log(n)\}\,\log(n$). See Corollary 4 for further details.

3 Introducing the Hessian sketch

In controlling the error with respect to the least-squares solution $x^{\mbox{\tiny{LS}}}$ the set of possible descent directions $\{x-x^{\mbox{\tiny{LS}}}\,\mid\,x\in\mathcal{C}\}$ plays an important role. In particular, we define the transformed tangent cone

Note that the error vector $\widehat{v}:\,=A(\widehat{x}-x^{\mbox{\tiny{LS}}})$ of interest belongs to this cone. Our approximation bound is a function of the quantities

where $u$ is a fixed unit-norm vector. These variables played an important role in our previous analysis of the classical sketch (2). The following bound applies in a deterministic fashion to any sketching matrix.

For random sketching matrices, Proposition 1 can be combined with probabilistic analysis to obtain high probability error bounds. For a given tolerance parameter $\rho\in(0,\frac{1}{2}]$, consider the “good event”

where the final inequality holds for all $\rho\in(0,1/2]$.

Thus, for a given family of random sketch matrices, we need to choose the projection dimension $m$ so as to ensure the event $\mathcal{E}{\rho}$ holds for some $\rho$. For future reference, let us state some known results for the cases of sub-Gaussian and ROS sketching matrices. We use $(c_{0},c_{1},c_{2})$ to refer to numerical constants, and we let $D=\dim(\mathcal{C})$ denote the dimension of the space $\mathcal{C}$. In particular, we have $D=d$ for vector-valued estimation, and $D=d_{1}d_{2}$ for matrix problems.

Our bounds involve the “size” of the cone $\mathcal{K}^{\mbox{\tiny{LS}}}$ previously defined (17), as measured in terms of its Gaussian width

where $g\sim N(0,I_{n})$ is a standard Gaussian vector. With this notation, we have the following:

For sub-Gaussian sketch matrices, given a sketch size $m>\frac{c_{0}}{\rho^{2}}\mathcal{W}^{2}(\mathcal{K}^{\mbox{\tiny{LS}}})$, we have

For randomized orthogonal system (ROS) sketches over the class of self-bounding cones, given a sketch size $m>\frac{c_{0}\,\log^{4}(D)}{\rho^{2}}\mathcal{W}^{2}(\mathcal{K}^{\mbox{\tiny{LS}}})$, we have

4 Iterative Hessian sketch

Despite the deficiency of the Hessian sketch itself, it serves as the building block for an novel scheme—known as the iterative Hessian sketch—that can be used to match the accuracy of the least-squares solution using a reasonable sketch dimension. Let begin by describing the underlying intuition. As summarized by the bound (20b), conditioned on the good event $\mathcal{E}(\rho)$, the Hessian sketch returns an estimate with error within a $\rho$-factor of $\|x^{\mbox{\tiny{LS}}}\|_{A}$, where $x^{\mbox{\tiny{LS}}}$ is the solution to the original unsketched problem. As show by Lemma 1, as long as the projection dimension $m$ is sufficiently large, we can ensure that $\mathcal{E}(\rho)$ holds for some $\rho\in(0,1/2)$ with high probability. Accordingly, given the current iterate $x^{t}$, suppose that we can construct a new least-squares problem for which the optimal solution is $x^{\mbox{\tiny{LS}}}-x^{t}$. Applying the Hessian sketch to this problem will then produce a new iterate $x^{t+1}$ whose distance to $x^{\mbox{\tiny{LS}}}$ has been reduced by a factor of $\rho$. Repeating this procedure ${N}$ times will reduce the initial approximation error by a factor $\rho^{N}$.

With this intuition in place, we now turn a precise formulation of the iterative Hessian sketch. Consider the optimization problem

where $x^{t}$ is the iterate at step $t$. By construction, the optimum to this problem is given by $\widehat{u}=x^{\mbox{\tiny{LS}}}-x^{t}$. We then apply to Hessian sketch to this optimization problem (23) in order to obtain an approximation $x^{t+1}=x^{t}+\widehat{u}$ to the original least-squares solution $x^{\mbox{\tiny{LS}}}$ that is more accurate than $x^{t}$ by a factor $\rho\in(0,1/2)$. Recursing this procedure yields a sequence of iterates whose error decays geometrically in $\rho$.

Formally, the iterative Hessian sketch algorithm takes the following form:

The following theorem summarizes the key properties of this algorithm. It involves the sequence $\{Z_{1}(S^{t}),Z_{2}(S^{t})\}_{t=1}^{N}$, where the quantities $Z_{1}$ and $Z_{2}$ were previously defined in equations (18a) and (18b). In addition, as a generalization of the event (20a), we define the sequence of “good” events

With this notation, we have the following guarantee:

The final solution $\widehat{x}=x^{{N}}$ satisfies the bound

Note that for any $\rho\in(0,1/2)$, then event $\mathcal{E}^{t}(\rho)$ implies that $\frac{Z_{2}(S^{t})}{Z_{1}(S^{t})}\leq\rho$, so that the bound (26b) is an immediate consequence of the product bound (26a).

Lemma 1 can be combined with the union bound in order to ensure that the compound event $\cap_{t=1}^{N}\mathcal{E}^{t}(\rho)$ holds with high probability over a sequence of $N$ iterates, as long as the sketch size is lower bounded as $m\geq\frac{c_{0}}{\rho^{2}}\mathcal{W}^{2}(\mathcal{K}^{\mbox{\tiny{LS}}})\log^{4}(D)+\log{N}$. Based on the bound (26b), we then expect to observe geometric convergence of the iterates.

In order to test this prediction, we implemented the IHS algorithm using Gaussian sketch matrices, and applied it to an unconstrained least-squares problem based on a data matrix with dimensions $(d,n)=(200,6000)$ and noise variance $\sigma^{2}=1$. As shown in Appendix D.2, the Gaussian width of $\mathcal{K}^{\mbox{\tiny{LS}}}$ is proportional to $d$, so that Lemma 1 shows that it suffices to choose a projection dimension $m\succsim\gamma d$ for a sufficiently large constant $\gamma$. Panel (a) of Figure 2 illustrates the resulting convergence rate of the IHS algorithm, measured in terms of the error $\|x^{t}-x^{\mbox{\tiny{LS}}}\|_{A}$, for different values $\gamma\in\{4,6,8\}$. As predicted by Theorem 2, the convergence rate is geometric (linear on the log scale shown), with the rate increasing as the parameter $\gamma$ is increased.

Assuming that the sketch dimension has been chosen to ensure geometric convergence, Theorem 2 allows us to specify, for a given target accuracy $\varepsilon\in(0,1)$, the number of iterations required.

Fix some $\rho\in(0,1/2)$, and choose a sketch dimension $m>\frac{c_{0}\log^{4}(D)}{\rho^{2}}\mathcal{W}^{2}(\mathcal{K}^{\mbox{\tiny{LS}}})$. If we apply the IHS algorithm for ${N}(\rho,\varepsilon):\,=1+\frac{\log(1/\varepsilon)}{\log(1/\rho)}$ steps, then the output $\widehat{x}=x^{{N}}$ satisfies the bound

with probability at least $1-c_{1}{N}(\rho,\varepsilon)e^{-c_{2}\frac{m\rho^{2}}{\log^{4}(D)}}$.

This corollary is an immediate consequence of Theorem 2 combined with Lemma 1, and it holds for both ROS and sub-Gaussian sketches. (In the latter case, the additional $\log(D)$ terms may be omitted.) Combined with bounds on the width function $\mathcal{W}(\mathcal{K}^{\mbox{\tiny{LS}}})$, it leads to a number of concrete consequences for different statistical models, as we illustrate in the following section.

One way to understand the improvement of the IHS algorithm over the classical sketch is as follows. Fix some error tolerance $\varepsilon\in(0,1)$. Disregarding logarithmic factors, our previous results on the classical sketch then imply that a sketch size $m\succsim\varepsilon^{-2}\>\mathcal{W}^{2}(\mathcal{K}^{\mbox{\tiny{LS}}})$ is sufficient to produce a $\varepsilon$-accurate solution approximation. In contrast, Corollary 1 guarantees that a sketch size $m\succsim\log(1/\varepsilon)\;\mathcal{W}^{2}(\mathcal{K}^{\mbox{\tiny{LS}}})$ is sufficient. Thus, the benefit is the reduction from $\varepsilon^{-2}$ to $\log(1/\varepsilon)$ scaling of the required sketch size.

It is worth noting that in the absence of constraints, the least-squares problem reduces to solving a linear system, so that alternative approaches are available. For instance, one can use a randomized sketch to obtain a preconditioner, which can then be used within the conjugate gradient method. As shown in past work , two-step methods of this type can lead to same reduction of $\varepsilon^{-2}$ dependence to $\log(1/\varepsilon)$. However, a method of this type is very specific to unconstrained least-squares, whereas the procedure described in this paper is generally applicable to least-squares over any compact, convex constraint set.

5 Computational and space complexity

Let us now make a few comments about the computational and space complexity of implementing the IHS algorithm using ROS sketches (e.g., such as those based on the fast Hadamard transform). For a given sketch size $m$, at iteration $t$, the IHS algorithm requires $\mathcal{O}(nd\log(m))$ basic operation for computing the data sketch $S^{t+1}A$ and also $\mathcal{O}(nd)$ operations to compute $A^{T}(y-Ax^{t})$. Consequently, if we run the algorithm for ${N}$ iterations, then the overall complexity is \mathcal{O}\big{(}(nd\log(m)+C(m,d))\,{N}\big{)}, where $C(m,d)$ is the complexity of solving the $m\times d$ dimensional problem in the update (24). The total space used scales as $\mathcal{O}(md)$.

If we want to obtain estimates with accuracy $\varepsilon$, then we need to perform $N\asymp\log(1/\varepsilon)$ iterations in total. Moreover, for ROS sketches, we need to choose $m\succsim\mathcal{W}^{2}(\mathcal{K}^{\mbox{\tiny{LS}}})\log^{4}(d)$. Consequently, it only remains to bound the Gaussian width $\mathcal{W}$ in order to specify complexities that depend only on the pair $(n,d)$.

For an unconstrained problem with $n>d$, the Gaussian width can be bounded as $\mathcal{W}^{2}(\mathcal{K}^{\mbox{\tiny{LS}}})\precsim d$, and the complexity of the solving the sub-problem (24) can be bounded as $d^{3}$. Thus, the overall complexity of computing an $\varepsilon$-accurate solution scales as $\mathcal{O}(nd\log(d)+d^{3})\log(1/\varepsilon)$, and the space required is $\mathcal{O}(d^{2})$.

Consequences for concrete models

In this section, we derive some consequences of Corollary 1 for particular classes of least-squares problems. Our goal is to provide empirical confirmation of the sharpness of our theoretical predictions, namely the minimal sketch dimension required in order to match the accuracy of the original least-squares solution.

second, choose a regression vector $x^{*}$ uniformly at random from the sphere $\mathcal{S}^{d-1}$

third, form the response vector $y=Ax^{*}+w$, where $w\sim N(0,\sigma^{2}I_{n})$ is observation noise with $\sigma=1$.

As discussed following Lemma 1, for this class of problems, taking a sketch dimension $m\succsim\frac{d}{\rho^{2}}$ guarantees $\rho$-contractivity of the IHS iterates with high probability. Consequently, we can obtain a $\varepsilon$-accurate approximation to the original least-squares solution by running roughly $\log(1/\varepsilon)/\log(1/\rho)$ iterations.

Now how should the tolerance $\varepsilon$ be chosen? Recall that the underlying reason for solving the least-squares problem is to approximate $x^{*}$. Given this goal, it is natural to measure the approximation quality in terms of $\|x^{t}-x^{*}\|_{A}$. Panel (b) of Figure 2 shows the convergence of the iterates to $x^{*}$. As would be expected, this measure of error levels off at the ordinary least-squares error

Consequently, it is reasonable to set the tolerance parameter proportional to $\sigma^{2}\frac{d}{n}$, and then perform roughly $1+\frac{\log(1/\varepsilon)}{\log(1/\rho)}$ steps. The following corollary summarizes the properties of the resulting procedure:

For some given $\rho\in(0,1/2)$, suppose that we run the IHS algorithm for ${N}=1+\lceil\frac{\log\sqrt{n}\;\frac{\|x^{\mbox{\tiny{LS}}}\|_{A}}{\sigma}}{\log(1/\rho)}\rceil$ iterations using $m=\frac{c_{0}}{\rho^{2}}d$ projections per round. Then the output $\widehat{x}$ satisfies the bounds

with probability greater than $1-c_{1}\,{N}\,e^{-c_{2}\frac{m\rho^{2}}{\log^{4}(d)}}$.

In order to confirm the predicted bound (28) on the error $\|\widehat{x}-x^{\mbox{\tiny{LS}}}\|_{A}$, we performed a second experiment. Fixing $n=100d$, we generated $T=20$ random least squares problems from the ensemble described above with dimension $d$ ranging over $\{32,64,128,256,512\}$. By our previous choices, the least-squares estimate should have error $\|x^{\mbox{\tiny{LS}}}-x^{*}\|_{2}\approx\sqrt{\frac{\sigma^{2}d}{n}}=0.1$ with high probability, independently of the dimension $d$. This predicted behavior is confirmed by the blue bars in Figure 3; the bar height corresponds to the average over $T=20$ trials, with the standard errors also marked. On these same problem instances, we also ran the IHS algorithm using $m=6d$ samples per iteration, and for a total of

Since $\|x^{\mbox{\tiny{LS}}}-x^{*}\|_{A}\asymp\sqrt{\frac{\sigma^{2}d}{n}}\approx 0.10$, Corollary 2 implies that with high probability, the sketched solution $\widehat{x}=x^{{N}}$ satisfies the error bound

for some constant $c^{\prime}_{0}>0$. This prediction is confirmed by the green bars in Figure 3, showing that $\|\widehat{x}-x^{*}\|_{A}\approx 0.11$ across all dimensions. Finally, the red bars show the results of running the classical sketch with a sketch dimension of $(6\times 4)d=24d$ sketches, corresponding to the total number of sketches used by the IHS algorithm. Note that the error is roughly twice as large.

2 Sparse least-squares

Under these conditions, the proof of Corollary 3 shows that a sketch size m\geq\gamma\,s\log\big{(}\frac{ed}{s}\big{)} suffices to guarantee geometric convergence of the IHS updates. Panel (a) of Figure 4 illustrates the accuracy of this prediction, showing the resulting convergence rate of the the IHS algorithm, measured in terms of the error $\|x^{t}-x^{\mbox{\tiny{LS}}}\|_{A}$, for different values $\gamma\in\{2,5,25\}$. As predicted by Theorem 2, the convergence rate is geometric (linear on the log scale shown), with the rate increasing as the parameter $\gamma$ is increased.

As long as n\succsim s\log\big{(}\frac{ed}{s}\big{)}, it also follows as a corollary of Proposition 2 that

with high probability. This bound suggests an appropriate choice for the tolerance parameter $\varepsilon$ in Theorem 2, and leads us to the following guarantee.

For the stated random ensemble of sparse linear regression problems, suppose that we run the IHS algorithm for ${N}=1+\lceil\frac{\log\sqrt{n}\;\frac{\|x^{\mbox{\tiny{LS}}}\|_{A}}{\sigma}}{\log(1/\rho)}\rceil$ iterations using m=\frac{c_{0}}{\rho^{2}}s\log\big{(}\frac{ed}{s}\big{)} projections per round. Then with probability greater than $1-c_{1}\,{N}\,e^{-c_{2}\frac{m\rho^{2}}{\log^{4}(d)}}$, the output $\widehat{x}$ satisfies the bounds

In order to verify the predicted bound (31) on the error $\|\widehat{x}-x^{\mbox{\tiny{LS}}}\|_{A}$, we performed a second experiment. Fixing n=100s\log\big{(}\frac{ed}{s}\big{)}. we generated $T=20$ random least squares problems (as described above) with the regression dimension ranging as $d\in\{32,64,128,256\}$, and sparsity $s=\lceil 2\sqrt{d}\rceil$. Based on these choices, the least-squares estimate should have error \|x^{\mbox{\tiny{LS}}}-x^{*}\|_{A}\approx\sqrt{\frac{\sigma^{2}s\log\big{(}\frac{ed}{s}\big{)}}{n}}=0.1 with high probability, independently of the pair $(s,d)$. This predicted behavior is confirmed by the blue bars in Figure 5; the bar height corresponds to the average over $T=20$ trials, with the standard errors also marked.

On these same problem instances, we also ran the IHS algorithm using ${N}=4$ iterations with a sketch size m=4s\log\big{(}\frac{ed}{s}\big{)}. Together with our earlier calculation of $\|x^{\mbox{\tiny{LS}}}-x^{*}\|_{A}$, Corollary 2 implies that with high probability, the sketched solution $\widehat{x}=x^{{N}}$ satisfies the error bound

for some constant $c_{0}\in(1,2]$. This prediction is confirmed by the green bars in Figure 5, showing that $\|\widehat{x}-x^{*}\|_{A}\succsim 0.11$ across all dimensions. Finally, the green bars in Figure 5 show the error based on using the naive sketch estimate with a total of $M={N}m$ random projections in total; as with the case of ordinary least-squares, the resulting error is roughly twice as large. We also note that a similar bound also applies to problems where a parameter constrained to unit simplex is estimated, e.g., in portfolio analysis and density estimation .

3 Matrix estimation with nuclear norm constraints

We now turn to the study of nuclear-norm constrained form of least-squares matrix regression. This class of problems has proven useful in many different application areas, among them matrix completion, collaborative filtering, multi-task learning and control theory (e.g., ). In particular, let us consider the convex program

where $R>0$ is a user-defined radius as a regularization parameter.

Suppose that we run the IHS algorithm for ${N}=1+\lceil\frac{\log\sqrt{n}\;\frac{\|X^{\mbox{\tiny{LS}}}\|_{A}}{\sigma}}{\log(1/\rho)}\rceil$ iterations using m={c_{0}}{\rho^{2}}r\big{(}d_{1}+d_{2}\big{)} projections per round. Then with probability greater than $1-c_{1}\,{N}\,e^{-c_{2}\frac{m\rho^{2}}{\log^{4}(d_{1}d_{2})}}$, the output $X^{N}$ satisfies the bound

We have also performed simulations for low-rank matrix estimation, and observed that the IHS algorithm exhibits convergence behavior qualitatively similar to that shown in Figures 3 and 5. Similarly, panel (a) of Figure 7 compares the performance of the IHS and classical methods for sketching the optimal solution over a range of row sizes $n$. As with the unconstrained least-squares results from Figure 1, the classical sketch is very poor compared to the original solution whereas the IHS algorithm exhibits near optimal performance.

3.2 Application to multi-task learning

To conclude, let us illustrate the use of the IHS algorithm in speeding up the training of a classifier for facial expressions. In particular, suppose that our goal is to separate a collection of facial images into different groups, corresponding either to distinct individuals or to different facial expressions. One approach would be to learn a different linear classifier ($a\mapsto\langle a,\,x\rangle$) for each separate task, but since the classification problems are so closely related, the optimal classifiers are likely to share structure. One way of capturing this shared structure is by concatenating all the different linear classifiers into a matrix, and then estimating this matrix in conjunction with a nuclear norm penalty .

More specifically, in order to verify the classification accuracy of the classifier obtained by IHT algorithm, we solved the original convex program, the classical sketch based on ROS sketches of dimension $m=100$, and also the corresponding IHS algorithm using ROS sketches of size $20$ in each of $5$ iterations. In this way, both the classical and IHS procedures use the same total number of sketches, making for a fair comparison. We repeated each of these three procedures for all choices of the radius $R\in\{1,2,3,\ldots,12\}$, and then applied the resulting classifiers to classify images in the test dataset. For each of the three procedures, we calculated the classification error rate, defined as the total number of mis-classified images divided by $n_{\mbox{\footnotesize{test}}}\times d_{\mbox{\footnotesize{task}}}$. Panel (b) of Figure 7 plots the resulting classification errors versus the regularization parameter. The error bars correspond to one standard deviation calculated over the randomness in generating sketching matrices. The plots show that the IHS algorithm yields classifiers with performance close to that given by the original solution over a range of regularizer parameters, and is superior to the classification sketch. The error bars also show that the IHS algorithm has less variability in its outputs than the classical sketch.

Discussion

In this paper, we focused on the problem of solution approximation (as opposed to cost approximation) for a broad class of constrained least-squares problem. We began by showing that the classical sketching methods are sub-optimal, from an information-theoretic point of view, for the purposes of solution approximation. We then proposed a novel iterative scheme, known as the iterative Hessian sketch, for deriving $\varepsilon$-accurate solution approximations. We proved a general theorem on the properties of this algorithm, showing that the sketch dimension per iteration need grow only proportionally to the statistical dimension of the optimal solution, as measured by the Gaussian width of the tangent cone at the optimum. By taking $\log(1/\varepsilon)$ iterations, the IHS algorithm is guaranteed to return an $\varepsilon$-accurate solution approximation with exponentially high probability.

In addition to these theoretical results, we also provided empirical evaluations that reveal the sub-optimality of the classical sketch, and show that the IHS algorithm produces near-optimal estimators. Finally, we applied our methods to a problem of facial expression using a multi-task learning model applied to the JAFFE face database. We showed that IHS algorithm applied to a nuclear-norm constrained program produces classifiers with considerably better classification accuracy compared to the naive sketch.

Both authors were partially supported by Office of Naval Research MURI grant N00014-11-1-0688, and National Science Foundation Grants CIF-31712-23800 and DMS-1107000. In addition, MP was supported by a Microsoft Research Fellowship.

Appendix A Proof of lower bounds

This appendix is devoted to the verification of condition (9) for different model classes, followed by the proof of Theorem 1.

We verify the condition for three different types of sketches.

showing that condition (9) holds with $\eta=1$.

showing that the condition holds with $\eta=1$.

Finally, suppose that we sample $m$ rows independently using a distribution $\{p_{j}\}_{j=1}^{n}$ on the rows of the data matrix that is $\alpha$-balanced (7). Letting $\mathcal{R}\subseteq\{1,2,\ldots,n\}$ be the subset of rows that are sampled, and let $N_{j}$ be the number of times each row is sampled. We then have

where $p_{\infty}=\max_{j\in[n]}p_{j}$. Consequently, as long as the row weights are $\alpha$-balanced (7) so that $p_{\infty}\leq\frac{\alpha}{n}$, we have

showing that condition (9) holds with $\eta=\alpha$, as claimed.

A.2 Proof of Theorem 1

Let $\{z^{j}\}_{j=1}^{M}$ be a $1/2$-packing of $\mathcal{C}_{0}$ in the semi-norm $\|\cdot\|_{A}$ and for a fixed $\delta\in(0,1/4)$, define $x^{j}=4\delta z^{j}$. We thus obtain a collection of vectors in $\mathcal{C}_{0}$ such that

Consider the multiway testing problem of determining the index $J$ based on observing $\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{5.67886pt}{0.43057pt}}{Y}$. With this set-up, a standard reduction in statistical minimax (e.g., ) implies that, for any estimator $x^{\dagger}$, the worst-case mean-squared error is lower bounded as

where the infimum ranges over all testing functions $\psi$. Consequently, it suffices to show that the testing error is lower bounded by $1/2$.

In order to do so, we first apply Fano’s inequality conditionally on the sketching matrix $S$ to see that

Computing the KL divergence for Gaussian vectors yields

where the final inequality uses the fact that $\|x^{j}-x^{k}\|_{A}\leq 8\delta$ for all pairs.

Combined with our previous bounds (35) and (36), we find that

Setting $\delta=\frac{\sigma^{2}\log(M/2)}{64\,\eta\,m}$ yields the lower bound (10).

Appendix B Proof of Proposition 1

Since $\widehat{x}$ and $x^{\mbox{\tiny{LS}}}$ are optimal and feasible, respectively, for the Hessian sketch program (16), we have

Adding these two inequalities and performing some algebra yields the basic inequality

Since $Ax^{\mbox{\tiny{LS}}}$ is independent of the sketching matrix and $A\Delta\in\mathcal{K}^{\mbox{\tiny{LS}}}$, we have

using the definitions (18a) and (18b) of the random variables $Z_{1}$ and $Z_{2}$ respectively. Combining the pieces yields the claim.

Appendix C Proof of Theorem 2

It suffices to show that, for each iteration $t=0,1,2,\ldots$, we have

The claimed bounds (26a) and (26b) then follow by applying the bound (39) successively to iterates $1$ through ${N}$.

For simplicity in notation, we abbreviate $S^{t+1}$ to $S$ and $x^{t+1}$ to $\widehat{x}$. Define the error vector $\Delta=\widehat{x}-x^{\mbox{\tiny{LS}}}$. With some simple algebra, the optimization problem (24) that underlies the update $t+1$ can be re-written as

where \widetilde{y}:\,=y+\Big{[}I-\frac{S^{T}S}{m}\Big{]}Ax^{t}. Since $\widehat{x}$ and $x^{\mbox{\tiny{LS}}}$ are optimal and feasible respectively, the usual first-order optimality conditions imply that

As before, since $x^{\mbox{\tiny{LS}}}$ is optimal for the original program, we have

Adding together these two inequalities and introducing the shorthand $\Delta=\widehat{x}-x^{\mbox{\tiny{LS}}}$ yields

Note that the vector $A(x^{\mbox{\tiny{LS}}}-x^{t})$ is independent of the randomness in the sketch matrix $S^{t+1}$. Moreover, the vector $A\Delta$ belongs to the cone $\mathcal{K}$, so that by the definition of $Z_{2}(S^{t+1})$, we have

Combining the two bounds (41a) and (41b) with the earlier bound (40) yields the claim (39).

Appendix D Maximum likelihood estimator and examples

In this section, we a general upper bound on the error of the constrained least-squares estimate. We then use it (and other results) to work through the calculations underlying Examples 1 through 3 from Section 2.2.

The accuracy of $x^{\mbox{\tiny{LS}}}$ as an estimate of $x^{*}$ depends on the “size” of the star-shaped set

The following result bounds the mean-squared error associated with the constrained least-squares estimate:

For any set $\mathcal{C}$ containing $x^{*}$, the constrained least-squares estimate (1) has mean-squared error upper bounded as

We provide the proof of this claim in Section D.3.

D.2 Detailed calculations for illustrative examples

In this appendix, we collect together the details of calculations used in our illustrative examples from Section 2.2. In all cases, we make use tof the convenient shorthand $\widetilde{A}=A/\sqrt{n}$.

By definition of the Gaussian width, we have

since the vector $\widetilde{A}(x-x^{*})$ belongs to a subspace of dimension $\operatorname{rank}(A)=d$. The claimed upper bound (11a) thus follows as a consequence of Proposition 2.

D.2.2 Sparse vectors: Example 2

The RIP property of order $8s$ implies that

a fact which we use throughout the proof. By definition of the Gaussian width, we have

Consequently, by the Sudakov-Fernique comparison , we have

where the final inequality standard results on Gaussian widths . All together, we conclude that

Combined with Proposition 2, the claimed upper bound (12a) follows.

In the other direction, a straightforward argument (e.g., ) shows that there is a universal constant $c>0$ such that \log M_{1/2}\geq c\,s\log\big{(}\frac{ed}{s}\big{)}, so that the stated lower bound follows from Theorem 1.

D.2.3 Low rank matrices: Example 3:

By definition of the Gaussian width, we have width, we have

Thus, by duality between the nuclear and operator norms, we have

so that the upper bound (15a) follows from Proposition 2.

D.3 Proof of Proposition 2

Throughout this proof, we adopt the shorthand $\varepsilon_{n}=\varepsilon_{n}(\mathcal{K}^{*})$. Our strategy is to prove the following more general claim: for any $t\geq\varepsilon_{n}$, we have

A simple integration argument applied to this tail bound implies the claimed bound (44) on the expected mean-squared error.

Since $x^{*}$ and $x^{\mbox{\tiny{LS}}}$ are feasible and optimal, respectively, for the optimization problem (1), we have the basic inequality

Introducing the shorthand $\Delta=x^{\mbox{\tiny{LS}}}-x^{*}$ and re-arranging terms yields

where $g\sim N(0,I_{n})$ is a standard normal vector.

For a given $u\geq\varepsilon_{n}$, define the “bad” event

The following lemma controls the probability of this event:

Returning to prove this lemma momentarily, let us prove the bound (45). For any $t\geq\varepsilon_{n}$, we can apply Lemma 2 with $u=\sqrt{t\varepsilon_{n}}$ to find that

If $\|\Delta\|_{A}<\sqrt{t\,\varepsilon_{n}}$, then the claim is immediate. Otherwise, we have $\|\Delta\|_{A}\geq\sqrt{t\,\varepsilon_{n}}$. Since $\Delta\in\mathcal{C}-x^{*}$, we may condition on $\mathcal{B}^{c}(\sqrt{t\varepsilon_{n}})$ so as to obtain the bound

Combined with the basic inequality (46), we see that

a bound that holds with probability greater than $1-e^{-\frac{nt\varepsilon_{n}}{2\sigma^{2}}}$ as claimed.

It remains to prove Lemma 2. Our proof involves the auxiliary random variable

We first claim that $\mathcal{B}(u)\subseteq\{V_{n}(u)\geq 2u^{2}\}$. Indeed, if $\mathcal{B}(u)$ occurs, then there exists some $z\in\mathcal{C}-x^{*}$ with $\|z\|_{A}\geq u$ and

Define the rescaled vector $\widetilde{z}=\frac{u}{\|z\|_{A}}z$. Since $z\in\mathcal{C}-x^{*}$ and $\frac{u}{\|z\|_{A}}\leq 1$, the vector $\widetilde{z}\in\operatorname{star}(\mathcal{C}-x^{*})$. Moreover, by construction, we have $\|\widetilde{z}\|_{A}=u$. When the inequality (47) holds, the vector $\widetilde{z}$ thus satisfies $|\frac{\sigma}{n}\sum_{i=1}^{n}g_{i}\>(A\widetilde{z})_{i}|\geq 2u^{2}$, which certifies that $V_{n}(u)\geq 2u^{2}$, as claimed.

The final step is to control the probabability of the event $\{V_{n}(u)\geq 2u^{2}\}$. Viewed as a function of the standard Gaussian vector $(g_{1},\ldots,g_{n})$, it is easy to see that $V_{n}(u)$ is Lipschitz with constant $L=\frac{\sigma u}{\sqrt{n}}$. Consequently, by concentration of measure for Lipschitz Gaussian functions, we have

References