Hardness of Low Rank Approximation of Entrywise Transformed Matrix Products
Authors: Tamas Sarlos
Xingyou Song
David Woodruff
Qiuyi Zhang
Venue: NeurIPS 2023 poster (2023)
Keywords: Low rank approximation
kernel methods
fine-grained complexity
TL;DR: We show new conditional hardness results for low rank approximation of entrywise transformed matrix products.

Abstract: Inspired by fast algorithms in natural language processing, we study low rank approximation in the entrywise transformed setting where we want to find a good rank $k$ approximation to $f(U \cdot V)$, where $U, V^\top \in \mathbb{R}^{n \times r}$ are given, $r = O(\log(n))$, and $f(x)$ is a general scalar function. Previous work in sublinear low rank approximation has shown that if both (1) $U = V^\top$ and (2) $f(x)$ is a PSD kernel function, then there is an $O(nk^{\omega-1})$ time constant relative error approximation algorithm, where $\omega \approx 2.376$ is the exponent of matrix multiplication. We give the first conditional time hardness results for this problem, demonstrating that both conditions (1) and (2) are in fact necessary for getting better than $n^{2-o(1)}$ time for a relative error low rank approximation for a wide class of functions. We give novel reductions from the Strong Exponential Time Hypothesis (SETH) that rely on lower bounding the leverage scores of flat sparse vectors and hold even when the rank of the transformed matrix $f(UV)$ and the target rank are $n^{o(1)}$, and when $U = V^\top$. Furthermore, even when $f(x) = x^p$ is a simple polynomial, we give runtime lower bounds in the case when $U \neq V^\top$ of the form $\Omega(\min(n^{2-o(1)}, \Omega(2^p)))$. Lastly, we demonstrate that our lower bounds are tight by giving an $O(n \cdot \text{poly}(k, 2^p, 1/\epsilon))$ time relative error approximation algorithm and a fast $O(n \cdot \text{poly}(k, p, 1/\epsilon))$ additive error approximation using fast tensor-based sketching. Additionally, since our low rank algorithms rely on matrix-vector product subroutines, our lower bounds extend to show that computing $f(UV)W$, for even a small matrix $W$, requires $\Omega(n^{2-o(1)})$ time.

Hardness of Low Rank Approximation of Entrywise Transformed Matrix Products
Tamas Sarlos Google Research
Xingyou Song Google Deepmind xingyousongggooglee.. com
David P. Woodruff Carnegie Mellon University dwoodruf@cs, cmu edu
stamas@google.com com
Qiuyi (Richard) Zhang Google Deepmind qiuyizogoogle.c com
Abstract
Inspired by fast algorithms natural language processing, we study low rank approximation in the entrywise transformed setting where we want to find good rank k approximation to fUUVV., where U, V E _ are given, r = (log(n)), and f(x) is a general scalar function. Previous work in sublinear low rank ap- proximation has shown that if both (1) U = VT and (2) f(x) is a PSD kernel function, then there is an Onk"~11 time constant relative error approximation algorithm, where z 2.376 is the exponent of matrix multiplication. We give the first conditional time hardness results for this problem, demonstrating that both conditions (1) and (2) are fact necessary for getting better than n² n=>-((1 time for relative еттоr low rank approximation for a wide class of functions. We give novel reductions from the Strong Exponential Time Hypothesis (SETH) that rely on lower bounding the leverage scores of flat sparse vectors and hold even when the rank of the transformed matrix f(UV) and the target rank are n (1) and when U = VT_ Furthermore, even when f(x) xº is a simple polynomial, bounds in the case when U * V of the form (nmin(n2-ot) we demonstrate that our lower bounds are tight by giving an O(n poly 2 1/c)) time relative eтror approximation algorithm and fast O poly(k: (k. p. 1/e)) addi- tive еггог approximation using fast tensor- tensor-based sketching. Additionally, since our low rank algorithms rely on matrix-vecter product subroutines, our lower bounds extend to show that computing (UV),, for even small matrix W, requires 2 time.
give runtime lower (2))))) Lastly,
Introduction
The central idea behind the classic problem of low rank approximation (LRA) is to approximate a given matrix with rank & matrix that preserves the important features of the original matrix, while being computationally more efficient and statistically more stable. One particular area of interest in LRA is the low rank decomposition of entrywise transformed matrices, where the entries of the original matrix transformed through non-linear function before being approximated, with various applications in kernel methods, self-attention, and likelihood computations [Choromanski al.) [022 Levy and Goldberg] 20014). In these settings. the Gram matrix of dot-product kernel, the attention module's 's product operator, or even the non-linear computation of activation in deep network can all be represented as fUUV), where the inputs are low dimensional matrices U, V. Note that while UV is low rank, fUV) may not be, with the rank blowup dependent on the choice of the transformation f.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

For entrywise transformed low rank approximation, we aim to find an approximately optimal rank k approximation, terms of relative Frobenius norm error, to the matrix A f(UV) given U E R and V grxd, where f a scalar transformation, In this paper, our goal is study the functions such that this task is solvable in subquadratic 0n)²-9) time and for case of presentation, we assume d in this section. This problem was studied in the distributed setting in Woodruff and Zhong (2166) for functions f(x) Ixel", where the goal was to minimize communication. It was also studied in the streaming setting, where was shown that On poly(k/c)) memory suffices to solve rank approximation with additive error [(el]U [I for the function f(r) = log([æ] + 1) in a single pass Jiang et al. 202T]. improving the earlier work of [Liang et al., 2020 In a related paper, |Han et al. 2020] present an algorithm for low-rank approximation of polynomial entrywise transforms; however, they make no comparison to the optimal rank h error.
In the setting when U = VT, certain choices of f(x) can surprisingly admit subgradratic relative error algorithms, such as when f represents a positive semi-definite (PSD) kernel [Musco and Musco 2017]. Moreover, for any PSD matrix of any rank, by sampling according to the leverage scores of the matrix square root of the kernel, recent work shows that relative error low rank decomposition is possible in (nz(k/c)=-) time, where w 2 2.373 is the exponent of matrix multiplication: see [Musco and Woodruff 2017 Bakshi et al. 2020]. Note that in many cases when U = VT, carefully choosing J will result in A J(VV) being PSD. In fact, when u; are all unit norm, if f admits Taylor expansion with non-negative coefficients, then A = f(UUT) is PSD matrix and admit aub--uadratic low rank approximations.
On the surface, this astounding algorithmic result seems to conflict with subquadratic lower bounds for basic linear algebraic primitives for many kernel functions, especially the exponential kernel used in the attention architecture. Specifically, such work looks at kernel matrix vector products and shows they take (()==(1)) time for relative and additive error approximations [Keles et al. 2023]. Furthermore, these bounds can be refined to hold in the restricted regime when the entries of the matrix are (2( (loog(n) [Alman and Song 2023 J. These conditional lower bounds, as with others for linear algebraic problems, are derived from hardness based on the Strong Exponential Time Hypothesis (SETH) [Lokshtanov et al., 2013).
The crucial observation to resolve this seeming contradiction is to recognize that the low rank approximation allows for an error bound that depends on the error of the best rank-k approximation. which can be much larger that the error tolerated in the previous lower bounds. Therefore, hardness for LRA is inherently different than hardness for approximating the entire matrix additively or than matrix- vector (MV) multiplication. In fact, we will show that hardness for MV multiplication is, in some sense, strictly easier to establish.
Some related lower bounds include the work of Backurs et al. [017]] that solving kernel Support Vector Machines (SVM), ridge regression, or Principal Component Analysis (PCA) problems to high accuracy or approximating kernel density estimates up to a constant factor for kernels with exponential tails, requires 2-1) time assuming SETH. Also,|Amman et al. [2020] show that for linear algebraic primitives for the Laplacian of a graph with weights given by kernel, most operations, such as c-approximate matrix- vector multiplication, are hard, although their hardness results assume log(1) (1/c) dependence when reducing to SETH.
1.1 Our Contributions
In the setting of entrywise-transfommee low rank approximation, we show that we cannot get sub- quadratic relative error LRA generally when either 1) U + VT even for PSD kernel functions or 2) for transformations that are approximately polynomials of [ar] even for constant degree. Therefore, without multiple strong structural assumptions on A, there is no subquadratic approximation algo- rithm, even when (log(()) and the rank of the transformed matrix fLV and the target rank are (1)) We emphasize these novel hardness results hold for large class of transformations and in fact generalize to hardness for matrix vector multiplication, which can be derived as a corollary, implying that computing (UVV) also requires ((n=2001)) time for many / that have not been studied before (see Theorem 3.2). On the positive side, for the polynomial activation f(x) = =², we provide an O(n poly(r", k, 1/e)) time algorithm for relative error LRA and O(n poly(p, k, 1/e)) time for additive error LRA via fast consor-based sketches (see Algorithm 2).

---

Upper Bounds
Prior Work
Complexity
Matrix Type
Task
Musco and Woodrutf 2017
O poly (k/e))
PSD. U VT
Relative LRA
Bakshi ot al: 2020
Oin h (E/e)=-1)
PSD. U
Relative LRA
This work
( (oly(2", 1/e]]
(x) o², U #V
Relative LRA
This work
O(n polyip. R. 1.1))
(()) 5², UU
Additive RA
Lower Bounds
Backurs al. 2017
2in <
Gaussian, U V
PCA, Regression
Alman et al. 2020
(M
Kernel I Laplacions U = V
MV product
Keles et all
2023
Rin
Exponential Kernlss, U - V
MV product
This work
Qin wayg
[(x) = a + O[ P ').U= V
Relative LRA, MV
This work
D(miniaa 20,, (019))
(x) 2².0 XV
Relative LRA,
Table 1: Overview of upper bounds for low rank approximation (LRA) and lower bounds for relative linear algebraic primitives, where the dependence on r is omitted. This table illustrates multiple separation results: 1) LRA strictly easier than matrix vector (MV) products for positive semidefinite kernels, 2) LRA is strictly easier for functions f()) that are kernels, even when considering low- degree polynomials of |z]]. 3) LRA is strictly easier when U = V for f(x) = xx for = (nlog(n)), which has a 20(P) lower bound when U * VT, 4) additive error LRA A is strictly easier than relative enror LRA. Our work extends to any function (ax that admits Taylor series dominated by |ar around x = 0. For example, this includes f(x) = log([z] + 1) = [2]] + ((1z)
We note that previous lower bound techniques do not apply for LRA since relative error LRA approximations of f(VV) are possible without approximating a matrix vector product, as shown by the (n))))()) time LRA algorithm for the popular exponential kernel
[(x) = exp(2 (3) by querying sublinear number of entries of f(UUT) [Musco and Woodruff 2017 Bakshi et al. 2020). Indeed, our novel lower bounds differ from previous reductions from Orthogonal Vectors Problem (OVP) by
explicitly creating LRA instances and exploiting linear algebraic structural properties of the column spaces of a low-rank tensored matrix.
Specifically, our reduction uses an structural property that if there is a pair of input vectors to OVP which are orthogonal, then under mild assumptions, this implies we can find a sparse vector in the column span of the low rank approximation of S(UV), where U, V are matrices containing the input vectors to OVP. Also, our algorithm (Algorithm T} relies on another critical observation that this sparse vector has uniform magnitude in the non-zero entries, each of which corresponds to a vector with an orthogonal vector pair. By exploiting the structure of this vector, we can lower bound the leverage score of each row that corresponds to non-zero entry of this sparse vector, and we can appeal to fast leverage score computation algorithms to find a small (20 superset of the support of the sparse vector. Finally, we can quickly check which pairs of entries in the superset correspond to vectors in the original OVP problem via brute-force force search, completing the reduction. The precise reduction is bit more technically involved, as we also need to allow for additive error for our applications.
In the setting when U * VT, we show novel lower bounds in the case when f(x) x" is a polynomial kernel function and admits a simple rank r2 decomposition. We show that we cannot do better than the naive decomposition and prove an Q(min(n. ,2((p) lower bound (see Theorem 3.3). This implies in the setting when (log(n)) and k = (1), that surprisingly there is a separation between the (1) U = V setting and the (2) U V setting: when f(x) = d2p, in the first setting, previous works show that there admits an 0(k approximation algorithm; however, in the second setting when U & VT, our lower bounds imply that there cannot exists truly subquadratic time algorithm. We emphasize that these exponential lower bounds in the degree for the polynomial kernel are the first of their kind to rule out truly subquadratic algorithms for p = (o(log(n)). even when f(x) = x is a simple polynomial.
We also give an O(n (oly(rP/e) time algorithm for relative error approximation, which agrees with our lower bounds (see Theorem 4.1,, showing that our lower bounds for the polynomial kernel when U * VT are in fact tight. In addition, we provide an O(n poly k. 1/c)) time algorithm for additive error LRA that avoids the inherent exponential dependence on p. highlighting a separation between additive and multiplicative error LRA when
((log((n). Specifically, we can achieve the easier additive error approximation with fast tensor-based sketching matrices (see Theorem 4.2) by applying comparable techniques to an independent subsequent work on polynomial- based transformers, although the main difference is that they do not output a rank h approximation and suffers a worse dependence p, 1/c due to their non-negativity guarantees [Kacham et al.
3

---

2023]. We remark that while the polynomial tensor matrix itself can be approximated relatively in O(n-poly(p) poly(p) runtime, we emphasize that the polynomial kemel, a product of two tensor matrices, must incur additive enror or suffer (2(2)) runtime for relative emor guarantees.
O[n
Finally, we observe that we can generalize our lower bounds for LRA by exploiting the fact that our algorithms for fast LRA reduce to matrix vector product! Therefore, our relative error algorithms are reductions that give MV lower bounds, which are automatically derived for a large number of entrywise transformed matrices, even when U/ = (see Theorem 4.3). This implies that in some sense, LRA an easier problem than MV approximation, implying that our lower bounds for LRA are stronger. We summarize the most relevant prior results and our contributions in Table
1.2 Related Work to Transformers and Natural Language Processing
A primary downstream application from our work is in the field of natural language, as it is common to compute similarity matrices f(UV) from token embeddings U, V. Such is the case for the popular use of Transformers (Vaswani al. 2017] and their efficient attention variants (Tay et al. 2022]. The standard attention mechanism consists of the operation softmax EVV WW, where in our notation, U, V, W are the "query", "key". and "value" matrices respectively. Simplifying and ignoring normalizations, this can be seen f(LVV) where f(x) = exp(x). For given f. one can potentially linearize the attention mechanism if there exist U** V* such that f(VV) = UV as then the order of matrix multiplication can be rearranged into U ()) which allows memory and runtime in (ndk).
In the unnormalized softmax case where f(x y) = exp(x T y). Choromanski et al. 2021 notes that exp(x T y) EC-N(o. exp (-==11=11
|=]
2 exp ( y
which thus allows defining
U*, V* as random feature matrices from sampled <i. ........ CR- Further use of kernel properties to improve the softmax approximation have been introduced in Choromanski et al. [0022] Likhosherstov et al. [2022]. More generally, for f such that fx=Ty) y) K(x, [z, y) admits kernel representation, one may consider using variants of Bochner's s theorem [Feller 1968) to provide approximations via random Fourier features (Rahimi and Recht 2007).
While works such as Tsai et al. nomial, exponential, and RBF kernels, so far the dominant paradigm, termed the class of Linear Transformers Katharopoulos et al. (202)). is to conveniently instead consider the reverse case, where one defines mapping IR - IR in order to define the kernel K(x,, y) = <(x) Φ(v). Unfortunately, usually this does not lead to closed-form f such that f { y) Φ(x) Φ(v), which can lack interpretabillity and compatibility with classic attention mechanisms.
2019 Kacham al.
2023) have experimented with linear, poly-
However, one may consider nonlinear functions/ I which do not admit a kernel; for example, the class of functions [x) log"( ((x:] 1) for >> 0 is used in Levy and Goldberg [2144) Li et al. (2015) to compute implicit word embeddings and corresponding generative models. Within this application domain, our work thus provides answers to the question: for which classes of functions f can one efficiently compute low-rank approximations U*, V* such that F(UV) se U*V*, especially when f(x 3) does not admit a kernel structure y))
2 Preliminaries
We let our implicit matrix A f(UV) be n x d where U E RR=x and V E Brxd and n 2 d without loss of generality, where R + R is scalar function. We will also assume that d = 2(1) and it is often the case in empirical settings that n. We say that a matrix A is positive semi-definite (PSD) if it is symmetric and only has non-negative eigenvalues; a function f is a kernel function if (UUUT is PSD for any matrix U. The i-th leverage score & of matrix B E Rad is equal to its sensitivity, o² Tell where B, the i-th row of B. The Khatri-Rao product of B E xd is C' E gandd" where the i-th row of C is B, tensored with itself p times. We use the standard notation that Ö(f) is O(f poly(log f)).
R
i.e. & sup (B. 2)²
T| [Bal
sup suPpecolpan(E)
For the approximate rank & LRA problem, we want to find U' E RRnk and V Rx such that || U'VI < +c) min
||A ~Q
DEReXA,
4

---

where I- denotes the Frobenius norm, unless otherwise specified and we denote [A]x = Ü- as some rank & matrix that minimizes the objective. Note this approximation approximation guarantee, but can be analogously defined for additive error. For stronger lower bounds, we consider a weakened version of this problem given by outputting only a best rank projection, specifically we want to find orthogonal W R || - AWWT| || S (1++)]1A- [A]xil-
relative error
Our lower bounds will rely on reductions from the conditional hardness of OVP and Max-IP, whose hardness comes from SETH. We observe, from our remarks, that some more restrictive structure can be placed on the input vector sets A, B to OVP without removing the hard instances.
Assumption 1. (Hardness of Orthogonal Vectors Problem (OVP) Williams] 2005)). Let A = {ai...., On and B (b1....,ba' _ba) be sets, where a; E {0, 13" are binary vectors for all E [d]. Any algorithm which given any input (A, B3 decides with constant probability if there at least one pair of vectors
i [m] = (1.2....,,JJ.
A and b E BB such that a² 0 requires (nd))-o(1) time, provided r = w(log n) and r = no(1), Observe that the lower bounds also hold when A == B is enforced.
Remark 2. Previous work [Williams 2005 Vassilevska Williams) 2015] has shown thar Assumption I
true for - d unless the Strong Exponential Time Hypothesis (SETH, Impagliazzo and Paturi
(200Ty) is false. Given this assumption, one can handle general n and d by a padding argument: if one could solve OVP with an algorithm A running in most (nd)1-C time for a constant C > and without loss of generality n 2 d, then one could solve the problem when |A| = [B] = n by splitring B into ( d) disjoint sets each of size at most d, and solving the problem on each disjoint contradicting the assumption in
set in total time less than (nd)²-C (n//d) = n²-C/d /d S <n²-C, the |A| |B| = n case.
Remark 3. In Assumption 7 we can enforce more restrictive structure on the inpul sets A, B without making the problem easier. Specifically, we can assume that there are at most n( distinct pairs with a E A and b E B for which a a²'b=00 0. Indeed, otherwise by sampling (nd)2-S(()) (-(1 pairs at random and checking if a = O, we would solve the OVP problem in (rd)) 1-2(1) T (nd)) 12((1) time, using that r = (1), This would contrudict Assumption 7
Assumption 4. (Hardness of Apx-Max-II, of Chen and Williams| (2019]). Two sets A n.d Problem, Definition 2. Remark 2.2, and Lemma 4.1 {oi,.... an and B = (b1,.... bn} are given, where a, by E {0, are binary vectors for all E [n]. with (log(n)). Let max b. Any aEA,bEB
algorithm which outputs a number m [m/100, m] with constant probability requires (nd))-o(1) time.
Remark 5. As in Remark Land and Remark 3 we can reduce the general n and d case the case d of previous work [Chen and Williams, 2019), and we can also enforce that there are at most n (1) pairs A and b E B for which b m/11mm11m0 100, as otherwise sampling would solve the problem in less than (nl) time.
3 LRA Runtime Lower Bounds
We show that the complexity of low rank approximation of an entrywise transformed matrix turns out to heavily depend on the type of entrywise transformation. Based on standard complexity assumptions, we show an (nd))-o(1)
n 1+((1) time lower bound for entrywise function f(x) = | for odd integers in Theorem 3.1 below, while we show a weaker 20(p) lower bound for any integer in Theorem 3.3 which becomes (nd) for p = O(log n). This is no accident, as we show a matching 20()) ( +o(1) upper bound in Theorem 4.1 Furthermore, our lower bounds extend to functions that are approximately odd-degree polynomials of (ar], even if the degree is small. Specifically, this includes the commonly used [(x) = log([ar] + 1) 2 |ar]-
Our proofs use variations of the Orthogonal Vectors Problem (OVP) to create two types of instances of implicit low rank approximation to create fast (oo(+)()) time algorithm to solve OVP (see Al gorithm 1) when given access to a low rank approximation algorithm with constant relative error and small constant additive error guarantees. Let A, B be the set of input vectors of OVP. Then, in our reduction, if there is a pair of input vectors which are orthogonal, then applying LRA, we can find either 1) a column span deviation of the low rank approximation or 2) sparse vector in the column span of the low rank approximation. Note that we not look at the exact column spans, but allow an additive error depending the additive error of the output.
5

---

In the latter case, we note that our reduction structure ensures that all entries in this vector have the same magnitude. Consequently, since the column span of the output has low rank, each such entry in the support of this sparse vector corresponds to a large leverage score, and so to find the support of this unknown sparse column we can use algorithms to compute or approximate all leverage scores given the low rank factorization of the entrywise transformed matrix. This enables us to find a small superset of the support of the sparse vector, and then brute-force which pairs of entries in the superset correspond to vectors in the original OVP problem that intersect. The overall time in the reduction is negligible compared to the time to solve OVP, and therefore it must that finding the factorization of the entrywise transformed matrix itself was expensive.
In the former case, it follows that the column spans of the output are necessarily far from each other in the two cases, and we can quickly check this and therefore solve the OVP problem. Together, it follows that the time for finding the factorization of the entrywise transformed matrix itself must have been large. While there are quantitative differences in the two cases of odd and even integers the proofs both follow this strategy.
Algorithm 1 OVP to LRA Reduction
Input: Sets A = fai...., an E {0,1}*} and B = (b)..... ..., bn E {0, 1}x}
Output: YES if there exists a Tb = and NO otherwise.
I: Choose E {-1, 1) uniformly random
2: Let U [A] c| E [x*x(++1) and V = [B cj R=x((+1) and = + +1
3: Let W = LRA(U, V, = [x(P) for constant relative error and a additive error
4: Compute U = U @ @ U gaxr" U tensored with itself p times and similarly for V
5: Let e
W/U' e,112 for i-th column >Calculate column distances of
U"V" to column span of W.
6: If any of r?: 1.01a, output YES
pa is any upper bound additive error of LRA 7: Compute 1/ 1/2-approximate leverage scores & of U" and let S {i [m]le; 1/ 1/(100n(()) FFinds the representative subset and the threshold is given by the OVP assumption.
8: Compute all dot products between a, and B for all i S. Output YES if there exists pair such that 0 by = 0 and NO otherwise
3.1 General Functions
[x=| for odd integers p
In the following theorem we will also allow for tiny amount of additive error, this will be useful for our later lower bound applications where we approximate other functions using a Taylor series and reduce from the following theorem.
Theorem 3.1. Suppose n E N. f(x) = |x| for an odd integer constant = O(log(n)), and r = (l)g(n)). There is positive integer k = n1()), such that for any approximation factor A 2 1 and any constant a any possibly randomized, algorithm which given any input U Rarr and V E rxx outputs W E Bink satisfying
|I(UU V)WWT- f(U (IU-VI?ES IIF(U V)]a -(U-VV +a,
with constant probability, requires (nd)1 (o()) time, under Assumption Further, this holds even
ifU VT. Here [f -V)); denotes the best rank-k approximation to foJ JU VV) in the Frobenius norm.
Proof. Suppose A and B are input sets to the OVP problem of Assumption I Dwth parameter = where s is the dimension of the points. We let the rows of U be the points in A, but we append one additional dimension, represented column vector to U that is chosen uniformly at random in {-1,1}" Similarly, the columns of V correspond to the points in B, but we append one additional row to V, which is equal to cT. Observe that if A = B, then necessarily U VT.
Let Case 1 be when there is no E A and E B with a 0, and Case 2 be when there is an E and a B with a²') = 0, and we are deciding whether we are in Case We claim that our algorithm (Algorithm T) will always output NO when we are in Case and will output YES with constant probability when we are Case 2. This clearly suffices to solve OVP with constant probability.

---

Notice that if there is no E A and b E B for which a ab = O, then for all A and b E B, we have a b 2 1. Consequently with probability 1 over the choice of column vector c, for all i [n] and j E [d. U.V 2 O, where U; is the i-th row of U and Vj is the j-th column of V. Consequently, f(U.V)) (U:V))P. Thus, f(U V) U"V", where each row of U" E gnxr" is the Khatri-Rao product of itself p times, and each column of V E " xd is the Khatri- Rao product of itself p times.
Observe that the rank of f(U V) is at most runk + n (l))) S (s + 1)P + n since by Remark 3 we have at most n n1()) dot product pairs that are zero, implying that S(U (U--)) V) can be obtained from U" V" by changing at most no(1) entries. It follows that there is a value k v n + (s + 1)P such that, in both cases, for any multiplicative approximation factor A 2 1, necessarily the output W satisfies [I V'W wWW.wT--- - f(U V)ll S a with constant probability.
Next, we compare the column span of W to that of U" V", In Case 1, there is pairwise dot product that is zero, so we have that f(U - -V) = U"V". Therefore, each column of U" V" has squared distance at most a to the column span of W.
Note that we can compute all squared distances of U"V" to the column span of W in (n)))())) time. This so that we can detect whether we are in case I, where U" has at squared distance at most a to the column span of W, or we are in case 2. Specifically, we compute the squared distance by the Pythagorean theorem and since W is orthogonal:
= V - e:13
and for each i, we can compute this in (() time, given the matrices we have precomputed in the previous paragraph. Thus, we can compute all squared distances up to additive 1/ poly(n) in n I+o(1)) time. If we see that any squared distance is larger than a + 1/ poly(n), we know we are in Case
Therefore, we may assume in what follows that each column of U" V" has squared distance at most a + 1/ 1-olyol from the column span of W. In this case we can also determine in n 1+a(1) time, with constant probability, whether we are in Case 2. To do so, we apply leverage score approximation algorithm [Clarkson and Woodruff 2013) to find a small representative subset A' C A such that |A'l = (1) and if there exists a E A, E B such that a T b = 0, then a E A'. Specifically, suppose there is an E A and a b E B for which a 6 = O, and suppose U; extends by one coordinate and V; extends b by one coordinate. Then with probability at least 1/2, V2VV = 11, and so ((UVV)) 1. Let us condition on this event. Then we have the following claim in this case.
Claim 6. Let a E A and be E B be such that a = O, and the corresponding row Uj and column Vj satisfy V, 1. Also, suppose each column U" has squared distance at most a++ poly (n) from the column span of U" Then, in n (l(l) time, we can find a representative subset A' C A of size (l))) such that a E A' with high probability.
Proof. (of Claim) First note that by Remark there can be at most (1(1) coordinates in the j-th column of J(U V) which differ from their corresponding coordinate value in the j-th column of U"V", since each such difference corresponds to pair of points problem for which a² b = O. Hence, since f(U V) has squared distance at most a + 1/ poly(n) from the column span of W, then consider the following vectors by taking the difference of the j-th column of (U V) and the j-th column of U V", we would have vectors v denoting the difference vector and being the residual vector of projected on the column span of U", with the following properties:
A and b B in the OVP
v contains at most ( values, each of value equal to 2,
|||| 5 a + 1/ poly(n)
U + e is in the column span of U"
Let S C [n]] be the set of of indices i for which (u+e)?
||tt|| 2 1/) Since a + 1/ poly(n) is at most a
constant strictly less than 2, it follows that for each i in the support of since i2 2, we have that
(w+e)²
2 111 Since also ||v + ell < o()), the support of is included in S and |S] S o()),
l+tel
7

---

Recall the definition of leverage scores l from Section 2 and that the sum of the n leverage scores of B E Rx is exactly equal to the rank of B. which is at most t. Since v + e is in the column span of U", the i-th leverage score of U" satisfies ((())) 2 (u + /)l/ ell2 = (LL/no())
Consequently, each coordinate i in S satisfies (,())) = ()1/no())) It is known Clarkson and Woodruff 2013 how, given matrix B, in ( log n. (1)) time, one can compute list E 1...... C with C (4,) for all with probability 1 - 1/m 100 Consequently, using that rank(U") I rp, given U one can find a superset T containing S for which [T] = O(r"n(())) in n ++(()) time, where recall == + and 5 = (1) and p is constant. Note that our bound on [T] follows since we just need to keep the leverage scores that are (11/n(()) and the sum of all leverage scores is at most r" = no(1), recalling r = 8 + and = o())
Continuing the proof of Theorem 3.1 by our claim, when we are in Case 2. with at least constant probability, we find in time On a representative subset A' of size ()))))). Given our subset A', one can then compute all pairs of dot products between the points in A' and the points in B in O( ds) time, and since we may assume 2 d without loss of generality, we can compute all such pairs in n +o(l) time. Lastly, if no such a, b exist, when are in Case I, this algorithm cannot err. Therefore, we have an algorithm that can solve OVP, which then violates Assumption T using that (nd))=o(1) n (+0(1)
Our lower bound techniques also apply to number of important function f that do not have the form of |x]| for an integer p. We show that our lower bounds for LRA in fact holds for any function
g(x) = f[xx)), where f(x) is a function that admits a Taylor expansion with dominant term of
x², for odd p. around 0. Of particular interest in natural language processing (Liang et al.] 2020
Jiang et al.. 2021] is the function g(x) = log(1 + [æ]),, where satisfies g = f(zz)), where
[ = x + ( )) and thus we can appeal to our lower bound with p =1
Theorem 3.2. Suppose n E N, r = (log(n)), and g(x) f((z)), where f admits a Taylor
expansion f(x)= = cg + 0(xx[t+)) for odd p (log(n)). There is a positive integer k: n1))) such that for any approximation factor A 2 L, any, possibly randomized, algoritium which given
inputs U E Baxx and V E Bexxd outputs W E R = satisfying
ll - V.WW - 9(U V)II & A. Ilg(U v)) g(U VVII
with constant probability, requires (nd)1-o(1) time, under Assumption 1
Further, this holds even if
U = VT.
3.2 Polynomials of All Integers p
We now give the proof for p-degree polynomials for all integers p. which is weaker than Theorem 3.1 for odd integers, but this is the first such lower bound for even integers. The weaker bound cannot be substantially improved since in this setting fz) = x", so f(UV) is in fact a kernel and matrix vector multiplication can be performed in O(rr)) time. For intuition why our lower bound techniques do not extend, recall that our previous reduction forces the absolute value operation to essentially alter entries of (UV) but only at entries of the original with zero dot product. Therefore, we can write as sum of a low rank matrix (from tensor product) and a sparse matrix, whose sparse entries now represent the OVP pairs. However, when p even, the absolute value operation leaves the entries unchanged and does not induce the additional sparse matrix. Therefore, our lower bounds for this setting rely on a slightly different variant of OVP, specifically quadratic lower bounds for finding the maximum dot product (Assumption 4).
Intuitively, our new reduction is as follows: Let OPT be the maximum inner product and by assumption, consider this small set of large inner product pairs, which represents a small number of entries in that are large in magnitude so that when you apply a threshold
OPT/100, the resulting matrix is sparse. Since fx) = x² amplifies the magnitude differences, it follows that (UV) is much closer relatively to an approximately sparse, and therefore low rank, matrix. Therefore an approximate low rank approximation (LRA) algorithm can recover this sparse low rank matrix well enough so that the span of the approximate matrix can be used, via leverage score computations, to solve the APX-Max-IP problem.
8

---

Theorem 3.3. Let U E Re and V E R be given with r = O(log(n)). Suppose f(x) x¹ for an integer p 2 L. There is positive integer k n such that for any approximation factor A 2 1, any possibly randomized, algorithm which outputs W B satisfying [|f(U-V)WW* f(U VIII; = A .V)] - f(U V) with constant probability for any constant A > 1, requires min((nd)1-op) ,22()) time, under Assumption 4
The lower bound for general functions f(x in Theorem 3.2 with odd-degree dominating term in its Taylor expansion works because we can scale down the entries in our input matrix so that our relative-error approximation guarantee is still preserved but the function is largely approximated by only the leading polynomial term. By a combination of the same taylor expansion argument with a slight more general version of Theorem 3.3 that handles small additive error, we note that we can extend our results to when p is even.
We note that our lower bounds do not generalize to the case when U = V because as mentioned in our intuitive introduction, our reduction hinges on the sparse low-rank structure of (UV)'. However, when U = V. the sparsity structure is broken the diagonal of the matrix is now larger than the maximum inner product of two different vectors, and this destroys the low rank structure. In some sense, the diagonal of our matrix is forced to include the dot product of u, with itself and this shifts the entire matrix by a large multiple of the identity, crucially removing the sparse + low-rank structure that we exploited in lower bound argument before. Indeed,
remark below, the resulting positive semidefinite (PSD) structure allows us to derive a fast LRA algorithm and our lower bound no longer holds in this setting.
Remark 7. Recall that when U VT, there is an nr(k//)==> subquadratic time algorithm when k = (n((1)), for (1 + rr-rattivee error approximation, but reducing the U & V case to the U V case requires blowing up k to k rP, Therefore, our lower bounds imply that the U & V case is strictly more difficult in terms of runtime in certain settings.
4 LRA Algorithms from Matrix Vector Products
In this section, we can show upper bounds for low rank approximation of entrywise transformed products that are Onn ((())()))) when f(z represents a kernel matrix but U * V Specifically, we focus on polynomial functions of the form f(x) = I" and demonstrate that our lower bounds are tight for relative error LRA. Note that we have shown that low rank approximation guarantees for PSD matrices when U V cannot translate to the case when U A V Still, we demonstrate that relative error low rank approximation possible in nl(+(()) time for polynomial kernels with even degree, although there will be an exponential dependence on for relative error low rank approximation. Our relative error algorithms are relatively standard and rely on low rank projections by using matrix vector products to perform randomized sketching to reduce our row or column dimension to (ppoy((k/e)) [Woodruff] 2014).
Theorem 4.1. Let U E R2x and V E Rexd Supposee f(x)= x" for even integer p 2 1 and l < r.. For any approximation factor > O, there is an algorithm that outputs U² E BRn and V' E Raxd satisfying IIU" VV -(U 1 S (1 + +c) [[ VV)x - J(U - V) with constant probability with runtime O( d)rPk/e poly(r²/e))
Proof. Note that f(U.V,) (U,V))P. Thus, we may rewrite f(U . V) = U where each row of U E RRxr"" is the Khatri-Rao product of itself p times, and each column of V e RRrxxd the Khatri-Rao product of itself p times. Note that the rank of U"\" is at most dº.
Let S be random Gaussian sketching matrix with Okk/)) rows, so we know that these matrices satisfy the ( (Ve/k, 9/10, I)-JL property WWodrufff 2014). Furthermore, let be a random Gaussian matrix with (minn(k) (k/e3, r rp/))) columns, so we know that it a (1 (c)) lg subspace embedding of the row space of SU" Then, by Theorem 47 of Woodruff [0114), the following is true with constant probability
"S"""""""""""""""">""""""""""""" S(UVI S (1+ II V)) f(U V
Finally we bound the runtime computing this product. Note that we may compute ST and V in nr? (k/e + min(k/ē² ((2e))))) time. Then, note that the remaining products can be computed in

---

poly(kr2 /c) poly(r9/e (c) time. And lastly, the rank k approximation follows from solving for the low rank approximation on the restricted subspace SU"V R of rank and letting Z = VU" RU kk where U is the orthonormal basis such that UU the projection onto the row space of SAR.
We also provide an additive error guarantee on the low rank approximation guarantees that depends polynomially on p. as opposed to the tight exponential dependence in the relative error setting. This implies that additive error LRA is strictly easier and can be accomplished in subquadratic time when p (2))og(n)). Note that lower bounds for additive error are not emphasized in this paper as the upper bound is already quite competitive and simply outputting requires Q * poly(p))
Our guarantees follow from censor-based sketching techniques that are applied along each tensor dimension to remove the exponential dependence on p. We note that this results in an additive error term that is on the order of the p-norms of U, V. This term is related the absolute error guarantees of (Han et al. 2020), which also depends on the product of the p-norms of Euclidean norms of the rows of U.V.
Algorithm 2 TensorSketch LRA
Input: Matrices U E Raxx, V E R > O even integer
Output: rank approximation of f(UV), where f(x) = |x||
DSee Ahle et al. 2020 2: Compute U" = U and V TV", where U" E Raxr" V E R xd are U, V tensored themselves p times.
I: Let T be a TensorSketch with m OppE 2) rows
3: Let S, R. be random Gaussian matrices with Okk/e) rows and O(p/e columns respectively
4: Compute an orthonormal basis of row span of SU" (SU"V""RR)
>This implies
PP (SU"V"RR
5: Compute the rank & decomposition: U'V' [U"'V"RP):
>[M] the best rank
approximation to M
6: Output U", V'P (SU"'V"RR) SU"V""
Theorem 4.2. Let U E Raxr and V E Brxd, Suppose fxx) = xx for an even integer 2 1 and k < rP. For any approximation factor e >>> 0, there is an algorithm (Algorithm 2) that outputs U' E Bank and V' Rexd satisfying ||U". V f(U V)l S (1 c) [[[ VV) 1ol=-(UU f(U -V) + e² with constant probability with runtime O((n + d) polyppllypppp
k, 1/e)), where the additive term is 3p 2p
given by 2 = =E - [|U,]22)) [ ||U|| 2p.2
2p.2
4.1 Lower Bounds on Matrix Multiplication
Our low rank approximation algorithms use a matrix vector multiplication subroutine f(UV)z; for (poly(k//e) different vectors x; their main dimensionality reduction technique for subquadratic time upper bounds. This can be used to directly translate lower bounds for low rank approximation to implicit matrix vector multiplication for wide range of scalar functions that extends previous work beyond the exponential function. We emphasize that these -etrox product lower bounds also imply that our LRA lower bounds non-trivial and significantly strengthen those provided by previous works.
Theorem 4.3 (LRA reduces to Matrix Vector Products). Let U,V, be as Theorem 3.2 Then, any possibly randomized algorithm that for any vector E T outputs J(U (UV)) up to 1 polyy() entrywise error, with constant probability, requires (((n)))-oti) time, under Assumption / This
holds even when U = VT.
10

---

References
Thomas D Ahle, Michael Kapralov, Jakob BT Knudsen, Rasmus Pagh, Ameya Velingker, David Woodruff, and Amir Zandich. Oblivious sketching of high-degree polynomial kemels. In Proceed- ings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 141-160. SIAM, 2020.
Josh Alman and Zhao Song. Fast attention requires bounded entries. arXiv preprint ar arXiv:2302.13214 2023.
Josh Alman, Timothy Chu, Aaron Schild. and Zhao Song. Algorithms and hardness for linear algebra on geometric graphs. In 2020 IEEE 6l Annual Symposium on Foundations of Computer Science (FOCS), pages 541-552. IEEE, 2020.
Arturs Backurs, Piotr Indyk, and Ludwig Schmidt. On the fine-grained complexity of empirical risk minimization: Kernel methods and neural networks. Advances Neural Information Processing Systems, 2017.
Ainesh Bakshi, Nadiia Chepurko, and David P Woodruff. Robust and sample optimal algorithms for PSD low rank approximation. In 2020 IEEE st Annual Symposium on Foundations of Computer Science (FOCS), pages 506-516. IEEE, 2020.
Lijie Chen and Ryan Williams. An equivalence class for orthogonal vectors. In Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 21-40. SIAM, 2019.
Krzysztof Marcin Choromanski, Valerii Likhosherstov, David Dohan, Xingyou Song, Andreca Gane, Tamás Sarlós, Peter Hawkins, Jared Quincy Davis, Afroz Mohiuddin, Lukasz Kaiser, David Ben- jamin Belanger, Lucy J. Colwell, and Adrian Weller. Rethinking attention with Performers. In 9rh International Conference on Learning Representations, ICLR 2021, 2021.
Krzysztof Marcin Choromanski, Han Lin, Haoxian Chen, Arijit Sehanobish, Yuanzhe Ma, Deepali Jain, Jake Varley, Andy Zeng, Michael S. Ryoo, Valerii Likhosherstov, Dmitry Kalashnikov, Vikas Sindhwani, and Adrian Weller. Hybrid random features. In The Tenth International Conference on Learning Representations, ICLR 2022, 2022.
Kenneth L. Clarkson and David P. Woodruff. Low rank approximation and regression in input sparsity time. In Symposium on Theory of Computing Conference, STOC 2013, pages 81-90. ACM, 2013.
William Feller. An Introduction to Probability Theory and Its Applications, volume 1. Wiley, January 1968. ISBN 0471257087.
Insu Han, Haim Avron, and Jinwoo Shin. Polynomial tensor sketch for element wise function of low-rank matrix. In International Conference on Machine Learning, pages 2020.
PMLR,
Russell Impagliazzo and Ramamohan Paturi. On the complexity of k-SAT. Journal of Computer and System Sciences, 62(2):367-375, 2001.
Yifei Jiang, Yi Li, Yiming Sun, Jiaxin Wang, and David Woodruff. Single pass entrywise-transffrrme low rank approximation. In International Conference on Machine Learning (ICML), pages 4982- 4991. PMLR, 2021.
Praneeth Kacham, Vahab Mirrokni, and Peilin Zhong. Polysketchformer: Fast transformers via sketches for polynomial kernels. arXiv preprint arXiv:2310.01655, 2023.
Angelos Katharopoulos, Apoorv Vyas, Nikolaos Pappas, and François Fleuret. Transformers are RNNs: Fast autoregressive transformers with linear attention. In Proceedings of the 37th Inter- national Conference on Machine Learning. ICML 2020, 13-18 July 2020, Virtual Event, volume 119 of Proceedings of Machine Learning Research, pages 5156-5166 PMLR, 2020. URL html
http://proceedingg mml.innggs.m.......mmmmmm.. wttp://nw.........ssisii...
Feyza Duman Keles, Pruthuvi Mahesakya Wijewardena, and Chinmay Hegde. On the computational complexity of self-attention. International Conference on Algorithmic Learning Theory (ALT), 2023.
11

---

Omer Levy and Yoav Goldberg. Neural word embedding as implicit matrix factorization. In Advances in Neural Information Processing Systems, 2014.
Shaohua Li, Jun Zhu, and Chunyan Miao. generative word embedding model and its low rank positive semidefinite solution. In Proceedings of the 2015 Conference Empirical Methods in Natural Language Processing, pages 59991660, 2015.
Yingyu Liang, Zhao Song, Mengdi Wang, Lin Yang, and Xin Yang. Sketching transformed matri- with applications natural language processing. In International Conference Artificial Intelligence and Statistics, pages 467-481. PMLR, 2020.
Valerii Likhosherstov, Krzysztof M Choromanski, Kumar Avinava Dubey, Frederick Liu, Tamas Sarlos, and Adrian Weller. Chef'' random tables: Non-trigonometric random features. In Advances in Neural Information Processing Systems, pages 34559 14559-34573, 2022.
Daniel Lokshtanov, Dániel Marx, Saket Saurabh, et al. Lower bounds based on the exponential time hypothesis. Bulletin of EATCS, 3(105), 2013.
Cameron Musco and Christopher Musco. Recursive sampling for the Nystrom method. In Advances in Neural Information Processing Systems, 2017.
Cameron Musco and David P Woodruff. Sublinear time low-rank approximation of positive semidefi- nite matrices. In 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS), pages 672-683. IEEE, 2017.
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. Advances in Neural Information Processing Systems, 2007.
Yi Tay, Mostafa Dehghani, Dara Bahri, and Donald Metzler. Efficient transformers: A survey. ACM Computing Surveys, 55(6):1-28, 2022.
Yao-Hung Hubert Tsai, Shaojie Bai, Makoto Yamada, Louis-Philippe Morency, and Ruslan Salakhut- dinov. Transformer dissection: An unified understanding for transformer's attention via the lens of kernel. In Kentaro Inui, Jing Jiang, Vincent Ng, and Xiaojun Wan, editors, Proceedings of the 2019 Conference Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing. EMNLP-IJCNLL 2019, Hong Kong, China, November 3-7, 2019, pages 4343-4352. Association for Computational Linguistics, 2019. doi: 10.18653/v1/DI1 1444-1443, URL htt.hhhttp/ooott org/10. 10t.001115///1111001
Virginia Vassilevska Williams. Hardness of easy problems: Basing hardness on popular conjectures such as the strong exponential time hypothesis (invited talk). In 10th International Symposium on Parameterized and Exact Computation (IPEC 2015). Schloss Dagstuhl-LecbniizZZenttn fuer Informatik, 2015.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. In Advances in Neural Information Processing Systems 2017, pages 5998-6008, 2017.
Ryan Williams. A new algorithm for optimal 2-constraint satisfaction and its implications. Theor. Comput. Sei., 348 2005.
David P Woodruff. Sketching tool for numerical linear algebra. Foundations and Trends@ in Theoretical Computer Science, 10(1-2):1-157, 2014.
David P Woodruff and Peilin Zhong. Distributed low rank approximation of implicit functions of matrix. In 2016 IEEE 32nd International Conference on Data Engineering (ICDE), pages 847-858. IEEE, 2016.
12

---

A Missing Lower Bound Proofs
Proof of Theorem 3.2 Let U, VT 1,) be an instance of entrywise transformed low rank approximation with WLOG, by scaling our solution, we can let Cp = 1. Let B = poly(n) be sufficiently large, and let U U/VB and V = V/vB. For each i, [ni], for B sufficiently large and using that r = no(1),
s(Ü,V,) = F(U.VV)) f(UU,V,//B) = UNV
U.V Bo+1
0
BP
IUNV
rP
M((,VV)) (
rp
B
Bp+1
BB+1
where h(x) | Hence, if Ü", V are such that
I10 V' g(ū V)l S Allg(Ü-V)) 100k-g(Ü.. -(UVV
then by the triangle inequality,
1 -' -Nū Vllrr
- lgg(Ü. .V) h(U 1IIl= + VAllo(t. V) glle 1l
S 0 Ba+1
nr
+ VAllo( V) - I(ū Vill..
Setting ='= VBB. Ú and V VB. V and scaling both sides by B, we have
nr
110" - V)llE S o
+ BVAIIb( V)Da - g(ū V)lle
B
nr
0
+ BVΣ V)Dx 1900 Ville
By
nrP
0
+ BVDIINÜ V)Da - -hO Ville
B
+
BV BVAAH(U V) - g(ū VIILF
Anr
0
BBDIII( VDDa 1 h(O VIIl=
Br
VAnre
o
VAIIM(T V) A(U M(UVV))l.
B"
where in the second inequality we used that [h( V)]x has rank whereas (g(ū V)]x is the best rank-k approximation to g(ū V) in Frobenius norm. Note that the third line is triangle inequality and the fourth line follows entrywise approximation bounds.
Now by AM-GM, using that a < + e implies that a² S 26 + 2c for a, b,c 2 0, we have
2p
|I' V h(U VI O
24IIb(U V) h(U -V
B
Thus, U² V a rank-k approximation h(U V) with additive error poly(n). by setting B to be large enough, and relative error 24. The total time to find U' and V is the same as the time to solve entrywise transformed low rank approximation with respect to the function on inputs U and V, and thus by Theorem 3.1 is at least (nd))-o(1) Note that applying the case U V Theorem 3.1 establishes this theorem when Ü = VT
Proof of Theorem 3.3 Suppose A and B are the input sets to the Apx-Max-IP, n.d Problem of As- sumption 4 with parameter = We let the rows of U be the points in A and we let the columns of V be the points in B. As in the proof of Theorem 3.1 we have fUT
so J(U V)= = U" V". By Remark 5 there are at most (TT) entries of f(U V) that are at least (m/100)p.
13

---

Hence, if we set & to be an appropriate value in n²(1,, then
|| V) -UU V)ll N (m/)mm/100)3)
since one possible rank-k approximation is just to zero out the at most n2 entries of f(U V) that are at least (m/ (0/1 By the correctness guarantee, with constant probability,
ll VOWWT f(U V)I S A (nd)(m/100)
O(nd)(m/100) 100)
(A.1)
We may assume WLOG that p 2 C' log n. for the moment, for a sufficiently large constant C' > 0 since a trivial lower bound for outputting an LRA is R(()), so the inclusion of the 20(p) term in the lower bound allows us to make this assumption. Therefore, we conclude that J(U V) is in the column span of W up to 1/poly(n) error. Consequently, consider the maximum entry m, which we can assume is at least 1, as otherwise all points in A would have disjoint support from those in B but this can be verified in ((nr) = n +11) time, contradicting the n (+8(11 lower bound for d = n ((1) in Assumption 4
It follows that if the maximum m occurs in the (i.j)-th entry of f(U VV U"V", then there exists a vector in the column span of U? for which the i-th entry is [m - (alppoly)) (n)) (1+ (o)y(n)) and all other entries are in the range i-m" /poly(pp))) m² /poly(n)] since (1 [1/] = 1/poly(n). Consequently, since the column span of W close to the column span of U up to 1/poly(n) error, we see that the column span of W also contains these almost-sparse vectors that have n2 large entries. As argued before similarly in Theorem 3.1 the rows of these large dot products forces the row leverage scores of W to be ((1). Since W has rank k = n2(11, in +11) time we can find the set S of (1) row leverage scores of W that are ((1). We can then explicitly compute all dot products between pairs of vectors in A n S and B in +o(() total time, at which point we can output the maximum dot product, which includes m. Therefore by Assumption the total time to find W is at least (nd)²-0 lo(())
Proof of Theorem 4.3 This proof follows by reducing low rank approximation to matrix multipli- cation by the same algorithm as in Theorem 4.1 Specifically, let us set k = n1) and = d and note that our low rank approximation algorithm's 's runtime is dominated by computing S (UV) and f(UV) R, where S, R are matrices with poly(k/e) rows and columns, respectively. In fact, those operations are the only terms that incur polynomial dependence on n.
Therefore, suppose there exists such a matrix multiplication algorithm that takes time O(n2-c) for some c. Then this would directly imply that computing both S (UV) R takes 0 n 2-2200 29) time as [S](UV)] S(V U ISST. USS" By our guarantees in Theorem 47 of Woodruff [2014] this
implies a constant relative error LRA in time O(n which contradicts Theorem 3.2
B Missing Upper Bound Proofs
Again, f(UVV)) (U:V))P so we may rewrite J(U V) U"V", where each row of U E B"xr" is the Khatri-Rao product of itself p times, and each column of V e e is the Khatri-Rao product of itself, times.
Proof of Theorem 4.2
Using Theorem 1 of Ahle et al. [2020]. we see that with probability 0.9, our approximate matrix product guarantees hold for U" V" such that if an m x r2 TensorSRHT matrix with m = ((p/e), then
| IV" -U e|
Therefore, use the approximate low rank sketches again to approximately solve:
min UV I TV" - -V
Specifically, let S be random Gaussian sketching matrix with O(k/c) rows. We know that these matrices satisfy the (Ve/k, 9/10 9/10, I)- property Woodruff 2014). Furthermore, let R be a random
14

---

Gaussian matrix with Om columns, so we know that is (1 + O(e)) & subspace embedding of the row space of SU"IT,, which has rank at most m. Let U" = U"IT" and let ITV V" Then, by Theorem 47 of Woodruff [014). the following is true with constant probability
(I"""""""""""""" U"" (1+ min |IV - U""V"ll
Finally we bound the runtime of computing this product. Note that we may compute U" and V in time (mp((m + r)) time. Then, computing SU"", v" R can be done in mm (k/e + m/<²) time. Lastly, the remaining products can be computed in poly(kp/e) and the final rank decomposition can be computed and can find U", V" such that
I - -"' (1+ min JUV - U"'V"">
U.
Now, let U* V* V* be the optimal rank decomposition of the original problem of f(UV), then by the guarantees of approximate matrix product and the triangle inequality, |U-V UUV" - U IVV "
=||U|||||||"||
Therefore, we conclude that
IIU'V' II "Iv s (1+ (1+e)| uU'''
-UVV
We end by rewriting | TT(U""") (UTUU E, 11 1 and similarly for || and then applying AM-GM and noting that (1+ = =1+0)()
15


Reviews:
[{"review_id": "71kRHVkZlQ", "reviewer": "Reviewer_9ux2", "summary": "For two matrices U,V one is often interested in computing f(UV) for some entrywise function f. Several papers in the past decade considered the problem of approximating f(UV) by a product AB with small inner-dimension (low-rank approximation). While a number of past papers presented clever positive results, this paper is about demonstrating various negative results. I am very far from the field of complexity and cannot evaluate soundedness of the claims. However, I should say that the negative results are probably best presented in a different venue (SODA etc), since this conference is mostly interested in doing things rather than in enumerating all things that cannot be done.", "questions": "n/a", "limitations": "n/a", "rating": 6, "confidence": 1, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "Key problem of contemp interest.", "weaknesses": "Complexity results are likely irrelevant to most attendees.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive feedback. We would like to address your main comment concerning “irrelevant complexity results to most attendees” and “[NeurIPS] is mostly interested in doing things rather than in enumerating all things that cannot be done”: \n\n* First, while our submission mainly focuses on impossibility results, we note that we also provide relative and additive error low rank approximation algorithms that match our lower bounds of $n * 2^p$ and $n * \\text{poly}(p)$, respectively.\n \n* Furthermore, all of our lower bounds are based on algorithmic reductions. We focus on the lower bounds as they present more novelty and difficulty than our upper bounds and in fact, stand out in a field that is dominated by positive results. \n\n* Given the large focus on Transformers currently, our paper’s resulting conclusions for the feasibility of efficient attention is very relevant and timely for NeurIPS.\n\n* More generally in the field of machine learning, we believe there have been multiple influential complexity (i.e., impossibility) results published at NeurIPS, to name a few out of many:\n  * Coarse Sample Complexity Bounds for Active Learning (Dasgupta, 2005)\n  * Sample Complexity of Testing the Manifold Hypothesis (Narayanan and Mitter, 2010)\n  * Sample Complexity of Episodic Fixed-Horizon Reinforcement Learning (Dann and Brunskill, 2015)\n  * Lower Bounds on the Robustness to Adversarial Perturbations (Peck et al, 2017)\n  * On the Complexity of Learning Neural Networks (Song, 2017)\n\n* Even more generally, we believe that lower bounds are important as they show the limitations of what is algorithmically possible, and so it makes sense for them to be presented along with upper bounds at the same venue for visibility, so that researchers know when to stop searching for better algorithms, without, e.g., making additional assumptions or changing the problem. We believe this viewpoint is shared by many at NeurIPS, and that the call for papers certainly includes such results. \n"}, {"review_id": "6nqVcGkEVY", "reviewer": "Reviewer_GCRp", "summary": "This paper studies the problem of low rank approximation of the form f(U.V) where the matrices U,V $\\in \\mathbb{R}^{n \\times r}$ and $r=O(\\log n)$ and $f(.)$ is applied entrywise. The main contributions of the paper are the following:\n\n1) The paper presents a $\\Omega(n^{2-o(1)})$ time lower bound for getting relative error rank $k$ approximation to $f(U.V)$ for $k=n^{o(1)}$ and $f(x)$ admits a Taylor series expansion of the form $f(x)=x^p+O(|x|^{p+1})$ when $p$ is odd and $p=O(\\log n)$. This lower bound holds for any $U, V$. Another $\\min(n^{2-o(1)},2^{\\Omega(p)})$ time lower bound, which only holds when $U \\neq V^T$, is presented for the case when $f(x)=x^p$ for any general integer $p \\geq 1$. The first lower bound proof for odd $p$ follows by a reduction of the Orthogonal Vectors Problem to the Low-Rank Approximation (LRA) problem of $f(U.V)$. The proof then proceeds by lower bounding the leverage score of rows corresponding to non-zero entries of a flat sparse vector in the column span of the low-rank approximation of $f(U.V)$. The second lower bound proof for a general $p$ follows by a reduction of the Approximate Maximum Inner Product problem respectively to the LRA problem.\n\n2) The paper presents a $O(n . poly(2^p,k,\\frac{1}{\\epsilon}))$ time algorithm for getting a $1+\\epsilon$ relative error low-rank approximation to $f(U.V)$ mathcing the lower bound. The algorithm uses random Gaussian sketching of the Khatri-Rao products of U and V followed by a low-rank approximation of the smaller sketches. The paper also presents a $O(n, poly(k,p, \\frac{1}{\\epsilon}))$ time algorithm for getting a low-ran approximation with additive $\\epsilon . L$ guarantee where $L$ depends on the $p$-norms of $U,V$. Thi algorithm proceeds by doing a tensor sketch of the Kahtri-Rao products of $U$ and $V$.\n\n3) Finally, the paper presents a lower bound of $\\Omega(n^{2-o(1)})$ on matrix vector products of the form $f(UV)z$ for a vector $z \\in \\mathbb{R}^n$ using the fact that the algorithms use a matrix vector product subroutine.", "questions": "Other than the points raised above, there doesn't seem to be any discussion on the lower bound for additive error so it is unclear how good the upper bound presented in Theorem 4.2 is. It might be useful to add a short discussion there.\n", "limitations": "I don't think there are any other limitations other than the ones mentioned above.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1) Previous subquadratic time algorithms for relative error low-rank approximation for $f(U.V)$ assume that $U=V^T$ and $f(x)$ is a PSD kernel function. The quadratic lower bounds for relative error low-rank approximation of $f(U.V)$ for polynomial $f(.)$ show that both these assumptions are necessary for a subquadratic algorithm. This is an interesting result of broader interest to the sketching community. The lower bound proof techniques via reductions from OVP and MIP problems and bounding the leverage score of flat sparse vectors are also novel.\n\n2) The paper is generally easy to follow for the most part.", "weaknesses": "1) Apart from the lower bounds of Theorem 3.1 and 3.3, the other results seem to be straightforward extensions of these results or, for the case of the upper bounds, applications of existing algorithms and results in the literature. \n\n2) Though the results are interesting, the intuition behind some of the theoretical results presented in the paper are unclear. \n\n     a) Specifically, it might be useful to explain more in detail why is there a separation in the lower bound between the $U=V^T$ and $U \\neq V^T$ cases for general $p$.\n\n     b) It might be good to explain why the lower bound for general $p$ extend to more general functions with Taylor series of the form $f(x)=x^p+O()$ like in the case for odd $p$. Where is the bottleneck for general $p$?\n\n2) Though this is primarily a theoretical paper, in the context of Neurips, some experimental evaluation of the algorithms for LRA might be useful for understanding how well the algorithms perform well in practice. \n\nThere are also some typos and confusing notations in the paper:\n1) In the line 260, $Z$ has dimension $n \\times s$. The $s$ here is the rank of $U''$ and is different from the dimension of the  vectors $a_i$ and $b_i$.\n2) In the statement of Theorem 4.1, on the RHS, it seems the $U,V$ in argument $\\min_{U,V}$ should be changed to different notation as this is different from the input to the problem.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful review of our paper - the following are our responses to the main concerns raised:\n\nAddressing Weaknesses:\n\n1. We emphasize that we presented our results in a way that highlights the novelty and difficulty of our lower bounds, which could be of general interest even in algorithm design. Our upper bounds, though novel to our knowledge, are indeed relatively mathematically straightforward compared to our lower bounds. However, their main purpose is to supplement the lower bounds, not necessarily in novelty or difficulty, but in showing that our derived lower bounds are in fact optimal. Therefore, including such upper bounds is important.\n\n2. Yes, we will add more intuition behind the theoretical results, which we mention below:\n\n  * There is a separation of algorithmic time complexity for relative low rank approximation (LRA) between (1) the $U = V^T$ setting and (2) the $U \\neq V^T$ setting. For the sake of illustration, let’s consider the case when $f(x) = x^{2p}$ for $p = \\Omega(\\log(n))$ sufficiently large and $k = O(1)$. In the setting 1, we have an $n * \\text{poly}(\\log n)$ time algorithm for constant factor LRA using the cited paper by Bakshi et al.,  2020. However, in setting 2, we provide a lower bound of nearly $n^2$. Therefore, we conclude that there is an inherent complexity separation in the two settings.\n\n  * The extension of the lower bound for general $p$ occurs because the Taylor expansion of the function is dominated by the $x^p$ term around $0$. Therefore, we can scale down the entries in our input matrix so that our relative-error approximation guarantee is still preserved but the function is dominated by the leading polynomial term (please see Theorem 3.2 for more details). In fact, this technique also generalizes to the general $p$ setting; however, it was not included for brevity and we will include this in the final version.\n\n3. Experimental evaluation: while there are many papers focusing on empirical evaluations, our focus is on our lower bounds as they require novel arguments involving leverage scores, and stand out in a field that is dominated by positive results. \n\nAddressing the additional questions:\n\n* “No discussion on lower bounds for additive error” - Lower bounds for additive error LRA are not emphasized because the upper bound is already quite small and is based on the optimality of TensorSketch in Ahle et al., 2020. Specifically note that when $k = \\text{poly}(p)$, simply outputting $U, V$ already takes $\\Omega(n * \\text{poly}(p))$ times and our upper bound on additive LRA is $O(n * \\text{poly}(p))$. Our upper bounds for additive error were used to highlight a separation between additive and multiplicative error. We will add this discussion in our final version.\n\nWe will use your detailed feedback in improving our paper and are happy to address any other questions you have. If our responses sufficiently address the weaknesses and questions that you raised, we kindly and respectfully ask if our rating could be improved. \n"}, {"review_id": "RP7WroHIwt", "reviewer": "Reviewer_ExZq", "summary": "The authors discuss a setting where one observes a scalar transformation of an r-rank matrix, f(UV) and wants to find the best k-rank approximation of it.\n\nIt is known that under various assumptions, if f(UU^t) is PSD then the problem is solvable in nearly linear time. It is also known that for various kernels f (eg Gaussian kernel), one needs quadratic time.\n\nThe authors discuss the cases where f are simple monomials, f(x)=x^p. They show a min(2^p,n^2) time lower bound based on SETH.  They also show a matching upper bound, among other results.\n\n\n", "questions": "See above.", "limitations": "See above.", "rating": 6, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "The paper seems novel and interesting. The connection with SETH is non-trivial and intruiguing.", "weaknesses": "I am failing to see the importance of choosing the x^p activation function (besides the log(1+|x|) example). Could the authors elaborate more on it?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive comments. \n\nWe address your comment regarding the “importance of choosing the $x^p$ activation function”: Note that we extend our base case results over $x^p$ to a broader class of functions (via a Taylor expansion) in Theorem 3.2. Specifically, our results also hold for any function $g(x) = f(|x|)$ where $f(x)$ is a function that admits a Taylor expansion with a dominant term of $x^p$ for odd $p$, around $0$. Examples include most activation functions, as they have a dominant linear term, such as $\\text{sigmoid}(x)$ or unnormalized softmax $\\exp(x)$, which behave like $1 + x$ around $0$.\n"}, {"review_id": "ir1xVkiQHM", "reviewer": "Reviewer_26kb", "summary": "This paper studies the problem of computing a low-rank approximation (LRA) for $f(UV)$, where $U\\in \\mathbb{R}^{n\\times r}$, $V\\in \\mathbb{R}^{r\\times n}$, and $f:\\mathbb{R}\\rightarrow \\mathbb{R}$ applies to each entry of $UV$. This problem has very important applications in deep learning, in particular for Transformers and natural language processing. The main results of this paper consist of two parts:\n\nFor lower bounds, they show that we cannot get subquadratic relative error LRA generally when either 1) $U\\ne V^\\top$ even for PSD kernel functions or 2) for $f$ that are approximately polynomials of $|x|$ even for a constant degree. \n\nFor upper bounds, they show that for $f(x) = x^p$, there is an $O(n\\cdot poly(r^p, k, 1/\\epsilon))$-time algorithm for relative error LRA and $O(n \\cdot poly(p, k, 1/\\epsilon))$-time for additive error LRA. \n\nFurthermore, their techniques can be generalized to obtain lower bounds for the implicit matrix-vector multiplications of the form $f(UV)z$.\n", "questions": "1. Line 157: define $B_i$.\n\n2. Do the lower bounds for OV (Assumption 1) and Max-IP (Assumption 4) hold for randomized algorithms? They may be based on (randomized) SETH.\n\n3. Line 262: $s\\leq r^p$.$s$ has already been used to represent the number of vectors in OV. I guess the $s$ here is different?\n\n4. Line 265: the derivation of the below equation via the Pythagorean theorem is unclear to me. And how does it imply that when we are in case 2, there is a column of $U’V’$ that is $\\alpha+1/poly(n)$ far from the column span of $U’’$?\n\n5. Line 287: how do you construct the vectors $v$ and $e$?\n\n6. Line 306: where does the $n^{o(1)}dr^ps$ come from? Shouldn’t it be $n^{o(1)}ds$?\n", "limitations": "n/a", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This paper significantly advances our understanding of LRA with very solid results. Prior to this work, most of the research focuses on the symmetric setting where $U=V^\\top$. This paper considers the more general, asymmetric setting. For relative error LRA and $f(x)=x^p$, their lower bound almost matches their upper bound in terms of the dependence on $p$. On the other hand, they also make progress in the symmetric setting. Their lower bound in the symmetric setting applies to a larger family of the function $f$ compared to previous results. Technically, they provide clean and novel reductions from the orthogonal vector problem (OV) to LRA via the fast computation of leverage scores. The reductions in this paper will be useful in proving new lower bounds for other numerical linear algebra problems. Additionally, this paper provides a clear review of many related works about LRA. \n\n", "weaknesses": "The algorithms proposed in this paper is quite straightforward given prior results  [Woodruff, 2014] and [Ahle et al., 2020]. Also, the lower bound for LRA with an even $p$ is weak. And the lower bounds do not depend on $r$ and $k$, the dimensions of the input matrices and output matrices, respectively. Moreover, there are some typos in the proofs that affects reading this paper. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive review and careful reading of our paper.\n\nWe address your questions below:\n\n1. $B_i$ is the $i$-th row of matrix $B$.\n\n2. Yes, the lower bounds do hold for randomized algorithms (and OVP is stated with constant probability). In general, by Yao’s minimax principle, deterministic lower bounds translate into lower bounds for randomized algorithms.\n\n3. Generally, $n$ is used to represent the number of vectors in $OV$, and $s$ is defined in Theorem 3.1/Algorithm 1 and is used to denote the dimension of the input vectors. We will clear up any confusing notation.\n\n4. Yes, we will clarify the notation and derivations. Generally, we are claiming that we can compute all squared distances of $U’V’$ to the column span of $U’'$ in $O(n^{1+o(1)})$ time. This is so that we can detect whether we are in case 1, where $U’V’$ has squared distance at most alpha to the column span of $U’’$, or we are in case 2 (the opposite case). To compute this, we can first compute orthogonal bases $Z$, $W$, respectively for $U’’$ and $U’$, and then we can compute the squared distance using the Pythagorean theorem, since the squared distance is equal to (total squared distance) - (projected squared distance).\n\n5. Let $v$ = distance of $j$-th column of $U’’V’$’ to $j$-th column of $U’V’$, and $e$ = residual vector of $v$ projected on the column span of U’’, such that $v + e$ is in the span.\n\n6. Yes, thanks that is indeed a typo from a previous version and it is fixed.\n"}]

Author rebuttal:

Meta-review: The reviewers agree unanimously that the paper is solid and worthy of acceptance. Please incorporate reviewer feedback into the final camera-ready version.
Decision: Accept (poster)