{"id":1719,"date":"2024-01-04T01:16:58","date_gmt":"2024-01-04T01:16:58","guid":{"rendered":"https:\/\/www.ethanepperly.com\/?p=1719"},"modified":"2024-01-28T02:56:45","modified_gmt":"2024-01-28T02:56:45","slug":"how-good-can-stochastic-trace-estimates-be","status":"publish","type":"post","link":"https:\/\/www.ethanepperly.com\/index.php\/2024\/01\/04\/how-good-can-stochastic-trace-estimates-be\/","title":{"rendered":"How Good Can Stochastic Trace Estimates Be?"},"content":{"rendered":"\n

I am excited to share that our paper XTrace: Making the most of every sample in stochastic trace estimation<\/em><\/a> has been published<\/a> in the SIAM Journal on Matrix Analysis and Applications<\/a>. (See also our paper on arXiv<\/a>.) <\/p>\n\n\n\n

Spurred by this exciting news, I wanted to take the opportunity to share one of my favorite results in randomized numerical linear algebra<\/a>: a \u201cspeed limit<\/em>\u201d result of Meyer<\/a>, Musco<\/a>, Musco<\/a>, and Woodruff<\/a> that establishes a fundamental limitation on how accurate any trace estimation algorithm can be.<\/p>\n\n\n\n

Let\u2019s back up. Given an unknown square matrix \"A\", the ](https:\/\/www.ethanepperly.com\/index.php\/2023\/01\/26\/stochastic-trace-estimation\/”>trace estimation problem<\/a> is to estimate the trace<\/a> of \"A\", defined to be the sum of its diagonal entries

  <\/span>   <\/span>\"\[\tr(A) \coloneqq \sum_{i=1}^n A_{ii}.\]\"<\/p>The catch? We assume that we can only access the matrix \"A\" through matrix\u2013vector products (affectionately known as \u201cmatvecs\u201d): Given any vector \"x\", we have access to \"Ax\". Our goal is to form an estimate \"\hat{\tr}\" that is as accurate as possible while using as few matvecs as we can get away with.<\/p>\n\n\n\n

To simplify things, let\u2019s assume the matrix \"A\" is symmetric<\/a> and positive (semi)definite<\/a>. The classical algorithm for trace estimation<\/a> is due to Gir<\/a>ard<\/a> and Hutchinson<\/a>, producing a probabilistic estimate \"\hat{\tr}\" with a small average (relative) error:

  <\/span>   <\/span>\"\[\expect\left[\frac{|\hat{\tr}-\tr(A)|}{\tr(A)}\right] \le \varepsilon \quad \text{using } m= \frac{\rm const}{\varepsilon^2} \text{ matvecs}.\]\"<\/p>If one wants high accuracy, this algorithm is expensive<\/em>. To achieve just a 1% error (\"\varepsilon=0.01\") requires roughly \"m=10,\!000\" matvecs!<\/p>\n\n\n\n

This state of affairs was greatly improved by Meyer<\/a>, Musco<\/a>, Musco<\/a>, and Woodruff<\/a>. Building upon previous<\/a> work<\/a>, they proposed the Hutch++ algorithm<\/a> and proved it outputs an estimate \"\hat{\tr}\" satisfying the following bound:

(1) <\/span>   <\/span>\"\[\expect\left[\frac{|\hat{\tr}-\tr(A)|}{\tr(A)}\right] \le \varepsilon \quad \text{using } m= \frac{\rm const}{\varepsilon} \text{ matvecs}.\]\"<\/p>Now, we only require roughly \"m=100\" matvecs to achieve 1% error! Our algorithm, XTrace<\/a>, satisfies the same error guarantee (1)<\/a> as Hutch++. On certain problems, XTrace can be quite a bit more accurate than Hutch++.<\/a><\/p>\n\n\n\n

The MMMW Trace Estimation “Speed Limit”<\/h2>\n\n\n\n

Given the dramatic improvement of Hutch++ and XTrace over Girard\u2013Hutchinson, it is natural to hope: Is there an algorithm that does even better than Hutch++ and XTrace? For instance, is there an algorithm satisfying an even slightly better error bound of the form

  <\/span>   <\/span>\"\[\expect\left[\frac{|\hat{\tr}-\tr(A)|}{\tr(A)}\right] \le \varepsilon \quad \text{using } m= \frac{\rm const}{\varepsilon^{0.999}} \text{ matvecs}?\]\"<\/p>Unfortunately not. Hutch++ and XTrace are essentially as good as it gets.<\/p>\n\n\n\n

Let’s add some fine print. Consider an algorithm for the trace estimation problem. Whenever the algorithm wants, it can present a vector \"x_i\" and receive back \"Ax_i\". The algorithm is allowed to be adaptive: It can use the matvecs \"Ax_1,\ldots,Ax_s\" it has already collected to decide which vector \"x_{s+1}\" to present next. We measure the cost of the algorithm in terms of the number of matvecs alone, and the algorithm knows nothing about the psd matrix \"A\" other what it learns from matvecs.<\/p>\n\n\n\n

One final stipulation:<\/p>\n\n\n\n

\n

Simple entries assumption.<\/strong> We assume that the entries of the vectors \"x_i\" presented by the algorithm are real numbers between \"-1\" and \"1\" with up to \"b\" digits after the decimal place.<\/p>\n<\/blockquote>\n\n\n\n

To get a feel for this simple entries assumption, suppose we set \"b=2\". Then \"(-0.92,0.17)\" would be an allowed input vector, but \"(0.232,-0.125)\" would not be (too many digits after the decimal place). Similarly, \"(18.3,2.4)\" would not be valid because its entries exceed \"1\". The simple entries assumption is reasonable as we typically represent numbers on digital computers by storing a fixed number of digits of accuracy.1<\/a><\/sup>We typically represent numbers on digital computers by floating point<\/em> numbers<\/a>, which essentially represent numbers using scientific notation like \"1.3278123 \times 10^{27}\". For this analysis of trace estimation, we use fixed point<\/em> numbers<\/a> like \"0.23218\" (no powers of ten allowed)!<\/span>\n\n\n\n

With all these stipulations, we are ready to state the “speed limit” for trace estimation proved by Meyer, Musco, Musco, and Woodruff:<\/p>\n\n\n\n

\n

Informal theorem (Meyer, Musco, Musco, Woodruff).<\/strong> Under the assumptions above, there is no trace estimation algorithm producing an estimate \"\hat{\tr}\" satisfying

  <\/span>   <\/span>\"\[\expect\left[\frac{|\hat{\tr}-\tr(A)|}{\tr(A)}\right] \le \varepsilon \quad \text{using } m= \frac{\rm const}{\varepsilon^{0.999}} \text{ matvecs}.\]\"<\/p><\/p>\n<\/blockquote>\n\n\n\n

We will see a slightly sharper version of the theorem below, but this statement captures the essence of the result.<\/p>\n\n\n\n

Communication Complexity<\/h2>\n\n\n\n

To prove the MMMW theorem, we have to take a journey to the beautiful subject of communication complexity<\/em><\/a>. The story is this. Alice and Bob are interested in solving a computational problem together. Alice has her input \"x\" and Bob has his input \"y\", and they are interested in computing a function \"f(x,y)\" of both their inputs.<\/p>\n\n\n\n

Unfortunately for the two of them, Alice and Bob are separated by a great distance, and can only communicate by sending single bits (0 or 1) of information over a slow network connection. Every bit of communication is costly. The field of communication complexity<\/strong> is dedicated to determining how efficiently Alice and Bob are able to solve problems of this form.<\/p>\n\n\n\n

The Gap-Hamming problem<\/a> is one example of a problem studied in communication complexity. As inputs, Alice and Bob receive vectors \"x,y \in \{\pm 1\}^n\" with \"+1\" and \"-1\" entries from a third party Eve. Eve promises Alice and Bob that their vectors \"x\" and \"y\" satisfy one of two conditions:

(2) <\/span>   <\/span>\"\[\text{Case 0: } x^\top y \ge\sqrt{n} \quad \text{or} \quad \text{Case 1: } x^\top y \le -\sqrt{n}. \]\"<\/p>Alice and Bob must work together, sending as few bits of communication as possible, to determine which case they are in.<\/p>\n\n\n\n

There’s one simple solution to this problem: First, Bob sends his whole input vector \"y\" to Alice. Each entry of \"y\" takes one of the two value \"\pm 1\" and can therefore be communicated in a single bit. Having received \"y\", Alice computes \"x^\top y\", determines whether they are in case 0 or case 1, and sends Bob a single bit to communicate the answer. This procedure requires \"n+1\" bits of communication.<\/p>\n\n\n\n

Can Alice and Bob still solve this problem with many fewer than \"n\" bits of communication, say \"\sqrt{n}\" bits? Unfortunately not. The following theorem of Chakrabati and Regev<\/a> shows that roughly \"n\" bits of communication are needed to solve this problem:<\/p>\n\n\n\n

\n

Theorem (Chakrabati\u2013Regev).<\/strong> Any algorithm which solves the Gap-Hamming problem that succeeds with at least \"2/3\" probability for every pair of inputs \"x\" and \"y\" (satisfying one of the conditions (2)) must take \"\Omega(n)\" bits of communication.<\/p>\n<\/blockquote>\n\n\n\n

Here, \"\Omega(n)\" is big-Omega<\/a> notation, closely related to big-O notation<\/a> \"\order(n)\" and big-Theta notation<\/a> \"\Theta(n)\". For the less familiar, it can be helpful to interpret \"\Omega(n)\", \"\order(n)\", and \"\Theta(n)\" as all standing for \u201cproportional to \"n\"\u201d. In plain language, the theorem of Chakrabati and Regev result states that there is no algorithm for the Gap-Hamming problem that much more effective than the basic algorithm where Bob sends his whole input to Alice (in the sense of requiring less than \"\order(n)\" bits of communication).<\/p>\n\n\n\n

Reducing Gap-Hamming to Trace Estimation<\/h2>\n\n\n\n

This whole state of affairs is very sad for Alice and Bob, but what does it have to do with trace estimation? Remarkably, we can use hardness of the Gap-Hamming problem to show there\u2019s no algorithm that fundamentally improves on Hutch++ and XTrace. The argument goes something like this:<\/p>\n\n\n\n

    \n
  1. If there were a trace estimation algorithm fundamentally better than Hutch++ and XTrace, we could use it to solve Gap-Hamming in fewer than \"\order(n)\" bits of communication.<\/li>\n\n\n\n
  2. But no algorithm can solve Gap-Hamming in fewer than \"\order(n)\" bits or communication.<\/li>\n\n\n\n
  3. Therefore, no trace estimation algorithm is fundamentally better than Hutch++ and XTrace.<\/li>\n<\/ol>\n\n\n\n

    Step 2 is the work of Chakrabati and Regev, and step 3 follows logically from 1 and 2. Therefore, we are left to complete step 1 of the argument.<\/p>\n\n\n\n

    Protocol<\/h3>\n\n\n\n

    Assume we have access to a really good trace estimation algorithm. We will use it to solve the Gap-Hamming problem. For simplicity, assume \"n\" is a perfect square. The basic idea is this:<\/p>\n\n\n\n