Year: 2021

Abstract:
Product quantization (PQ) coupled with a space rotation, is widely used in modern approximate nearest neighbor (ANN) search systems to significantly compress the disk storage for embeddings and speed up the inner product computation. Existing rotation learning methods, however, minimize quantization distortion for fixed embeddings, which are not applicable to an end-to-end training scenario where embeddings are updated constantly. In this paper, based on geometric intuitions from Lie group theory, in particular the special orthogonal group SO(n), we propose a family of block Givens coordinate descent algorithms to learn rotation matrix that are provably convergent on any convex objectives. Compared to the state-of-the-art SVD method, the Givens algorithms are much more parallelizable, reducing runtime by orders of magnitude on modern GPUs, and converge more stably according to experimental studies. They further improve upon vanilla product quantization significantly in an end-to-end training scenario.

text: Product quantization (PQ) coupled with a space rotation, is widely used in modern approximate nearest neighbor (ANN) search systems to significantly compress the disk storage for embeddings and speed up the inner product computation. Existing rotation learning methods, however, minimize quantization distortion for fixed embeddings, which are not applicable to an end-to-end training scenario where embeddings are updated constantly. In this paper, based on geometric intuitions from Lie group theory, in particular the special orthogonal group SO(n), we propose a family of block Givens coordinate descent algorithms to learn rotation matrix that are provably convergent on any convex objectives. Compared to the state-of-the-art SVD method, the Givens algorithms are much more parallelizable, reducing runtime by orders of magnitude on modern GPUs, and converge more stably according to experimental studies. They further improve upon vanilla product quantization significantly in an end-to-end training scenario.

heading: 1 INTRODUCTION
text: Search index is the core technology to enable fast information retrieval in various modern computational systems, such as web search, e-commerce search, recommendation and advertising in the past few decades. As a traditional type of search index, inverted index (Dean, 2009), which maps terms to documents in order to retrieve documents by term matching, have been the mainstream type of search index for decades, thanks to its efficiency and straightforward interpretability. Recently, with the advent of deep learning era, embedding indexes coupled with approximate nearest neighbor (ANN) search algorithms, have established as a promising alternative to search index (Zhang et al., 2020) and recommendation index (Covington et al., 2016; Huang et al., 2020; Li et al., 2019), in part due to its learnable representations and efficient ANN algorithms.
The main idea of embedding indexes is to encode users (or queries) and items in a latent vector space, and represent their semantic proximity in terms of inner product or cosine similarity. Embedding indexes enjoy a few appealing characteristics: a) the embeddings can be learned to optimize downstream retrieval task of interests, and b) items can be efficiently retrieved within tens of milliseconds. The latter leverages decades of algorithmic innovations, including a) maximum inner product search (MIPS) or approximate nearest neighbors (ANN), such as locality sensitive hashing (LSH) (Datar et al., 2004), hierarchical navigable small world graphs (HNSW) (Malkov & Yashunin, 2020), space indexing by trees (Bentley, 1975; Dasgupta & Freund, 2008; Muja & Lowe, 2014; Bernhardsson, 2018) or graphs (Harwood & Drummond, 2016; Iwasaki, 2015), and b) state-of-the-art product quantization (PQ) based approaches (Jegou et al., 2010; Johnson et al., 2019; Guo et al., 2020; Wu et al., 2017). In particular, PQ based indexes and its variants have regularly claimed top spots on public benchmarks such as GIST1M, SIFT1B (Jegou et al., 2010) and DEEP10M (Yandex & Lempitsky, 2016). Moreover, along with the open-sourced libraries Faiss (Johnson et al., 2019) and ScaNN (Guo et al., 2020), PQ based embedding indexes are widely adopted in many industrial systems (Huang et al., 2020; Zhang et al., 2020).

heading: 1.1 PRODUCT QUANTIZATION (PQ) BASED INDEXES
text: Since Jegou et al. (2010) first introduced PQ and Asymmetric Distance Computation (ADC) from signal processing to ANN search problem, there have been multiple lines of research focused on improving PQ based indexes. We briefly summarize the main directions below.
Coarse Quantization, also referred to as inverted file (IVF), is introduced in the original work of Jegou et al. (2010) to first learn a full vector quantization (referred to as coarse quantization) by k-means clustering and then perform product quantization over the residual of the coarse quantization. This enables non-exhaustive search of only a subsets of the clusters, which allows ANN to retrieve billions of embeddings in tens of milliseconds. More work has been done in this direction, including Inverted Multi-Index (IMI) (Lempitsky, 2012);
Implementation Optimization efforts have mainly been spent on the computation of ADC, including using Hamming distance for fast pruning (Douze et al., 2016), an efficient GPU implementation of ADC lookup (Johnson et al., 2019), and SIMD-based computation for lower bounds of ADC (André et al., 2015). Our proposed method is fully compatible with all these implementation optimizations, since we are solely focused on the rotation matrix learning algorithms.
Rotation Matrix Learning is introduced to reduce the dependencies between PQ subspaces, since PQ works best when the different subspaces are statistically independent, which may not be true in practice. As early as 2013, Optimized PQ (OPQ) (Ge et al., 2013), ITQ (Gong et al., 2013) and Cartesian k-means (Norouzi & Fleet, 2013) all propose the idea of alternating between learning (product) vector quantization and that of the rotation matrix; the latter can be formulated as the classic Orthogonal Procrustes problem (Schönemann, 1966) that has a closed form solution in terms of singular value decomposition (SVD), which is also widely used in CV (Levinson et al., 2020).
Cayley Transform (Cayley, 1846) recently inspires researcher to develop end-to-end algorithm to learn rotation matrix in various neural networks, especially unitary RNN (Helfrich et al., 2018; Casado & Martínez-Rubio, 2019) and PQ indexes (Wu et al., 2017). Formally, Cayley transform parameterizes a d× d rotation matrix R by R = (I −A)(I +A)−1, where A is a skew-symmetric matrix (i.e., A = −A>). Note that the above parameterization of R is differentiable w.r.t. the ( d 2
) parameters of A. Thus, it allows for end-to-end training of rotation matrix R, as long as the gradient of R can be obtained. However, as we will discuss in Section 2.4, these Cayley transform based methods can not be easily parallelizable to modern GPUs thus inefficient compared to our proposed method. Another drawback of Cayley transforms is numerical instability near orthogonal matrices with -1 eigenvalues. This has been partially addressed in (Lezcano-Casado & Martínez-Rubio, 2019) using approximate matrix exponential methods. However the latter incurs higher computation cost per step, due to the higher order rational functions involved. The question of minimal Euclidean dimension to embed SO(n) has been shown in (Zhou et al., 2019) to be of order O(n2). By contrast, our proposed GCD methods only require O(n) dimension of local parameterization.

heading: 1.2 TRAINABLE INDEXES
text: Despite its various advantages, a notable drawback of embedding indexes is the separation between model training and index building, which results in extra index building time and reduced retrieval accuracy. Typically, for a large industrial dataset with hundreds of millions of examples, it takes hours to build the index, and recall rates drop by at least 10% (Zhang et al., 2021). This spurred a new trend of replacing embedding indexes with jointly learnable structural indexes. Researchers argue that product quantization based embedding indexes are suboptimal, as it can not be learned jointly with a retrieval model. To overcome this drawback, researchers propose a few learnable index structures, such as tree-based ones (Zhu et al., 2018; 2019), a K-D matrix one (Gao et al., 2020) and an improved tree-based one (Zhuo et al., 2020). Though these alternative approaches have also shown improved performance, they often require highly specialized approximate training techniques, whose complexities unfortunately hinder their wide adoptions.
On the other hand, learning PQ based indexes and retrieval model in an end-to-end scenario is also challenging, in large part due to the presence of non-differential operators in PQ, e.g., arg min. Recently, however, Zhang et al. (2021) proposes an end-to-end training method which tackles the non-differentiability problem by leveraging gradient straight-through estimator (Bengio et al., 2013), and it shows improved performance over the separately trained embedding indexes. However, learning
a rotation matrix in an end-to-end training scenario remains unaddressed. As alluded to earlier, all existing methods of learning rotation matrix are based on alternative minimization procedure with an expensive non-parallelizable SVD step in each iteration. Thus it is not applicable to various back-propagation based techniques in neural network learning.
Note that the orthonormality constraint on the rotation matrix is critical to ensure pairwise distances are preserved between query and item embeddings. In end to end training, however, one could relax this condition by using an L2 regularization loss of the form ‖XX> − I‖22. We have not explored this latter approach here, since we aim for applicability of the method even for standalone PQ.

heading: 1.3 OUR CONTRIBUTIONS
text: In this paper, we overcome the last hurdle to fully enable end-to-end training of PQ based embedding index with retrieval models, by leveraging mathematical studies of decomposition of orthogonal group by Hurwitz (1963) as early as 1963. Based on intuitions about maximal tori in compact Lie groups (Hall, 2015), we propose an efficient Givens coordinate gradient descent algorithm to iteratively learn the rotation matrix – starting from an identity matrix and applying a set of maximally independent (mutually commuting) Givens block rotations at each iteration. We open sourced the proposed algorithm, which can be easily integrated with standard neural network training algorithms, such as Adagrad and Adam.
Our contribution can be summarized as follows.
• Methodologically, we change the landscape of learning rotation matrix in approximate nearest neighbor (ANN) search from SVD based to iterative Givens rotation based, in order to be applicable to end-to-end neural network training.
• Algorithmically, we propose a family of Givens coordinate block descent algorithms with complexity analysis and convergence proof of the least effective variant, GCD-R.
• Empirically, we prove that for fixed embeddings the proposed algorithm is able to converge similarly as the existing rotation matrix learning algorithms, and for end-to-end training the proposed algorithm is able to learn the rotation matrix more effectively.

heading: 2 METHOD

heading: 2.1 REVISIT TRAINABLE PQ INDEX
text: Figure 1 illustrates an overview of an embedding indexing layer inside a typical deep retrieval twotower model. Formally, the indexing layer defines a full quantization function T : Rm×n −→ Rm×n that maps a batch of m input embedding X in n-dimensional space to output embeddings T (X), which can be decomposed into two functions: a product quantization function φ : Rm×n −→ Rm×n and a rotation function with an orthonormal matrix R. The indexing layer multiplies the input embedding X by a rotation matrix R and the product quantized embedding φ(XR) by its inverse matrix R−1, which is equal to its transpose R> for orthonormal rotation matrix. Formally,
T (X) = φ(XR)R>.
We omit the definition of product quantization function φ here, since it is not the focus of this paper. Interested readers can find the full definitions in (Jegou et al., 2010) or (Zhang et al., 2021).
The final loss function includes two parts, a retrieval loss and a quantization distortion loss, formally,
L(X) = Lret(T (X)) + (1/m)‖XR− φ(XR)‖2 (1)
where Lret denotes any retrieval loss – typically by softmax cross entropy loss or hinge loss (Covington et al., 2016; Zhang et al., 2020; Huang et al., 2020). Note that straight-foward optimization of L is intractable, since the quantization function φ is not differentiable, which, however, can still be tackled by gradient straight-through estimator (Bengio et al., 2013). Then the only remaining unaddressed problem is how to iteratively learn the rotation matrix R, since existing methods only solve the case where the input vectors X are fixed (Ge et al., 2013; Jegou et al., 2010). We will present our approach in the following sections.

heading: 2.2 GIVENS ROTATION
text: Before we present the algorithm of learning a rotation matrix, let us first present a few mathematical preliminaries that serve as foundations to our algorithm. For simplicity of exposition, the embedding dimension n = 2m will always be an even number.
Definition 1. Let Rn×n denote the set of n×n matrices. An orthogonal groupO(n) = {A ∈ Rn×n : AA> = In} is the group of all distance-preserving transformations of the Euclidean space Rn. The special orthogonal group, denoted SO(n), consists of all orthogonal matrices with determinant 1.
It is straight-forward to check that axioms of group theory are satisfied by O(n) and SO(n) (see any introductory text on Lie groups, e.g., (Hall, 2015)). In fact they form Lie sub-groups of the group of n× n invertible matrices GL(n), both of which are manifolds embedded naturally in Rn2 . Definition 2. A Givens rotation, denoted as Ri,j(θ), defines a rotational linear transformation around axes i and j, by an angle of θ ∈ [0, 2π). Formally,
Ri,j(θ) := i j  Ii−1 cos θ − sin θ i
Ij−i sin θ cos θ j
In−j
where Ik denotes the k× k identity matrix. In other words, the Givens rotation matrix Rij(θ) can be regarded as an identity matrix with four elements at (i, i), (i, j), (j, i), (j, j) replaced.
Hurwitz (1963) showed the following important result as early as in 1963. The essential idea is similar to Gram-Schmidt process, but specialized to rotational matrices.
Lemma 1. Every element in SO(n) can be written as a product ∏
1≤i<j≤nRi,j(θi,j), though the decomposition is not unique.
In our case, to iteratively learn the rotation, we want to take θ to be small, such as uniformly in (− , ), at each step. While most existing work considers a single Givens rotation at each step, to our best knowledge no published work has considered multiple rotations in one steps, namely with a step rotation composed of n/2 independent rotations:
R = Rσ(1),σ(2)(θ1,2)Rσ(3),σ(4)(θ3,4) . . . Rσ(n−1),σ(n)(θn−1,n),
where σ is a randomly chosen permutation, σ(i) is the value at position i, and θi,j’s are again independent and uniform in (− , ). The fact that these rotation axes are mutually disjoint ensures that their product can be computed in a parallel fashion. Also the convergence proof of the coordinate descent algorithm relies on estimating the amount of descent in the tangent space, which in turns requires that the rotations commute with one another. This is naturally satisfied by the choice of rotations above.
Lemma 2. The n/2 families of Givens rotations Rσ(1),σ(2)(θ1,2), . . . , Rσ(n−1),σ(n)(θn−1,n) generate a maximally commuting subgroup within SO(n).
This follows from the observation that a maximal torus of SO(n) is precisely given by simultaneous rotations along n/2 pairs of mutually orthogonal 2d-planes (see (Bröcker & Tom Dieck, 2013) Chapter 5). For odd n = 2m+ 1, simply replace n/2 by (n− 1)/2. Note that both groups SO(n) and GL(n) are non-commutative for all n > 2, just as matrix multiplication is in general non-commutative. Unlike GL(2), however, SO(2) is commutative, since it consists of the 1-parameter family of planar rotations.

heading: 2.3 BLOCK COORDINATE DESCENT ALGORITHM
text: The overall objective function, as shown in Eq. (1), needs to be optimized with respect to R and other parameters in T , formally,
minimize R∈SO(n),T
L(X) = Lret(T (X)) + (1/m)||XR− φ(XR)||2, (2)
where we omit other parameters in T since they are not the focus of this paper, though clearly optimizing those parameters needs special techniques such as gradient straight-through estimator (Bengio et al., 2013). Note that, in earlier PQ work, R is either fixed to be In, the identity matrix (Chen et al., 2020), or separately optimized in a non-differentiable manner, e.g., using singular value decomposition (Ge et al., 2013), which is equivalent to optimizing only the second term in Eq. (2) with respect to only R and fixing all other parameters.
These earlier approaches do not apply to our case, since we would like to optimize the rotation matrix iteratively and jointly with other parameters. The naive update rule R← R− α∇RL, however, does not preserve the requirement that R ∈ SO(n). A straight-forward remedy is to project the update back to SO(n), such as with SVD or the tangent space exponential map, both of which however require essentially a diagonalization procedure, which is prohibitively expensive. Instead we take a direct approach to ensure that R stays on SO(n) after every iteration.
Consider the Hurwitz parameterization in Lemma 1, since θij ∈ R, we can treat it as Euclidean optimization, thus a step of gradient descent looks like R← ∏ 1≤i<j≤nRij(θij + α∇θijL). However, full Hurwitz product requires Ω(n2) sparse matrix multiplications, which is too expensive. Furthermore the various θij are non-trivially correlated. To address both, we consider selecting a subset of Givens rotation Rij at each step, according to the directional derivative as follows. Proposition 1. Given an objective function L, the (un-normalized) directional derivative d dθ ∣∣ θ=0 L(XRij(θ)), 1 ≤ i < j ≤ n, is given by 〈R′ij(0), X>∇L(X)〉 = Trace(X>∇L(X)R′ij(0)>) = (∇L(X)>X −X>∇L(X))ij . where (.)ij denotes the element at i-th row and j-th column, and 〈·, ·〉 is the Frobenius inner product.
Proof. Consider the map h : R→ Rn×n, h(θ) = XRij(θ). By linearity, h′(0) = XR′ij(0). Thus, the chain rule gives
d
dθ ∣∣∣∣ θ=0 L(XRij(θ)) = 〈∇L(X), XR′ij(0)〉.
We would like to group ∇L(X) and X together. Using the fact that 〈u, v〉 = u>v, the right hand side can be written as
∇L(X)>XR′ij(0) = (X>∇L(X))>R′ij(0) = 〈X>∇L(X), R′ij(0)〉. (3)
From Definition 2 of Givens rotation Rij(θ), we see that
R′ij(0) = i j  0i−1 cos′(0) − sin′(0) i
0j−i sin′(0) cos′(0) j
0n−j
= i j  0i−1 0 −1 i
0j−i 1 0 j
0n−j
.
Algorithm 1 Bipartite matching in greedy Givens coordinate descent (GCD-G). 1: Input: ( n 2 ) directional derivatives gij := ddθL(RRij(θ)), 1 ≤ i < j ≤ n.
2: Output: n/2 disjoint coordinate pairs {(i`, j`) : 1 ≤ ` ≤ n/2}. 3: Let S := {1, . . . , n} be the set of all coordinate axes. 4: for 1 ≤ ` ≤ n/2 do 5: Among indices in S, pick (i`, j`), i` < j`, that maximizes ‖gij‖2. 6: Remove i` and j` from S: S → S \ {i`, j`}. 7: end for
It follows that 〈R′ij(0), Y 〉 = (Y > − Y )ij for any square matrix Y , which completes the proof by plugging in Eq. (3). Finally to normalize the directional derivative, it needs to be divided by ‖R′ij(0)‖ = √ 2.
Given the above directional derivatives for each Rij(θ), we have a few options to select a subset of n/2 disjoint (i, j) pairs for Rij , listed below in order of increasing computational complexities, where GCD stands for Givens coordinate descent.
• GCD-R randomly chooses a bipartite matching and descents along the n/2 Givens coordinates, similar to stochastic block coordinate gradient descent (Beck & Tetruashvili, 2013). This can be done efficiently by first shuffling the set of coordinates {1, . . . , n}.
• GCD-G sorts the ij pairs according to the absolute value of its directional derivatives in a decreasing order, and greedily picks pairs one by one if they forms a bipartite matching. Since this is our main proposed algorithm, we list the details in Algorithm 1.
• GCD-S finds a bipartite matching with maximum sum of edge weights, giving the steepest descent. The fastest exact algorithm (Kolmogorov, 2009) for maximal matching in a general weighted graph, has a O(n3) running time, which is impractical for first-order optimization.
With any of the above options to pick n/2 disjoint ij pairs, we can present the full algorithm in Algorithm 2. Moreover, it is worth mentioning that picking non-disjoint ij pairs neither guarantees theoretical convergence as shown in Section 2.5, nor works well in practice (see Section 3.1).

heading: 2.4 COMPLEXITY ANALYSIS
text: One iteration of Givens coordinate descent, as shown in Algorithm 2, consists of three major steps: a) computing directional derivatives A in O(n3) but parallelizable to O(n); b) selecting n/2 disjoint pairs in O(n) by random, O(n2 log n) by greedy, or O(n3) by the exact non-parallelizable blossom algorithm (Kolmogorov, 2009); c) applying rotation in O(2n2) FLOPs by a matrix multiplication between R and an incremental rotation, which is a sparse matrix with 2n nonzero elements. Thus, the sequential computational complexities of one iteration for the three proposed algorithms are O(n3 + n+ 2n2) = O(n3) for GCD-R, O(n3 + n2 log n+ n2 + 2n2) = O(n3) for GCD-G, and O(n3 + n3 + 2n2) = O(n3) for GCD-S.
Though the time complexity for SVD and Cayley transform’s gradient computation is also O(n3), the derivative computation in the proposed GCD algorithms can be fully parallelized under modern GPUs. Thus GCD-R and GCD-G have parallel runtimes of O(n2) and O(n2 log n) respectively. In contrast, the SVD computation can not be easily parallelized for general dense matrices (Berry et al., 2005), due to its internal QR factorization step (Golub & Van Loan, 1996). So is the Cayley transform gradient computation, since matrix inversion by solving a linear system can not be parallelized.

heading: 2.5 CONVERGENCE ANALYSIS
text: In this section we consider gradient descent of a loss function L : SO(n)→ R. Note that SO(n) is not a convex subset of its natural embedding space Rn2 because convex combinations of two distinct orthogonal matrices are never orthogonal.
Algorithm 2 One iteration of Givens coordinate descent algorithm for rotation matrix learning.
1: Input: current n× n rotation matrix R, loss function L : Rm×n → R, learning rate λ. 2: Get a batch of input X ∈ Rm×n and run forward pass using X ′ = XR. 3: Run backward pass to get the gradient G = ∇RL(XR) and compute A := G>R−R>G. 4: for 1 ≤ i < j ≤ n do 5: Gather Aij as the directional derivative gij := 1√2 d dθ |θ=0L(RRij(θ)) (see Proposition 1). 6: end for 7: Select n/2 disjoint pairs {(i`, j`) : l} by random, greedy (see Algo. 1) or steepest options. 8: R← R ∏n/2 `=1Ri`,j`(−λgi`,j`).
As usual, convergence results require some assumptions about the domain and the function. Before stating them, first we introduce some definitions, specialized to the space SO(n) of orthogonal matrices with determinant 1.
For all below lemmas, theorem and corollary, see Appendix for a complete proof.
Definition 3. The tangent space of SO(n) ⊂ Rn2 at the identity In is given by the space of antisymmetric matrices so(n) := {g ∈ Rn×n : g + g> = 0}. That of GL(n) is simply gl(n) := Rn2 . Definition 4. A geodesic on SO(n) is a curve of the form ξ(t) = G∗ exp(tg) for a starting point G∗ ∈ SO(n) and a direction given by some g ∈ so(n); exp is the matrix exponential. Definition 5. A function L : SO(n) → R is called geodesically convex if for any geodesic ξ, L ◦ ξ : [0, 1]→ R is a convex function.
We make two assumptions on L in order to prove global convergence results for the random Givens coordinate descent method (GCD-R) on SO(n). The convergence of GCD-G and GCD-S follows by simple comparison with GCD-R at each step. Assumption 1. The objective function L is geodesically convex on the level set of the manifold S(G0) := {G ∈ SO(n) : L(G) < L(G0)} and has a unique global minimum at G∗. Assumption 2. (Global Lipschitz) For some η > 0 and all G ∈ S(G0), g ∈ so(n):
L(G exp(g)) ≤ L(G) + 〈∇L(G), g〉+ η 2 ‖g‖2. (4)
Assumption 2 leads to the following Lemma Lemma 3. For learning rate α = 1η , the descent is sufficiently large for each step:
E[L(Gk)− L(Gk+1)|Gk] ≥ 1
(n− 1)η ‖∇L(Gk)‖2.
Lemma 3.4-3.6 of (Beck & Tetruashvili, 2013) combined with the above lemma and Assumption 1 finally yield the following Theorem: Theorem 1. Let {Gk}k≥0 be the sequence generated by the GCD-R method, and let Dn := nπ be the diameter of SO(n). With the two assumptions above,
L(Gk)− L(G∗) ≤ 1
k(D2n(n− 1)η)−1 + (L(G0)− L(G∗))−1 .
In particular, the value of L(Gk) approaches that of the minimum, L(G∗), at a sub-linear rate.
By verifying Assumptions 1 and 2 for individual functions F , we have Corollary 1. Let the joint rotational product quantization learning objective be of the form
F̃ (R, C) = F ( R>
D⊕ i=1 ui
) s.t. ui ∈ Ci ⊂ Rn/D and ‖Ci‖ = K,
where ⊕D
i=1 ui := (u1, . . . , uD) stands for vector concatenation. Consider the specialization to a convex F , such as, F (x) = 〈w, x〉 (dot product) or F (x) = 〈w,x〉‖w‖‖x‖ (cosine similarity). If we fix the value of w ∈ Rn as well as the choice of quantization components ui, gradient descent applied to the map L : R 7→ F̃ (R, C) converges to the global minimum, provided the latter is unique.

heading: 3 EXPERIMENTS
text: We evaluate the convergence performance on fixed embeddings in Section 3.1, and the overall embedding retrieval performance on end-to-end trainable embedding indexes in Section 3.2. We implement the proposed GCD algorithm in Tensorflow 1.15, and perform all experiments in a single machine with 4 Nvidia 16G V100 GPU cards.

heading: 3.1 FIXED EMBEDDINGS
text: Dataset and Setup. To make sure the proposed GCD algorithms practically converge, we first evaluate, for a given set of embeddings from SIFT1M (Lowe, 2004), whether the GCD algorithms can converge similarly as the original OPQ algorithm (Ge et al., 2013), which alternates between k-means for PQ centroid learning and SVD projection for optimal rotation. Our GCD algorithms simply replace the latter SVD step by a number (5 in this experiment) of Givens coordinate descent iterations with learning rate λ=0.0001 as shown in Algorithm 2. Similarly, Cayley method (Helfrich et al., 2018; Casado & Martínez-Rubio, 2019) also replaces the SVD steps by a number of Cayley parameter updates. All methods are evaluated by quantization distortion, which is the square of L2 distance between the original embedding and quantized embedding, averaged over all examples.
Figure 2a shows the comparative results between OPQ, Cayley method and the proposed GCD algorithms. Specifically, the proposed GCD-G and GCD-S can converge as well as the original OPQ algorithm, as the convergence curves are similar. The previously widely used method, Cayley transform, does not converge as fast as the proposed GCD algorithms, potentially due to the varying scales of the transformed gradient. It also shows the results of two other alternatives: overlapping GCD-G and overlapping GCD-R, which do not enforce the (i, j) pairs to be disjoint as we discussed in Section 2.2. The comparison shows the necessity of enforcing the disjoint selection, since otherwise GCD-G does not converge well, though it does not affect GCD-R much. But GCD-R in general does not work as well as GCD-G, which indicates the necessity of picking the steeper descent directions. Moreover, we can observe that GCD-G converges similarly as GCD-S, which indicates that a greedy bipartite matching algorithm as given in Algorithm 1 may be good enough in practice.
Figures 2b and 2c further compare the convergence performance between GCD-G and OPQ, by averaging over 10 runs and varying the size of data. We can make two important observations: a) GCD-G converges much more stably than OPQ, as the curves show significantly lower variances; b) GCD-G converges better for smaller size of data, e.g., for 10% of the data. It indicates that GCD-G potentially works better in the stochastic descent scenario where the batch size is usually small.

heading: 3.2 END-TO-END TRAINABLE EMBEDDING INDEXES
text: Dataset and Setup. In this section, we evaluate GCD-S and GCD-R performances in an industrial search click log data where a user input query is used to retrieve items, and two public datasets of MovieLens (Harper & Konstan, 2015) and Amazon Books (He & McAuley, 2016) where user historical behaviors are used to retrieve next behavior. The industrial dataset contains 9,989,135 training examples with 1,031,583 unique queries and 1,541,673 unique items, which is subsampled about 2% from our full dataset, but enough for training a reasonably well search retrieval model. We evaluate the effectiveness of the learned rotation matrix by quantization distortion, precision@k
(p@k) and recall@k (r@k) metrics, which are standard retrieval quality metrics. In our experiments, we choose k=100. Specifically, for each query, we retrieve from embedding indexes a set of top k items, which is then compared with the ground truth set of items to calculate the precision and recall. We perform ANOVA test to evaluate the p-value for all comparisons. We implement a two-tower retrieval model (Zhang et al., 2021) with cosine scoring, hinge loss of margin 0.1 and embedding size 512. At the beginning, a number of 10,000 warmup steps are performed before applying the indexing layer as shown in Figure 1. Then a number of 8,192 examples are collected to warm start the PQ centroids and rotation matrix with 200 iterations of OPQ algorithm. The baseline method uses this initial rotation matrix and stops updating it, while the other methods continue updating it.
As shown in Figure 3, in all three datasets, the proposed GCD algorithms and Cayley method clearly outperform the baseline in terms of quantization distortion, which consequentially translates into the improvement in retrieval quality metrics, measured by precision@100 and recall@100 as shown in Table 1. However, the Amazon Book dataset is an exception. We suspect the dataset may have subtle distribution, thus the small quantization distortion reduction does not result in improvements on retrieval quality metrics. Finally, we observe that the three proposed GCD algorithms consistently with theoretical analysis – GCD-R and GCD-S can be regarded as the lower bound and upper bound of GCD methods, with GCD-G in between. The GCD-S with steepest Givens rotations picked at each step, generally shows the best performance in both quantization distortion and retrieval metrics, especially better than the widely known Cayley method in literature.

heading: 4 CONCLUSION
text: In this paper, we have proposed a novel method called Givens coordinate descent (GCD) to iteratively learn a rotation matrix, in the context of end-to-end trainable PQ based embedding indexes. Based on a well studied mathematical result that any orthogonal matrix can be written as a product of Givens rotations, the proposed GCD algorithm picks n/2 maximally commuting Givens rotations to update the learning rotation matrix by their directional derivatives. We also provide complexity analysis and prove convergence of the proposed algorithm under convexity and other mild assumptions. Experimental results show that the proposed method converges more stably than previous SVD based and Cayley transform method, and it significantly improves retrieval metrics in the end-to-end training of embedding indexes. Our algorithm is also applicable to other scenarios in any neural network training where learning a rotation matrix is potentially helpful.

heading: A CONVERGENCE RESULTS
text: In this section we prove the results stated in Section 2.5. Lemma 1. For learning rate α = 1L , the descent is sufficiently large for each step:
E[F (xk)− F (xk+1)|xk] ≥ 1
2(n− 1)L ‖∇F (xk)‖2.
Proof. Recall that each step of RGCGD is given by xk+1 = xk · G where G =∏n/2 i=1 gσ(2i−1),σ(2i)(θi), where σ ∈ Sn is a uniformly random permutation. Furthermore since the gσ(2i−1),σ(2i)’s in the product are along disjoint axes, they commute with one another:
gσ(2i−1),σ(2i)(θi)gσ(2j−1),σ(2j)(θj) = gσ(2j−1),σ(2j)(θj)gσ(2i−1),σ(2i)(θi), (5)
we can write them as gi,j(θ) = exp( θ√2 (ei ∧ ej)), where the √
2 factor accounts for the fact that ‖ei ∧ ej‖2 = 2 under Euclidean norm. Thus let’s define s ∈ so(n) := TISO(n) by
s := n/2∑ i=1 θi√ n eσ(2i−1) ∧ eσ(2i),
xk+1 = Rxk(αs) = xk · exp(αs) = xk · exp(α n/2∑ i=1 θi√ 2 eσ(2i−1) ∧ eσ(2i)).
In fact, the θi’s are projections of the full gradient∇F̃ (x), for F̃ (x) := F (xk · x),
θi,j := −〈∇F̃ (I), 1√ 2 ei ∧ ej〉, for 1 ≤ i < j ≤ n.
So from (4), and orthogonality of 1√ 2 eσ(2i−1) ∧ eσ(2i), we have
F (xk+1) ≤ F (xk)− α‖s‖2 + Lα2
2 ‖s‖2. (6)
Now choose α = 1L . We easily get
F (xk+1)− F (xk) ≤ − 1
2L ‖s‖2. (7)
Finally recall that s ∈ so(n) is a random vector with randomness given by the uniform σ ∈ Sn. By orthogonality and the definition of θi above, we have
E‖s‖2 = 1 n! ∑ σ∈Sn n/2∑ i=1 θ2σ(2i−1),σ(2i) = 1 n− 1 ∑ 1≤i<j≤n θ2i,j = ‖∇F̃ (0)‖2 n− 1 = ‖∇F (xk)‖2 n− 1 . (8)
Taking expectation of both sides of (7):
E[F (xk+1)|xk] ≤ F (xk)− 1
2(n− 1)L ‖∇F (xk)‖2. (9)
Next we recall Lemma 3.5 from Beck & Tetruashvili (2013) Lemma 2. Let {Ak}k≥0 be a non-negative sequence of real numbers satisfying Ak −Ak+1 ≥ ηA2k, k = 0, 1, . . ., and A0 ≤ 1mη . for positive η and m. Then
Ak ≤ 1
η(k +m) , for k = 0, 1, . . . (10)
Ak ≤ 1
ηk , for k = 1, 2, . . . (11)
This along with Assumption 1 leads to global convergence result below (combining Lemma 3.4 and Theorem 3.6 of Beck & Tetruashvili (2013)): Theorem 1. Let {xk}k≥0 be the sequence generated by the RGCGD method. With the two assumptions above,
F (xk)− F (x∗) ≤ 1
2kR2(x0)(n− 1)L+ (F (x0)− F (x∗))−1 . (12)
In particular, the value of F (xk) approaches that of the minimum, F (x∗), at a sub-linear rate.
Proof. We can write xk = Rx∗(s), where s = γ′(0) and γ : [0, 1]→ S(x0) ⊂ SO(n) is a constant speed geodesic path from x∗ to xk. This means that ‖γ′(t)‖ is constant for all t. By geodesic convexity, F ◦ γ : [0, 1]→ R is convex, therefore
F (xk)− F (x∗) ≤ (F ◦ γ)′(0)〈∇F (xk), γ′(0)〉.
Since γ′(0) = d(xk, x∗), the geodesic distance between xk and x∗, by the constant speed parametrization of γ, Cauchy-Schwarz inequality then implies
F (xk)− F (x∗) ≤ ‖∇F (xk)‖d(xk, x∗) ≤ R(x0)‖∇F (xk)‖.
This combined with Lemma 1 gives
E[F (xk)− F (xk+1|xk] ≥ 1
2R2(x0)(n− 1)L (F (xk)− F (x∗))2.
Now consider the Taylor expansion around the global minimum x∗. Since the gradient vanishes, the Lipschitz condition (4) simplifies to
F (xk)− F (x∗) ≤ L‖s‖2 2 ≤ L 2 R2(x0).
For the last inequality, we again use γ′(0) = d(xk, x∗) ≤ R(x0) by definition of the latter. So letting Ak = E[F (xk) − F (x∗)], η = 12R2(x0)(n−1)L , and m = 1 η(F (x0)−F (x∗)) , Lemma 2 implies that
F (xk)− F (x∗) ≤ 1
(k +m)η =
1
2kR2(x0)(n− 1)L+ (F (x0)− F (x∗))−1 .
Corollary 1. Let the joint rotational product quantization learning objective be given by
L(w,R, C) = F (w)(R D⊕ i=1 ui) s.t. ui ∈ Ci ⊂ Rn/D and ‖Ci‖ = K. (13)
Consider the specialization of w ∈ Rn and F (w)(x) = 〈w, x〉 (dot product) or F (w)(x) = 〈w,x〉‖w‖‖x‖ (cosine similarity). If we fix the value of w ∈ W as well as the choice of quantization components ui, gradient descent applied to the map F̃ : R 7→ L(w,R,C) converges to the global minimum, provided the latter is unique.
Proof. In the dot product case, F̃ (R) = 〈w,Rx〉. If there is a unique global minimizer, it can be given by UV ⊥ from the following singular value decomposition:
wx⊥ = USV ⊥.
First we verify Assumption 1, namely that the following function is convex:
t 7→ F̃ (G0 exp(tA)) = 〈w,G0 exp(tA)x〉.
Indeed this follows immediately from diagonalizing A, which yields a linear combination of convex functions of t.
For assumption 2, recall we need to show F̃ (G0 exp(A)) ≤ F̃ (G0) + 〈∇F̃ (G0), A〉+ L2 ‖A‖ 2, for some L > 0. By differentiating 〈w,G0 exp(tA)x〉 at t = 0, we see that
d dt |t=0F̃ (G0 exp(tA)) = 〈∇F̃ (G0), A〉
F̃ (G0 exp(A))− F̃ (G0) = 〈w,G0(exp(A)− I)x〉 = 〈w,G0(A+ 1
2 A2 + . . .)x〉
= 〈∇F̃ (G0), A〉+ 1
2 〈∇F̃ (G0), A2) ∞∑ k=2 1 k! 〈w,Akx〉.
For ‖A‖ small, the last term can be bounded by e‖A‖ − 1− ‖A‖, which is O(‖A‖2). For larger ‖A‖, say > C, we can choose L large enough so that LC 2
2 > supF .
The cosine similarity case follows immediately from the dot product case since the denominator is a constant for fixed w and ⊕D i=1 ui, thanks to the fact that ‖Rx‖ = ‖x‖.
Corollary 2. Under the two assumptions 1 and 2, the steepest Givens coordinate gradient descent (SGCGD) method also converges to its unique global minimum.
Proof. By the choice of coordinates, each step of SGCGD is dominated by the corresponding step of RGCGD based at the same point x ∈ SO(n), in other words, the estimate (12) continues to hold.

heading: B EMPIRICAL RUNNING TIME
text: Figure 4 shows the comparison of per step runtime cost between GCD methods and Cayley transform for standalone QP training, with batch size chosen to be 1. In Figure 4a, for a fair comparison, we only compare the GCD-R method with Cayley transform, since we have not implemented GCD-G and GCD-S fully in GPU operators yet. Note that the GCD-G method needs a sorting operator implemented in GPU (see step 5 in Algorithm 1), which is possible but needs a customized GPU operator implemented in addition to the original TensorFlow release. On the other hand for Cayley transform, while the gradient propagation and matrix multiplication are all computed on GPU, the matrix inversion operation unfortunately cannot be parallelized, thus effectively runs on CPU. As expected from the complexity analysis in Section 2.4, GCD-R slows down quadratically in the embedding dimension n, whereas Cayley incurs an asymptotic O(n3) cost per step. To further explore the empirical computation cost in a completely fair setup, we compare the methods all in CPU. As shown in Figure 4b, GCD-G and GCD-R run much faster than the Cayley method, with about 5 times and 50 times faster running speed, respectively, which roughly reflect the total number of floating point computations in each method.

Review summaries:
This paper proposes a block coordinate descent algorithm for rotation learning. The algorithm is based on Lemma 1 and Theorem 1. The rotation matrix on SO(n) is decomposed into diverse simple Givens rotation matrices. Then the optimized variable is converted into these Givens rotation matrices so that the rotation matrix is always on SO(3) and the projection is not required anymore. The authors also discuss how to select the coordinate, including random strategy, greedy strategy, and steepest strategy. Different from the existing work, it considers multiple Given rotations matrices in one step. 

Current rotation learning methods are trying to minimize quantization distortion for fixed embeddings, which are not applicable to an end-to-end training scenario where embeddings are getting updated constantly. Therefore, this paper tries to address this issue to fully enable end-to-end training of Product Quantization (PQ) based embedding index with retrieval models, by using mathematical studies of the decomposition of orthogonal group. They proposed a family of block Givens coordinate descent algorithms to learn rotation matrices that are provably convergent on any convex objectives by leveraging geometric intuitions from Lie group theory. Authors claimed that their algorithms are much more parallelizable, reducing runtime by orders of magnitude on modern GPUs, and converge more stably according to experimental studies in comparison to the state-of-the-art SVD method.   Their main contributions can be summarized as follows:  - Changing the landscape of learning rotation matrix in approximate nearest neighbor (ANN) search from SVD based to iterative Givens rotation-based, to be applicable to end-to-end neural network training.  - Proposing a family of Givens coordinate block descent algorithm with complexity analysis and convergence proof.  - Proves that for the fixed embedding, their algorithm shows similar convergence result as the existing rotation matrix learning algorithms. Therefore, their proposed algorithm is able to learn the rotation matrix more effectively for the end-to-end training.

This paper proposes to learn rotation matrix by Givens coordinate descent algorithms in the context of minimizing the quantization distortion for efficient storage. The proposed family of Givens coordinate descent algorithms are based on geometric intuitions of the special orthogonal group and are provably convergent on any convex objectives. The experiments show that the proposed algorithms are much time efficient and lead to performance improvement in an end-to-end training of embedding indexes.

The paper address rotation matrix learning during product quantization in modern ANN embedding search systems. The main contribution is addressing rotation matrix learning via gradient descent of small rotation updates. The approach relies on the decomposition of any small rotation matrix into a product of Given rotations, so that partial derivatives can be obtained in parallel, but the product causes the computation of the rotation matrix itself to be slow, in O(n^2) matrix multiplications, hence the need to select a subset of coordinates and do coordinate descent. Experiments on product quantization show a marginal improvement in results over existing approaches OPQ and Cayley.