Sparse Quadratic Optimisation over the Stiefel Manifold with Application to Permutation Synchronisation
Year: 2021

Abstract:
We address the non-convex optimisation problem of finding a sparse matrix on the Stiefel manifold (matrices with mutually orthogonal columns of unit length) that maximises (or minimises) a quadratic objective function. Optimisation problems on the Stiefel manifold occur for example in spectral relaxations of various combinatorial problems, such as graph matching, clustering, or permutation synchronisation. Although sparsity is a desirable property in such settings, it is mostly neglected in spectral formulations since existing solvers, e.g. based on eigenvalue decomposition, are unable to account for sparsity while at the same time maintaining global optimality guarantees. We fill this gap and propose a simple yet effective sparsity-promoting modification of the Orthogonal Iteration algorithm for finding the dominant eigenspace of a matrix. By doing so, we can guarantee that our method finds a Stiefel matrix that is globally optimal with respect to the quadratic objective function, while in addition being sparse. As a motivating application we consider the task of permutation synchronisation, which can be understood as a constrained clustering problem that has particular relevance for matching multiple images or 3D shapes in computer vision, computer graphics, and beyond. We demonstrate that the proposed approach outperforms previous methods in this domain.

heading: 1 Introduction
text: We are interested in optimisation problems of the form
arg max U∈St(m,d)
f(U) with f(U) = tr(UTWU) and St(m, d) := {X ∈ Rm×d : XTX = Id} , (1)
where Id is the identity matrix of dimension d, W ∈ Rm×m and the set St(m, d) denotes the Stiefel manifold (m ≥ d). Throughout the paper, w.l.o.g. we consider a maximisation formulation and we assume that W is a symmetric and positive semidefinite matrix (see Lemma 1 in Sec. 2). Despite the non-convexity of Problem (1), problems of this form can be solved to global optimality based on the eigenvalue decomposition, i.e. by setting U∗ = Vd, where Vd ∈ Rm×d denotes an orthogonal basis of eigenvectors corresponding to the d largest eigenvalues of W , see Lemma 2. The fact that we can efficiently find global optima of Problem (1) makes it a popular relaxation formulation of various difficult combinatorial problems. This includes spectral relaxations [43, 26, 17] of the NP-hard quadratic assignment problem (QAP) [35], spectral clustering [33, 44], or spectral permutation synchronisation [34, 39, 31]. Yet, a major difficulty of such approaches is to discretise the (relaxed) continuous solution in order to obtain a feasible solution of the original (discrete) problem. For example, in the QAP one optimises over the set of permutation matrices, in clustering one optimises over the set of matrices with rows comprising canonical basis vectors, and in permutation synchronisation one optimises over multiple permutation matrices that are stacked into a large block
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
matrix. Often such combinatorial sets can be characterised by Stiefel matrices that are sparse – for example, the set of (signed) d×d permutation matrices can be characterised by matrices in St(d, d) that have exactly d non-zero elements, whereas any other element in St(d, d) that has more than d non-zero elements is not a signed permutation matrix.
The most common approach to obtain sparse solutions in an optimisation problem is to integrate an explicit sparsity-promoting regulariser. However, a major hurdle when considering quadratic optimisation problems over the Stiefel manifold is that the incorporation of sparsity-promoting terms is not compatible with solutions based on eigenvalue decomposition, so that desirable global optimality guarantees are generally no longer maintained.
Instead, we depart from the common path of integrating explicit regularisers and instead exploit the orthogonal-invariance present in Problem (1) to assure sparsity. To be more specific, for any orthogonal matrix Q ∈ O(d) := St(d, d) it holds that f(U) = tr(UTWU) = tr(UTWUQQT ) = tr((UQ)TW (UQ)) = f(UQ). Hence, if Ū ∈ St(m, d) is a solution to Problem (1), so is ŪQ ∈ St(m, d) for any Q ∈ O(d). The subspace im(ŪQ) that is spanned by ŪQ for a given Ū ∈ St(m, d) (and arbitrary Q ∈ O(d)) is equal to the subspace im(Ū). Motivated by this observation we utilise orthogonal-invariance in order to find a solution U∗ ∈ {U ∈ St(m, d) : im(U) = im(Ū)} such that U∗ is sparse. To this end, we build upon the additional degrees of freedom due to Q ∈ O(d), which allows to rotate a given solution Ū ∈ St(m, d) to a sparser representation U∗ = ŪQ ∈ St(m, d) – most notably, while remaining a globally optimal solution to Problem (1).
Main contributions. We summarise our main contributions as follows: (i) For the first time we propose an algorithm that exploits the orthogonal-invariance in quadratic optimisation problems over the Stiefel manifold while simultaneously accounting for sparsity in the solution. (ii) Despite its simplicity, our algorithm is effective as it builds on a modification of the well-established Orthogonal Iteration algorithm for finding the most dominant eigenspace of a given matrix. (iii) Our algorithm is guaranteed to converge to the dominant subspace with the same convergence rate as the original Orthogonal Iteration algorithm, and our solution constitutes a global optimiser of Problem (1). (iv) We experimentally confirm the efficacy of our approach in the context of the permutation synchronisation problem.

heading: 2 Preliminaries & Related Work
text: In this section we clarify our assumptions, introduce additional preliminaries, and provide reference to related work. Let λ1, λ2, . . . , λm be the eigenvalues of W ∈ Rm×m ordered decreasingly. We impose the following assumption on W :
Assumption 1 (Separated eigenspace). We assume that λd > λd+1.
Throughout the paper we also assume that W is symmetric and positive semidefinite (p.s.d.), which, however, is not a restriction as the following straightforward result indicates:
Lemma 1 (Symmetry and positive semidefiniteness). In Problem (1), if W is not symmetric and not p.s.d. there is an equivalent optimisation problem (i.e. with the same optimisers) where W has been replaced by a symmetric p.s.d. matrix W̃ .
Proof. See Appendix.
Next, we define the notion of a dominant subspace and convergence of matrix sequences to such a subspace.
Definition 1 (Dominant invariant subspace). The d-dimensional dominant invariant subspace (or dominant subspace in short) of the matrix W ∈ Rm×m is defined as the subspace im(Vd) ⊆ Rm, where Vd ∈ St(m, d) is the matrix whose columns are formed by the d eigenvectors corresponding to the d largest eigenvalues of W .
Definition 2 (Convergence). We say that a sequence of matrices {Ut} converges to the dominant subspace of W if limt→∞ ‖VdV Td Ut − Ut‖ = 0.
There exists a close relation between the dominant subspace of W and solutions of Problem (1):
Lemma 2 (Solution to Problem (1)). Problem (1) is solved for any matrix U∗ that forms an orthogonal basis for the d-dimensional dominant subspace of W , i.e. U∗ ∈ {U ∈ St(m, d) : im(U) = im(Vd)}.
Proof. See Appendix.
Basic algorithms for computing eigenvectors. In order to find the dominant subspace of W we can use algorithms for finding eigenvectors. The Power method [21] is an efficient way for finding the single most dominant eigenvector of a given matrix, i.e. it considers the case d = 1. It proceeds by iteratively updating a given initial v0 ∈ St(m, 1) based on the update vt+1 ←Wvt/‖Wvt‖. In order to find the d most dominant eigenvectors, one can consider the Orthogonal Iteration algorithm [21], which generalises the Power method to the case d > 1. The algorithm proceeds by repeatedly computing Vt+1Rt+1 ←WVt based on the (thin) QR-decomposition of WVt, where Vt ∈ St(m, d) and Rt ∈ Rd×d is upper triangular. For t→∞ the sequence {Vt} converges (under mild conditions on the intial V0) to the dominant invariant subspace as long as |λd| > |λd+1|, see Thm. 8.2.2 in [21]. Sparse Stiefel optimisation. There are numerous approaches for addressing general optimisation problems over the Stiefel manifold, including generic manifold optimisation techniques (e.g. [3, 14]), or Stiefel-specific approaches (e.g. [30, 46]). In the following we will focus on works that consider sparse optimisation over the Stiefel manifold that are most relevant to our approach. In order to promote sparse solutions, sparsity-inducing regularisers can be utilised, for example via the minimisation of the (non-convex) `p-‘norm’ for 0 ≤ p < 1, or the (convex) `1-norm [36]. However, the non-smoothness of such regularisers often constitutes a computational obstacle.
The optimisation of non-smooth functions over the Stiefel manifold has been considered in [16], where the sum of a non-convex smooth and a convex non-smooth function is optimised via a proximal gradient method. In [27], the authors consider the optimisation of a sum of non-smooth weakly convex functions over the Stiefel manifold using Riemannian subgradient-type methods. Yet, in practice often differentiable surrogates of non-smooth sparsity-promoting terms are considered [40, 29, 38, 15]. Instead of minimising `p-‘norms’ with 0 ≤ p ≤ 1, on the Stiefel manifold one may instead choose the maximisation of `p-norms with p > 2, such as the `3-norm [47], or the `4-norm [49]. Further motivation for sparsity promoting higher-order norms in this context can be found in [37, 50, 28].
Optimisation over the Stiefel manifold has a close connection to optimisation over the Grassmannian manifold [3]. There are numerous approaches for Grassmannian manifold optimisation (e.g. [19, 18]), including sparse optimisation via an `1-norm regulariser [45] and the optimisation of non-convex and non-smooth objective functions via a projected Riemannian subgradient method [52]. In our case, due to the rotation-invariance of the objective f(U) = f(UQ) for any Q ∈ O(d) in Problem (1), finding the dominant subspace of W could also be posed as an optimisation problem over the Grassmannian manifold. However, we are not only interested in identifying this subspace (which can for example be done via the Orthogonal Iterations algorithm [21], or via Grassmann-Rayleigh Quotient Iterations [4]), but we want to find a specific choice of coordinates for which the representation of the subspace is sparse.

heading: 3 Proposed Method for Sparse Quadratic Optimisation over the Stiefel
text: Input: W ∈ Rm×m, U0 ∈ Rm×d Output: U∗ Initialise: t← 0 repeat
Ut+1Rt+1 ←WUtZ(Ut) // unique QR-decomposition until convergence U∗ ← Ut+1 Algorithm 1: Overview of our proposed algorithm. In terms of convergence properties our algorithm is equivalent to the Orthogonal Iteration algorithm and thus produces a U∗ in the dominant subspace of W . However, our modification introduces the matrix Z(Ut) with the purpose of promoting sparsity of the solution (see Sec. 3.2 how we choose Z(Ut)).
In this section we introduce our Algorithm 1 that can be seen as a modification of the popular Orthogonal Iteration algorithm. The main difference is that we iteratively weigh the matrix Ut ∈ St(m, d) by a matrix Z(Ut) ∈ Rd×d, where the matrix function Z maps onto the set of full rank matrices. The purpose of the matrix Z(Ut) (see Sec. 3.2 for our specific choice) is to promote sparsity. Intuitively, we characterise sparsity as having few elements that are large, whereas most elements are close to zero, which is formally expressed in (4).
First, we focus on the overall interpretation of our approach: We want to ensure that we retrieve a U∗ that is a global maximiser of Problem (1), which, in some relaxed sense, is also sparse. However, we cannot in general augment the objective function in Problem (1) with an explicit sparsity regulariser and expect that the respective solution is still a maximiser of the original objective. Instead we steer the solution towards being more sparse by weighing our matrix of interest Ut with the matrix Z(Ut).

heading: 3.1 Convergence
text: We start by ensuring that, under the assumption that Zt(Ut) is full rank, the sequence {Ut} generated by Algorithm 1 converges to the dominant subspace of W . This can straightforwardly be shown by using the following result.
Lemma 3. Consider the two algorithms:
Algorithm A: Vt+1Rt+1 ←WVtXt, Algorithm B: V̄t+1R̄t+1 ←WV̄t,
where Xt ∈ Rd×d is full rank, and for each of the two left-hand sides above the two matrices in the product are obtained by the unique (thin) QR-factorisation of the corresponding right-hand side, where Rt+1 and R̄t+1 are upper triangular with positive diagonal. Now, if Vt = V̄tQ̄t ∈ St(m, d), where Q̄t ∈ O(d), then Vt′ is equal to V̄t′ up to rotation from the right for all t′ > t.
Proof. It holds that Vt+1Rt+1 = WVtXt = WV̄tQ̄tXt = V̄t+1(R̄t+1Q̄tXt) = V̄t+1Q̃t+1R̃t+1, where Q̃t+1 and R̃t+1 are the two matrices in the QR-factorisation of the matrix (R̄t+1Q̄tXt). This means that Vt+1 = V̄t+1Q̃t+1. This is due to the uniqueness of the QR-factorisation, so that also Rt+1 = R̃t+1. Now the results follows readily by using induction.
Thus, if Zt(Ut) is full rank for each t, we can identify our Algorithm 1 with Algorithm A in Lemma 3, and conclude that the sequence {Ut} converges to the dominant subspace of W as t→∞. This is because Algorithm B in Lemma 3 corresponds to the Orthogonal Iteration algorithm, which is known to converge to the dominant subspace of W [21].
We continue by investigating the behaviour at the limit where the columns of Ut are in the dominant subspace of W . Let us consider the update Ut+1Rt+1 ← WUtZ(Ut) of Algorithm 1, and further assume that the d largest eigenvalues ofW are equal, which is for example the case for synchronisation problems assuming cycle consistency (cf. Sec. 4).
Lemma 4. Assume that the d largest eigenvalues of the p.s.d. matrix W are all equal to one1 and strictly larger than the other eigenvalues, and assume that the columns of Ut ∈ St(m, d) are contained in the dominant subspace of W . Provided Z(Ut) is full rank, it holds that
WUtZ(Ut) = UtZ(Ut) = Ut+1Rt+1 = UtQt+1Rt+1,
where Ut+1Rt+1 is the unique QR-factorisation of UtZ(Ut) and Qt+1Rt+1 is the unique QRfactorisation of Z(Ut).
Proof. Since all the d largest eigenvalues of W are equal to one and the columns of Ut are in the dominant subspace it holds that WUt = UtUTt Ut = Ut, which explains the first equality. While the second equality follows by definition, the third equality remains to be proven. Let Qt+1R̃t+1 be the unique QR-decomposition of Z(Ut). We want to show that Ut+1 = UtQt+1 and Rt+1 = R̃t+1. From the QR-decomposition of Z(Ut) it follows that UtZ(Ut) = UtQt+1R̃t+1, while from the
1If they are all equal to a value λ1 6= 1, we can w.l.o.g. consider λ−11 W in place of W , since the optimiser of Problem (1) is invariant to scaling W .
QR-decomposition of UtZ(Ut) it follows that UtZ(Ut) = Ut+1Rt+1. The uniqueness of the QRdecomposition implies UtQt+1 = Ut+1 and R̃t+1 = Rt+1.
Hence, under the assumptions in Lemma 4, the columns of Ut′ span the dominant subspace of W for any t′ ≥ t. The update in Algorithm 1 reduces to the form
Qt+1Rt+1 ← Z(Ut), (2) Ut+1 ← UtQt+1, (3)
i.e. updating Ut to obtain Ut+1 simplifies to multiplication of Ut with the Q-matrix from the QRfactorisation of Z(Ut).
3.2 Choosing Z(Ut) to Promote Sparsity
We now turn our attention to the choice of Z(Ut) in Algorithm 1. We know that, as long as Z(Ut) is full rank in each iteration t, the columns of Ut converge to the dominant subspace of W . In addition to our objective function f in Problem (1), we now introduce the secondary objective
g(U) = m∑ i=1 d∑ j=1 (Uij) p = tr(UT (U.∧(p−1))), (4)
where for p being a positive integer the notation U.∧p means to raise each element in U to the power p. For U ∈ St(m, d) and p larger than 2, the maximisation of g promotes sparsity in a relaxed sense. To be specific, if p is odd, a few larger elements (with value close to 1) lead to a larger value of g compared to many smaller elements (with value close to 0), so that sparsity and non-negativity are simultaneously promoted. Analogously, if p is even, a few elements with value closer to ±1 lead to a larger value of g compared to many smaller elements close to 0.
Let us consider the maximisation of
g(UQ) = tr(h(U,Q)) (5)
with respect to Q, where h(U,Q) = (UQ)T ((UQ).∧(p−1)), Q ∈ Od, and U ∈ St(m, d). This means that we want to rotate U by Q in such a way that the secondary objective g is maximised after the rotation. By applying the rotation from the right, we ensure that the columns U span the same space after the rotation, but a higher objective value is achieved for the secondary objective g in (4).
The extrinsic gradient with respect toQ (forQ relaxed to be in Rd×d) of g(UQ) atQ = I is p·h(U, I). Thus the (manifold) gradient of g(UQ) at Q = I is given by the projection of the extrinsic gradient at Q = I onto the tangent space, which reads
∇Qg(UQ)|Q=I = p(h(U, I)− hT (U, I)). (6) For a small step size α, we can perform a first-order approximation of gradient ascent if we choose Z = I+αp(h(U, I)−hT (U, I)), and then update U in terms of the Q-matrix of the QR-factorisation of Z. Now, let us assume that we investigate the behaviour at the limit under the assumptions in Lemma 4, i.e. that the d largest eigenvalues of W are equal. In this context the Q-matrix of the QR-factorisation is a retraction [3] and serves as a first-order approximation of the exponential map. With that, the updates become
Qt+1Rt+1 ← I + αt(h(Ut, I)− hT (Ut, I)), (7) Ut+1 ← UtQt+1. (8)
This choice of the matrix Z(Ut) ensures that it is full rank, since adding the identity matrix and a skew-symmetric matrix results in a full rank matrix. Hence, for sufficiently small step size α, we iteratively move in an ascent direction on Od by utilising the QR-retraction.

heading: 4 Application to Permutation Synchronisation

heading: 4.1 Permutation Synchronisation in a Nutshell
text: Permutation synchronisation is a procedure to improve matchings between multiple objects [23, 34], and related concepts have been utilised to address diverse tasks, such as multi-alignment [7, 5,
25], multi-shape matching [22, 24, 20], multi-image matching [51, 42, 9, 13, 12], or multi-graph matching [48, 6, 41], among many others. Permutation synchronisation refers to the process of establishing cycle consistency in the set of pairwise permutation matrices that encode correspondences between points in multiple objects. In computer vision there is a typical application where the points are feature descriptors or key-points, and the objects are images, as shown in Fig. 2.
Let k denote the number of objects, where each object i contains mi points. For 1p being a pdimensional vector of all ones, and vector inequalities being understood in an element-wise sense, let Pij ∈ Pmimj := {X ∈ {0, 1}mi×mj : X1mj ≤ 1mi ,1TmiX ≤ 1 T mj} be the partial permutation matrix that represents the correspondence between the mi points in object i and the mj points in object j. In the case of bijective permutations, the set of pairwise permutations P := {Pij}ki,j=1 is said to be cycle-consistent if for all i, j, ` it holds that Pi`P`j = Pij .
We define the set of partial permutation matrices with full row-rank as Pmid := {X ∈ Pmid : X1d = 1mi}, where d denotes the total number of distinct points across all objects. Cycle consistency is known to be equivalent to the existence of so-called object-to-universe matchings U := {Pi ∈ Pmid}ki=1 such that for all i, j we can write Pij = PiPTj (see [23, 34] for details). The object-touniverse characterisation of cycle consistency is also valid for the case of non-bijective (i.e. partial) permutations (see [42, 8] for details).
Given the noisy (i.e. not cycle-consistent) set of pairwise permutations P = {Pij}ki,j=1, permutation synchronisation can be phrased as the optimisation problem
arg max {Pi∈Pmid} ∑ i,j tr(PTijPiP T j ) ⇔ arg max P∈U tr(PTWP ) , (9)
where form := ∑
imi we define the set U := Pm1d× . . .×Pmkd ⊂ Rm×d, the (m×d)-dimensional block matrix P = [PT1 , . . . , P T k ]
T , and the block matrix W := [Pij ]ij ∈ Rm×m in order to allow for a compact matrix representation of the objective.

heading: 4.2 Proposed Permutation Synchronisation Approach
text: The core idea of existing spectral permutation synchronisation approaches [34, 39, 31] is to replace the feasible set U in Problem (9) with the Stiefel manifold St(m, d), so that we obtain an instance of Problem (1). We utilise Algorithm 1 to obtain a (globally optimal) solution U∗ ∈ St(m, d) of this spectral formulation. In our algorithm we choose p = 3 to promote sparsity and non-negativity in the resulting U∗ via the function g(U) = tr(UT (U.∧2)) = ∑m i=1 ∑d j=1 U 3 ij . With that, in addition to U∗ being an orthogonal matrix, it contains few large elements that are close to 1 and many smaller elements that are close to 0. As such, we can readily project the matrix U∗ onto the set U in terms of a Euclidean projection. The Euclidean projection is given by projecting each of the k blocks of U∗ = [U∗T1 , . . . , U ∗T k ]
T individually onto the set of partial permutations, i.e. Pi = projPmid(U ∗ i ) =
arg maxX∈Pmid tr(XTU∗i ), see e.g. [8], which amounts to a (partial) linear assignment problem [32, 11] that we solve based on the efficient implementation in [10].

heading: 4.3 Experimental Results
text: We experimentally compare our proposed approach with various methods for permutation synchronisation and perform an evaluation on both real and synthetic datasets. In particular, our comparison includes two existing spectral approaches, namely MATCHEIG [31] and SPECTRAL [34], where for the latter we use the efficient implementation from the authors of [51]. In addition, we also compare against the alternating minimisation method MATCHALS [51], and against the non-negative matrix factorisation approach NMFSYNC [8]. To emphasise that the methods MATCHEIG and MATCHALS do not guarantee cycle consistency (but instead aim to improve the initial matchings), in all plots we show results of respective methods as dashed lines. We use the fscore to measure the quality of obtained multi-matchings, which is defined as the fraction f = 2·p·rp+r , where p and r denote the precision and recall, respectively. All experiments are run on a Macbook Pro (2.8 GHz quad core i7, 16 GB RAM), where for = 10−5 we use f(Ut)/f(Ut+1) ≥ 1− as convergence criterion in Algorithm 1, and a step size of αt = ‖h(Ut, I)− hT (Ut, I)‖−1∞ in (7).
Real data. In this experiment we use the CMU house image sequence [1] comprising 111 frames within the experimental protocol of [34]. We generate a sequence of permutation synchronisation problem instances with a gradually increasing number of objects k by extracting respective pairwise matchings for k objects. To this end, we vary k from 20 to 111 and sample the pairwise matchings evenly-spaced from the 111× 111 pairwise matchings. Quantitative results in terms of the fscore, the objective value of Problem (9), and the runtimes are shown in Fig. 1, in which the individual problem instances vary along the horizontal axis. We can see that our proposed method dominates other approaches in terms of the fscore and the objective value (the reported objective values are divided by k2 to normalise the scale in problems of different sizes), while being among the fastest. Note that we do not report the objective value for MATCHEIG and MATCHALS, since they do not lead to cycle-consistent matchings so that the obtained solution does not lie in the feasible set of Problem (9). Qualitative results of the matching between one pair of images for k = 111 are shown in Fig. 2. As expected, our proposed approach clearly outperforms the SPECTRAL baseline, since our method is guaranteed to converge to the same same subspace that is spanned by the spectral solution, while at the same time providing a sparser and less negative solution that is thereby closer to the feasible set U. Synthetic data. We reproduce the procedure described in [8] for generating synthetic instances for the synchronisation of partial permutations. Four different parameters are considered for generating a problem instance: the universe size d, the number of objects k that are to be matched, the observation rate ρ, and the error rate σ (see [8] for details). One of these parameters varies in each experimental setting, while the others are kept fixed. Each individual experiment is repeated 5 times with different random seeds. In Fig. 3 we compare the performance of MATCHEIG [31], MATCHALS [51], SPECTRAL [34], NMFSYNC [8] and OURS. The first row shows the fscore, and the second row the respective runtimes. Note that we did not run MATCHALS on the larger instances since the runtime is prohibitively long. The methods MATCHEIG and MATCHALS do not guarantee cycle consistency, and are thus shown as dashed lines. It can be seen that in most settings our method obtains superior performance in terms of the fscore, while being almost as fast as the most efficient methods that are also based on a spectral relaxation (MATCHEIG and SPECTRAL).

heading: 5 Discussion and Limitations
text: Our algorithm has several favourable properties: it finds a globally optimal solution to Problem (1), it is computationally efficient, it has the same convergence behaviour as the Orthogonal Iteration algorithm [21], and promotes (approximately) sparse solutions. Naturally, since our solution is globally optimal (with respect to Problem (1)) and thus converges to im(Vd), the amount of sparsity that we can achieve is limited by the sparsest orthogonal basis that spans im(Vd). In turn, we rather interpret sparsity in some looser sense, meaning that there are few elements that are large, whereas most elements are close to zero. Furthermore, since the secondary sparsity-promoting problem (with objective in (4)) is generally non-convex, we cannot guarantee that we attain its global optimum.
Since permutation synchronisation was our primary motivation it is the main focus of this paper. In the case of permutations, the assumption that the eigenvalues corresponding to the Vd are equal (as we explain in Sec. 3.2) is reasonable, since this must hold for cycle-consistent bijective matchings. However, broadening the scope of our approach to other problems in which this assumption may
not be valid will require further theoretical analysis and possibly a different strategy for choosing the matrix Z(Ut). Studying the universality of the proposed method, as well as analysing different
convergence criteria and different matrix functions Z(Ut) are open problems that we leave for future work.
The proposed method comprises a variation of a well-known iterative procedure for computing the dominant subspace of a matrix. The key component in this procedure is the QR-factorisation, which is differentiable almost everywhere and is therefore well-suited to be applied within a differentiable programming context (e.g. for end-to-end training of neural networks).

heading: 6 Conclusion
text: We propose an efficient algorithm to find a matrix that forms an orthogonal basis of the dominant subspace of a given matrix W while additionally promoting sparsity and non-negativity. Our procedure inherits favourable properties from the Orthogonal Iteration algorithm, namely it is simple to implement, converges almost everywhere, and is computationally efficient. Moreover, our method is designed to generate sparser solutions compared to the Orthogonal Iteration algorithm. This is achieved by rotating our matrix of interest in each iteration by a matrix corresponding to the gradient of a secondary objective that promotes sparsity (and optionally non-negativity).
The considered problem setting is relevant for various spectral formulations of difficult combinatorial problems, such as they occur in applications like multi-matching, permutation synchronisation or assignment problems. Here, the combination of orthogonality, sparsity and non-negativity is desirable, since these are properties that characterise binary matrices such as permutation matrices. Experimentally we show that our proposed method outperforms existing methods in the context of partial permutation synchronisation, while at the same time having favourable a runtime.
Broader Impact
The key contribution of this paper is an effective and efficient optimisation algorithm for addressing sparse optimisation problems over the Stiefel manifold. Given the fundamental nature of our contribution, we do not see any direct ethical concerns or negative societal impacts related to our work.
Overall, there are numerous opportunities to use the proposed method for various types of multimatching problems over networks and graphs. The problem addressed is highly relevant within the fields of machine learning (e.g. for data canonicalisation to faciliate efficient learning with nonEuclidean data), computer vision (e.g. for 3D reconstruction in structure from motion, or image alignment), computer graphics (e.g. for bringing 3D shapes into correspondence) and other related areas.
The true power of permutation synchronisation methods appears when the number of objects to synchronise is large. With large quantities of data becoming increasingly available, the introduction of our scalable optimisation procedure that can account for orthogonality – while promoting sparsity and non-negativity – is an important contribution to the synchronisation community on the one hand. Moreover, on the other hand, it has the potential to have an impact on more general nonconvex optimisation problems, particularly in the context of spectral relaxations of difficult and large combinatorial problems.

heading: Acknowledgement
text: JT was supported by the Swedish Research Council (2019-04769). We thank Marvin Eisenberger for helpful feedback.

Review summaries:
In this paper the authors propose a modification of a well-known algorithm, in order to maximize a quadratic function over the Stiefel manifold, promoting sparsity at the same time. The base method is not a gradient-retraction over a manifold algorithm (although some tools are used in the middle), but a QR type of method, specifically the Orthogonal Iteration algorithm. The modification is the addition of a right multiplicative matrix, which is selected in order to promote sparsity. This is done for a maximization problem, and the $l_p$ norm (with $p>2$) is used to that end. The method is tested for the permutation synchronization problem.  

The paper is concerned with pursuing sparse optimal solutions to optimization problems of the form max Trace(U^T W U) [*] where U is an orthonormal matrix, that is U lies in the Stiefel manifold. This problem is meaningful since it typically arises as a relaxation over the Stiefel manifold of combinatorial problems such as spectral clustering and permutation synchronization. It this latter problem that occupies the application aspect of the paper. A global (not necessarily sparse) solution to [*] can be obtained in closed form from the eigendecomposition of W. Algorithmically, this can be computed by the power iteration algorithm. The paper proposes a variation of the power iteration algorithm (Algorithm 1), which includes a sparsity promoting step, forcing the basis matrix U to be close to a power of U. Permutation synchronization experiments with synthetic data and on the CMU "house" sequence for finding correspondences in two views,  show that this approach gives better permutation synchronization than several other existing alternatives. 

This paper proposes a novel algorithm to solve the permutation synchronization problem which can be formulated as minimizing a quadratic function over the Stiefel manifold. The algorithm is modified from the Orthogonal Iteration method with sparsity promotion. The proposed algorithm is shown to perform better than existing methods in terms of both efficiency and accuracy. Moreover, the limitation  of the work  is also discussed.