{"id":1896,"date":"2024-11-04T03:56:55","date_gmt":"2024-11-04T03:56:55","guid":{"rendered":"https:\/\/www.ethanepperly.com\/?p=1896"},"modified":"2024-12-24T17:50:05","modified_gmt":"2024-12-24T17:50:05","slug":"note-to-self-hanson-wright-and-trace-estimation-with-random-vectors-on-the-sphere","status":"publish","type":"post","link":"https:\/\/www.ethanepperly.com\/index.php\/2024\/11\/04\/note-to-self-hanson-wright-and-trace-estimation-with-random-vectors-on-the-sphere\/","title":{"rendered":"Note to Self: Hanson\u2013Wright and Trace Estimation with Random Vectors on the Sphere"},"content":{"rendered":"\n

Tyler Chen<\/a> has just released an excellent monograph<\/a> on the Lanczos method<\/a>, one of the most powerful algorithms for extracting information from a symmetric matrix \"A\". I would definitely recommend checking it out.<\/p>\n\n\n\n

One of the subjects touched on by this book is stochastic trace estimation<\/a>, a topic I\u2019ve written about several<\/a> times<\/a> before<\/a> on this blog. In his book, Tyler uses stochastic trace estimates based on random vectors on the sphere (a great choice<\/a>), and he provides a simple, but often pessimistic, mathematical analysis of this approach. To complement Tyler\u2019s amazing book, I felt it was a good opportunity to share a sharper analysis that I came up with using a more sophisticated argument. I would not be surprised if this argument appears elsewhere in the literature, but I am not aware of a source.<\/p>\n\n\n\n

Let \"x,x_1,\ldots,x_m\" be random vectors in \"\real^n\" drawn uniformly at random from the sphere of radius \"\sqrt{n}\",1<\/a><\/sup>Informally, we define the uniform distribution on the sphere to be such that the probability that \"x\" belongs to a set \"U\" is the ratio of the (hyper-)surface area of \"U\" to the (hyper-)surface area of the full sphere. Formally, the uniform distribution on the sphere can be defined multiple ways. First, treating the sphere as a subset of \"\real^n\", we can define the uniform distribution as the \"(n-1)\"-dimensional Hausdorff measure<\/a>, normalized so that the sphere has unit measure. Second, observe that the sphere is acted upon in a natural way by the group of orthogonal matrices. The uniform distribution is the unique probability measure which is invariant to this group action.<\/span> let \"A\" be a real symmetric matrix, and define the quadratic form

  <\/span>   <\/span>\"\[q \coloneqq x^\top Ax\]\"<\/p>and the Girard\u2013Hutchinson trace estimator<\/a>

  <\/span>   <\/span>\"\[\hat{\tr}\coloneqq \frac{1}{m}\sum_{i=1}^m x_i^\top Ax_i.\]\"<\/p>Both \"q\" and \"\hat{\tr}\" are unbiased<\/em> estimates for the trace of \"A\", \"\expect[q] = \expect[\hat{\tr}] = \tr A\". The goal of this post will be to bound the probability of these quantities being much smaller or larger than the trace of \"A\".<\/p>\n\n\n\n

Hanson\u2013Wright on the Sphere<\/h2>\n\n\n\n

Observe that \"\hat{\tr}\" is an average of \"m\" independent copies of the random variable \"q\". Therefore, we will begin our analysis with the quadratic form \"q\" and discuss the Girard\u2013Hutchinson estimate \"\hat{\tr}\" at the end of the post.<\/p>\n\n\n\n

Centering<\/h3>\n\n\n\n

Begin by introducing the centered<\/em> matrix

  <\/span>   <\/span>\"\[\overline{A} = A - \frac{\tr A}{n}I.\]\"<\/p>The transformation \"A\mapsto\overline{A}\" has the effect of shifting \"A\"\u2019s eigenvalues to have mean zero. Consequently<\/a>, since the trace is the sum of the eigenvalues, \"\overline{A}\" has trace zero.<\/p>\n\n\n\n

Now, rewrite \"q\" in terms of the centered matrix \"\overline{A}\":

  <\/span>   <\/span>\"\[q = x^\top A x = x^\top \overline{A} x + \frac{\tr A}{n} \cdot x^\top I x = x^\top \overline{A} x + \tr A.\]\"<\/p>In the final equality, we use the fact that \"x\" is on the sphere of radius \"\sqrt{n}\" so that \"x^\top I x = \norm{x}^2 = n\". Rearranging, we see that the error \"\overline{q}\coloneqq q-\tr A\" satisfies

  <\/span>   <\/span>\"\[\overline{q} = q-\tr A = x^\top \overline{A}x.\]\"<\/p><\/p>\n\n\n\n

We have shown that the error \"\overline{q}\" depends only on the centered<\/em> matrix \"\overline{A}\" rather than the original matrix \"A\". This observation is at the root of the reason that quadratic forms with vectors on the sphere have smaller tail probabilities than other distributions like the Gaussian distribution. Indeed, \"\overline{A}\" is smaller than \"A\" when measured using the Frobenius norm<\/a> \"\norm{\overline A}_{\rm F} \le \norm{A}_{\rm F}\", and the spectral norm<\/a> is never much larger \"\norm{\overline{A}} \le 2\norm{A}\". These observations will be important later.<\/p>\n\n\n\n

Laplace Transform Method and Comparison to Standard Gaussian Vectors<\/h3>\n\n\n\n

To derive exponential concentration inequalities, we make use of the Laplace transform method. For a random variable \"z\", define the cumulant generating function (cgf)

  <\/span>   <\/span>\"\[\xi_z(\theta) = \log (\expect [\exp(\theta z)]).\]\"<\/p>Bounds on the cgf immediately yield bounds on the random variable taking large values, so it will behoove us to estimate this quantity for \"\overline{q}\".<\/p>\n\n\n\n

To estimate the cgf of \"\overline{q}\", we will make use of a comparison between Gaussian vectors and random vectors on the sphere. Recall that a standard Gaussian vector<\/a> \"g\" is a vector with independent standard Gaussian<\/a> entries. Standard Gaussian random vectors are spherically symmetric, meaning they can be written in the form \"g = a\cdot x\", where \"x\" is a uniform random vector on the sphere and \"a\" is a scaling factor.2<\/a><\/sup>Concretely, \"a\sqrt{n}\" is a \"\chi\" random variable<\/a> with \"n\" degrees of freedom.<\/span> The scaling factor \"a\" is independent of \"x\" and, in this case, satisfies \"\expect[a^2] = 1\".3<\/a><\/sup>Indeed, the expected squared length of a standard Gaussian random vector is \"\expect[\norm{g}^2] = \sum_{i=1}^n \expect[g_i^2] = n\". But also, \"\expect[\norm{g}^2] = \expect[a^2] \cdot \expect[\norm{x}^2] = \expect[a^2] \cdot n\". Therefore, \"\expect[a^2]=1\".<\/span>\n\n\n\n

Using this relationship, the error \"\overline{q}\" can be connected to the standard Gaussian random vector \"g\" as follows

  <\/span>   <\/span>\"\[g^\top \overline{A} g = a^2 \cdot x^\top \overline{A} x = a^2 \cdot \overline{q}.\]\"<\/p>Intuitively, we might expect that multiplying the random error \"\overline{q}\" by a scaling factor of mean one should only have the effect of increasing<\/em> the variability of \"\overline{q}\" and thus increasing its probability of taking large values. Based on this intuition, we might conjecture that the cgfs are related

  <\/span>   <\/span>\"\[\xi_{\overline{q}}(\theta) \le \xi_{g^\top A g}(\theta) \quad \text{for all } \theta \in \real.\]\"<\/p>Indeed, this conjecture is true.<\/p>\n\n\n\n

\n

Proposition<\/strong> (Scaling hurts). Let \"z\" be a random variable and let \"b\" be a random variable with expectation \"1\", independent of \"z\". Then

  <\/span>   <\/span>\"\[\xi_{z}(\theta) \le \xi_{bz}(\theta) \quad \text{for all } \theta \in \real.\]\"<\/p><\/p>\n<\/blockquote>\n\n\n\n

To prove this result, let \"\expect_b\" and \"\expect_z\" denote expectations take over the randomness in \"b\" and \"z\" separately. The total expectation is<\/a> \"\expect[\cdot] = \expect_z[\expect_b[\cdot ]]\". Begin with the right-hand side and invoke Jensen’s inequality<\/a> over \"b\":

  <\/span>   <\/span>\"\[\xi_{bz}(\theta) = \xi_z(\theta) = \log (\expect_z [\expect_b [\exp(\theta bz)]]) \ge \log (\expect_z [\exp(\theta \expect_b[b]z)]) = \xi_z(\theta).\]\"<\/p>In the last line, we use the hypothesis \"\expect[b] = 1\" and the definition of the cgf.<\/p>\n\n\n\n

Having established this proposition, we conclude that \"\xi_{\overline{q}}(\theta) \le \xi_{g^\top \overline{A} g}(\theta)\" for all \"\theta\".<\/p>\n\n\n\n

Completing the Analysis<\/h3>\n\n\n\n

Finally, invoking a standard bound for the cgf of \"g^\top \overline{A}g\" for a trace-zero matrix \"\overline{A}\", we obtain

  <\/span>   <\/span>\"\[\xi_{\overline{q}}(\theta) \le \xi_{g^\top \overline{A} g}(\theta) \le \frac{\theta^2 \norm{\overline{A}}_{\rm F}^2}{1 - 2\theta \norm{\overline{A}}}.\]\"<\/p>A cgf bound of this form immediately yields the following tail bound (see Boucheron<\/a>, Lugosi<\/a>, and Massart\u2019s<\/a> Concentration Inequalities<\/em><\/a> page 29 and Exercise 2.8):

  <\/span>   <\/span>\"\[\prob \{ x^\top A x - \tr(A) \ge t \} \le \exp \left( - \frac{t^2/2}{2 \norm{\overline{A}}_{\rm F}^2 + 2 \overline{\norm{\overline{A}}t}} \right).\]\"<\/p>One can also control the lower tail event \"x^\top A x - \tr(A) \le -t\" by instantiating this result with \"-A\". Union bounding<\/a> over the upper and lower tails gives the symmetric bound

  <\/span>   <\/span>\"\[\prob \{ |x^\top A x - \tr(A)| \ge t \} \le 2\exp \left( - \frac{t^2/2}{2 \norm{\overline{A}}_{\rm F}^2 + 2 \overline{\norm{\overline{A}}t}} \right).\]\"<\/p><\/p>\n\n\n\n

Is This Bound Good?<\/h3>\n\n\n\n

Tail bounds for quadratic forms such as \"x^\top A x\" or \"g^\top A g\" are known as Hanson\u2013Wright inequalities<\/em><\/a>. For vectors on the sphere, our bound shows the tail probabilities decay like \"\exp(-O(t^2/\norm{\overline{A}_{\rm F}}^2))\" for small \"t > 0\" and \"\exp(-O(t/\norm{\overline{A}}))\" for large \"t>0\". This pattern of results “subgaussian scaling for small \"t\", subexponential scaling for large \"t\"” is typical for Hanson\u2013Wright inequalities.<\/p>\n\n\n\n

The Hanson\u2013Wright inequalities for vectors on the sphere exhibit faster decay rate than for standard Gaussian vectors. Indeed, a standard Gaussian random vector \"g\" obeys a Hanson\u2013Wright inequality of the form

  <\/span>   <\/span>\"\[\prob \{ g^\top A g - \tr A \ge t \} \le \exp \left( -\frac{t^2/2}{C\norm{A}_{\rm F}^2 + c t \norm{A}}\right)\]\"<\/p>for positive constants \"C,c > 0\". For vectors on the sphere, \"\norm{A}_{\rm F}^2\" and \"\norm{A}\" have been replaced by \"\norm{\overline{A}}_{\rm F}^2\" and \"\norm{\overline{A}}\" which are always smaller (and sometimes much smaller). The smaller tail probabilities for \"x^\top A x\" versus \"g^\top A g\" is a real phenomenon, not an artifact of the mathematical analysis.<\/p>\n\n\n\n

As another way to evaluate the quality of these bounds, recall that a Gaussian random variable with mean zero and variance \"v\" has tail probabilities roughly of size \"\exp(-t^2/2v)\". For small \"t\", our Hanson\u2013Wright inequality for vectors on the sphere has the same form with \"v = 2\norm{\overline{A}}_{\rm F}^2\". The true variance of \"q = x^\top A x\" is

  <\/span>   <\/span>\"\[\Var(q) = \frac{n}{n+2} \cdot 2\norm{\overline{A}}_{\rm F}^2,\]\"<\/p>which is close to \"v = 2\norm{\overline{A}}_{\rm F}^2\" for large \"n\". (See this post of mine<\/a> for a derivation of the variance formula.) Thus, even up to the constants, the Hanson\u2013Wright inequality we’ve derived for random vectors on the sphere is nearly as good as one could possibly hope (at least for small \"t\").<\/p>\n\n\n\n

Girard\u2013Hutchinson Estimator on the Sphere<\/h2>\n\n\n\n

Using our results for one quadratic form, tail bounds for the Girard\u2013Hutchinson estimator \"\overline{\tr} = m^{-1} \sum_{i=1}^m x_i^\top A x_i\" easily follow. Indeed, the cgf is additive for independent random variables and satisfies the identity \"\xi_{cz}(\theta) = \xi_z(c\theta)\" for constant $c, so

  <\/span>   <\/span>\"\[\xi_{\overline{\tr}}(\theta) = \sum_{i=1}^m \xi_{x_i^\top Ax_i}(\theta/m) \le m \cdot \frac{(\theta/m)^2 \norm{\overline{A}}_{\rm F}^2}{1 - 2(\theta/m)\norm{\overline{A}}}.\]\"<\/p>This cgf bound leads to tail bounds<\/p>\n\n\n\n

  <\/span>   <\/span>\"\begin{align*}\prob \{ \overline{\tr} - \tr(A) \ge t \} &\le \exp \left( - \frac{mt^2/2}{2\norm{\overline{A}}_{\rm F}^2 + 2m \overline{\norm{\overline{A}}t}} \right), \\\prob \{ |\overline{\tr} - \tr(A)| \ge t \} &\le 2\exp \left( - \frac{mt^2/2}{2\norm{\overline{A}}_{\rm F}^2 + 2m \overline{\norm{\overline{A}}t}} \right).\end{align*}\"<\/p><\/p>\n","protected":false},"excerpt":{"rendered":"

Tyler Chen has just released an excellent monograph on the Lanczos method, one of the most powerful algorithms for extracting information from a symmetric matrix . I would definitely recommend checking it out. One of the subjects touched on by this book is stochastic trace estimation, a topic I\u2019ve written about several times before onRead more<\/a><\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[9,14],"tags":[],"class_list":["post-1896","post","type-post","status-publish","format-standard","hentry","category-note-to-self","category-trace-estimation"],"_links":{"self":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/1896","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/comments?post=1896"}],"version-history":[{"count":14,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/1896\/revisions"}],"predecessor-version":[{"id":2051,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/1896\/revisions\/2051"}],"wp:attachment":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/media?parent=1896"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/categories?post=1896"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/tags?post=1896"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}