Faster Linear Algebra for Distance Matrices

Our dataset will be the $n$ points $X=\{x_{1},\ldots,x_{n}\}\subset\mathbb{R}^{d}$. For points in $X$, we denote $x_{i}(j)$ to be the $j$th coordinate of point $x_{i}$ for clarity. For all other vectors $v$, $v_{i}$ denotes the $i$th coordinate. We are interested in matrices of the form $A_{i,j}=f(x_{i},x_{j})$ for $f:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ which measures the similarity between any pair of points. $f$ might not necessarily be a distance function but we use the terminology “distance function” for ease of notation. We will always explicitly state the function $f$ as needed. $w\approx 2.37$ denotes the matrix multiplication constant, i.e., the exponent of $n$ in the time required to compute the product of two $n\times n$ matrix [AW21].

Our work can be understood as being part of a long line of classical works on the matrix free or implicit model as well as the active recent line of works on the matrix-vector query model. Many widely used linear algebraic algorithms such as the power method, the Lanczos algorithm [Lan50], conjugate gradient descent [S⁺94], and Wiedemann’s coordinate recurrence algorithm [Wie86], to name a few, all fall into this paradigm. Recent works such as [MM15, BHSW20, BCW22] have succeeded in precisely nailing down the query complexity of these classical algorithms in addition to various other algorithmic tasks such as low-rank approximation [BCW22], trace estimation [MMMW21], and other linear-algebraic functions [SWYZ21b, RWZ20]. There is also a rich literature on query based algorithms in other contexts with the goal of minimizing the number of queries used. Examples include graph queries [Gol17], distribution queries [Can20], and constraint based queries [ES20] in property testing, inner product queries in compressed sensing [EK12], and quantum queries [LSZ21, CHL21].

Most prior works on query based models assume black-box access to matrix-vector queries. While this is a natural model which allows for the design non-trivial algorithms and lower bounds, it is not always clear how such queries can be initialized. In contrast, the focus of our work is not on obtaining query complexity bounds, but rather complementing prior works by creating an efficient matrix-vector query for a natural class of matrices.

Most work on subquadratic algorithms for distance matrices have focused on the problem of computing a low-rank approximation. [BW18, IVWW19] both obtain an additive error low-rank approximation applicable for all distance matrices. These works only assume access to the entries of the distance matrix whereas we assume we also have access to the underlying dataset. [BCW20] study the problem of computing the low-rank approximation of PSD matrices with also sample access to the entries of the matrix. Their results extend to low-rank approximation for the $\ell_{1}$ and $\ell_{2}^{2}$ distance matrices in addition to other more specialized metrics such as spherical metrics. Table 2 lists the runtime comparisons between their results and ours.

Practically, the algorithm of [IVWW19] is the easiest to implement and has outstanding empirical performance. We note that we can easily simulate their algorithm with no overall asymptotic runtime overhead using $O(\log n)$ vector queries. Indeed, their algorithm proceeds by sampling rows of the matrix according to their $\ell_{2}^{2}$ value and then post-processing these rows. The sampling probabilities only need to be accurate up to a factor of two. We can acquire these sampling probabilities by performing $O(\log n)$ matrix-vector queries which sketches the rows onto dimension $O(\log n)$ and preserves all row-norms up to a factor of two with high probability due to the Johnson-Lindenstrauss lemma [JL84]. This procedure only incurs an additional runtime of $O(T\log n)$ where $T$ is the time required to perform a matrix-vector query.

The paper [ILLP04] shows that the exact $L_{1}$ distance matrix can be created in time $O(n^{(w+3)/2})\approx n^{2.69}$ in the case of $d=n$, which is asymptotically faster than the naive bound of $O(n^{2}d)=O(n^{3})$. In contrast, we focus on creating an (entry-wise) approximate distance matrices for all values of $d$.

We also compare to the paper of [ACSS20]. In summary, their main upper bounds are approximation algorithms while we mainly focus on exact algorithms. Concretely, they study matrix vector products for matrices of the form $A_{i,j}=f(\|x_{i}-x_{j}\|_{2}^{2})$ for some function $f:\mathbb{R}\rightarrow\mathbb{R}$. They present results on approximating the matrix vector product of $A$ where the approximation error is additive. They also consider a wide range of $f$, including polynomials and other kernels, but the input to is always the $\ell_{2}$ distance squared. In contrast, we also present exact algorithms, i.e., with no approximation errors. For example one of our main upper bounds is an exact algorithm when $A_{i,j}=\|x_{i}-x_{j}\|_{1}$ (see Table 1 for the full list). Since it is possible to approximately embed the $\ell_{1}$ distance into $\ell_{2}^{2}$, their methods could be used to derive approximate algorithms for $\ell_{1}$, but not the exact ones. Furthermore, we also study a wide variety of other distance functions such as $\ell_{\infty}$ and $\ell_{p}^{p}$ (and others listed in Table 1) which are not studied in Alman et al. In terms of technique, the main upper bound technique of Alman et al. is to expand $f(\|x_{i}-x_{j}\|_{2}^{2})$ and approximate the resulting quantity via a polynomial. This is related to our upper bound results for $\ell_{p}^{p}$ for even $p$ where we also use polynomials. However, our results are exact, while theirs are approximate. Our $\ell_{1}$ upper bound technique is orthogonal to the polynomial approximation techniques used in Alman et al. We also employ polynomial techniques to give upper bounds for the approximate $\ell_{\infty}$ distance function which is not studied in Alman et al. Lastly, Alman et al. also focus on the Laplacian matrix of the weighted graph represented by the distance matrix, such as spectral sparsification and Laplacian system solving. In contrast, we study different problems including low-rank approximations, eigenvalue estimation, and the task of initializing an approximate distance matrix. We do not consider the distance matrix as a graph or consider the associated Laplacian matrix.

It is also easy to verify the “folklore” fact that for a gram matrix $AA^{T}$, we can compute $AA^{T}v$ in time $O(nd)$ if $A\in\mathbb{R}^{n\times d}$ by computing $A^{T}v$ first and then $A(A^{T}v)$. Our upper bound for the $\ell_{2}^{2}$ function can be reduced to this folklore fact by noting that $\|x-y\|_{2}^{2}=\|x\|_{2}^{2}+\|y\|_{2}^{2}-2\langle x,y\rangle$. Thus the $\ell_{2}^{2}$ matrix can be decomposed into two rank one components due to the terms $\|x\|_{2}^{2}$ and $\|y\|_{2}^{2}$, and a gram matrix due to the term $\langle x,y\rangle$. This decomposition of the $\ell_{2}^{2}$ matrix is well-known (see Section $2$ in [DPRV15]). Hence, a matrix-vector query for the $\ell_{2}^{2}$ matrix easily reduces to the gram matrix case. Nevertheless, we explicitly state the $\ell_{2}^{2}$ upper bound for completeness since we also consider all $\ell_{p}^{p}$ functions for any integer $p\geq 1$.

There have also been works on faster algorithms for approximating a kernel matrix $K$ defined as the $n\times n$ matrix with entries $K_{i,j}=k(x_{i},x_{j})$ for a kernel function $k$. Specifically for the polynomial kernel $k(x_{i},x_{j})=\langle x_{i},x_{j}\rangle^{p}$, recent works such as [ANW14, AKK⁺20, WZ20, SWYZ21a] have shown how to find a sketch $K^{\prime}$ of $K$ which approximately satisfies $\|K^{\prime}z\|_{2}\approx\|Kz\|_{2}$ for all $z$. In contrast, we can exactly simulate the matrix-vector product $Kz$. Our runtime is $O(nd^{p})$ which has a linear dependence on $n$ but an exponential dependence on $p$ while the aforementioned works have at least a quadratic dependence on $n$ but a polynomial dependence on $p$. Thus our results are mostly applicable to the setting where our dataset is large, i.e. $n\gg d$ and $p$ is a small constant. For example, $p=2$ is a common choice in practice [CHC⁺10]. Algorithms with polynomial dependence in $d$ and $p$ but quadratic dependence in $n$ are suited for smaller datasets which have very large $d$ and large $p$. Note that a large $p$ might arise if approximates a non-polynomial kernel using a polynomial kernel via a taylor expansion. We refer to the references within [ANW14, AKK⁺20, WZ20, SWYZ21a] for additional related work. There is also work on kernel density estimation (KDE) data structures which upon query $y$, allow for estimation of the sum $\sum_{x\in X}k(x,y)$ in time sublinear in $|X|$ after some preprocessing on the dataset $X$. For widely used kernels such as the Gaussian and Laplacian kernels, KDE data structures were used in [BIMW21] to create a matrix-vector query algorithm for kernel matrices in time subquadratic in $|X|$ for input vectors which are entry wise non-negative. We refer the reader to [CS17, BCIS18, SRB⁺19, BIW19, CKNS20] for prior works on KDE data structures.

We first consider the case that all points $x\in X$ are from $\{0,1\}^{d}$. We first claim the existence of a polynomial $T$ with the following properties. Indeed, the standard Chebyshev polynomials satisfy the following lemma, see e.g., see Chapter 2 in [SV⁺14].

Now note that $\|x-y\|_{\infty}$ can only take on two values, $0$ or $1$. Furthermore, $\|x-y\|_{\infty}=0$ if and only if $\|x-y\|_{2}^{2}=0$ and $\|x-y\|_{\infty}=1$ if and only if $\|x-y\|_{2}^{2}\geq 1$. Therefore, $\|x-y\|_{\infty}=0$ if and only if $T(\|x-y\|_{2}^{2}/d)=0$ and $\|x-y\|_{\infty}=1$ if and only if $|T(\|x-y\|_{2}^{2}/d)-1|\leq\varepsilon$. Thus, we have that

for all entries $A_{i,j}$ of $A$. Note that $T(\|x-y\|_{2}^{2}/d)$ is a polynomial with $O((2d)^{t})$ monomials in the variables $x(1),\ldots,x(d)$. Consider the matrix $B$ satisfying $B_{i,j}=T(\|x_{i}-x_{j}\|_{2}^{2}/d)$. Using the same ideas as our upper bound results for $f(x,y)=\langle x,y\rangle^{p}$, it is straightforward to calculate the matrix vector product $By$ (see Section A.2). To summarize, for each $k\in[n]$, we write the $k$th coordinate of $By$ as a polynomial in the $d$ coordinates of $x_{k}$. This polynomial has $O((2d)^{t})$ monomials and can be constructed in $O(n(2d)^{t})$ time. Once constructed, we can evaluate the polynomial at $x_{1},\ldots,x_{n}$ to obtain all the $n$ coordinates of $By$. Each evaluation requires $O((2d)^{t})$ resulting in an overall time bound of $O(n(2d)^{t})$.

We now consider the case that all points $x\in X$ are from $\{0,\ldots,M\}^{d}$. Our argument will be a generalization of the previous section. At a high level, our goal is to detect which of the $M+1$ possible values in $\{0,\ldots,M\}$ is equal to the $\ell_{\infty}$ norm. To do so, we appeal to the prior section and design estimators which approximate the indicator function $``\|x-y\|_{\infty}\geq i"$. By summing up these indicators, we can approximate $\|x-y\|_{\infty}$.

Our estimators will again be designed using the Chebyshev polynomials. To motivate them, suppose that we want to detect if $\|x-y\|_{\infty}\geq i$ or if $\|x-y\|_{\infty}<i$. In the first case, some entry in $x-y$ will have absolute value value at least $i$ where as in the other case, all entries of $x-y$ will be bounded by $i-1$ in absolute value. Thus if we can boost this ‘signal’, we can apply a polynomial which performs thresholding to distinguish the two cases. This motivates considering the functions of $\|x-y\|_{k}^{k}$ for a larger power $k$. In particular, in the case that $\|x-y\|_{\infty}\geq i$, we have $\|x-y\|_{k}^{k}\geq i^{k}$ and otherwise, $\|x-y\|_{k}^{k}\leq di^{k-1}$. By setting $k\approx\log(d)$, the first value is much larger than the latter, which we can detect using the ‘threshold’ polynomials of the previous section.

We now formalize our intuition. It is known that appropriately scaled Chebyshev polynomials satisfy the following guarantees, see e.g., see Chapter 2 in [SV⁺14].

Given $x,y\in\mathbb{R}^{d}$, our estimator will first try to detect if $\|x-y\|_{\infty}\geq i$. Let $T_{1}$ be a polynomial from Lemma 3.4 with $t=O(M^{k})$ for $k=O(M\log(Md))$ and assuming $k$ is even. Let $T_{2}$ be a polynomial from Lemma 3.4 with $t=O(\sqrt{d}\log(M/\varepsilon))$. Our estimator will be

If coordinate $j$ is such that $|x(j)-y(j)|\geq i$, then

and so $T_{1}$ will evaluate to a value very close to $1$. Otherwise, we know that

by our choice of $k$, which means that $T_{1}$ will evaluate to a value close to $0$. Formally,

will be at least $1/d$ if there is a $j\in[d]$ with $|x(j)-y(j)|\geq i$ and otherwise, will be at most $1/(10d)$. By our choice of $T_{2}$, the overall estimate output at least $1-\varepsilon$ in the first case and a value at most $\varepsilon$ in the second case.

The polynomial which is the concatenation of $T_{2}$ and $T_{1}$ has $O\left(\left(dk\cdot\text{deg}(T_{1})\right)^{\text{deg}(T_{2})}\right)=(dM)^{O(M\sqrt{d}\log(Md))}$ monomials, if we consider the expression as a polynomial in the variables $x(1),\ldots,x(d)$. Our final estimator will be the sum across all $i\geq 1$. Following our upper bound techniques for matrix-vector products for polynomial, e.g. in Section A.2, and as outlined in the prior section, we get the following overall query time: