Sparse Quadratic Optimisation over the Stiefel Manifold with Application to Permutation Synchronisation
Authors: Florian Bernard
Daniel Cremers
Anders Johan Thunberg
Venue: NeurIPS 2021 Poster (2021)
Keywords: Stiefel manifold
quadratic optimisation
permutation synchronisation
sparsity
multi-matching
correspondence problems
manifold optimisation
QR decomposition
orthogonal iteration algorithm

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.

Sparse Quadratic Optimisation over the Stiefel Manifold with Application to Permutation Synchronisation
Florian Bernard TU Munich, University of Bonn
Daniel Cremers TU Munich
Johan Thunberg Halmstad University
Abstract
We address the non-convex optimisation problem of finding sparse matrix the Stiefel manifold (matrices with mutually orthogonal columns of unit length) that maximises (or minimises) a quadratic objective function. Optimisation problems the Stiefel manifold occur for example spectral relaxations of various combinato- rial problems, such graph matching, clustering, or permutation synchronisation. Although sparsity is a desirable property in such settings, is mostly neglected spectral formulations since existing solvers, e.g. based on eigenvalue decom- position, are unable to account for sparsity while at the same time maintaining global optimality guarantees. We fill this gap and propose simple yet effective sparsity- prommttng modification of the Orthogonal Iteration algorithm for finding the dominant eigenspace of matrix. By doing so, we can guarantee that our method finds 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 3D shapes in computer vision, computer graphics, and beyond. We demonstrate that the proposed approach outperforms previous methods this domain.
Introduction
We interested in optimisation problems of the form
arg f(U) with f(U) tr(UutTUU WU) and St(m. = (X E madd x La}. (I)
UESt(m,d)
where Ia is the identity matrix dimension d, W E Rmxm and the set St(m, d) denotes the Stiefel manifold (m d). Throughout the paper, w.l.o.g. consider maxinisation formulation and we assume that W is symmetric and positive semidefinite matrix (see Lemma in Sec. 2
Despite the non-convexity of Problem (I). problems of this form can be solved global optimality based on the eigenvalue decomposition, i.e. by setting U = Va. where Va E R denotes an orthogonal basis of eigenvectors corresponding to the largest eigenvalues of W, see Lemma The fact that we can efficiently find global optima of Problem (I) makes it popular relaxation formulation of various difficult combinatorial problems. This includes spectral relaxations of the NP-hard quadratic assignment problem (QAP) (35), spectral clustering
43
26 I7
[3
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 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 (NcurIPS 2021).

---

matrix. Often such combinatorial sets can be characterised by Stiefel matrices that are sparse for example, the set of (signed) d'dd permutation matrices can be characterised by matrices St(d. d) that have exactly d non-7ero elements, whereas any other element in St d) that has more than d non-zero elements is not signed permutation matrix.
The most common approach to obtain sparse solutions in an optimisation problem is to integrate an explicit pparsty-promotiig regulariser. However, a major hurdle when considering quadratic optimisation problems over the Stiefel manifold is that the incorporation of sparsity-roomoting terms not compatible with solutions based on eigenvalue decomposition, so that desirable global optimality guarantees are generally no longer maintained.
Instead, depart from the common path of integrating explicit regularisers and instead exploit the aathogonal-invariinn present in Problem (I to assure sparsity. To be more specific, for any orthogonal matrix Q O(d) (t,d, d) holds that f(U) =
tr(U"
tr((U))
( QU))) = f(UQ). Hence, if U E St(m, d) is a solution to Problem ① so is UQ E St(m, d) for any Q E (d). The subspace im(UQ) that is spanned by UQ for given U St(m, d)
(and arbitrary Q E ((d) is equal to the subspace im(U) Motivated by this observation we utilise orthogonal-invariance in order find a solution U E {U E 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 E ((d), which allows to rotate a given solution U E St d) sparser representation U UQ E St d) most notably, while remaining a globally optimal solution to Problem
main contributions follows: (i) For the first time propose an algorithm that exploits the aathogonal-invarinnee 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 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 global optimiser of Problem ① (iv) We experimentally confirm the efficacy of our approach in the context of the permutation synchronisation problem.
Main contributions. We summarise
2 Preliminaries & Related Work
In this section we clarify our assumptions, introduce additional preliminaries, and provide reference to related work. Let X Ag..... .,, Am be the eigenvalues of W mmm ordered decreasingly. We impose the following assumption on W:
Assumption 1 (Separated eigenspace). We assume that Xa Ad+1-
Throughout the paper we also assume that W is symmetric and positive semidefinite (p.s.d.), which, however, is not restriction as the following straightforward result indicates:
Lemma (Symmetry and positive semidefiniteness). In Problem T if 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 E [mmxm is defined as the subspace im(Va) C R", where Va E 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 (U,} converges to the dominant subspace of W loo IIVaVT U.] O.
There exists close relation between the dominant subspace of W and solutions of Problem |

---

Lemma 2 (Solution to Problem (T). Problem ① is solved for any matrix U that forms an orthog- onal basis for the d-dimensional dominant subspace of W. i.e. U {U St(m, d) iim(U) im(Va)).
Proof. See Appendix.
Basic algorithms for computing eigenvectors. In order to find the dominant subspace of W we can algorithms for finding eigenvectors. The Power method [21] is 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 given initial vo St,m, 1) based on the update 1t+11 t W ot / In order to find the d most dominant eigenvectors. one can consider the Orthogonal Iteration algorithm [2I]. which generalises the Power method to the case d computing Vi+1 R++1 + WV, based on the (thin) QR OR-decomposition of WV, where Vi E St(m, d) and R2 E Rdxd is upper triangular. For f. - oo the sequence (V)) converges (under mild conditions on the intial Vi) to the dominant invariant subspace as long as Aal |14+1l, see Thm. 8.2.2 in (21)
The algorithm proceeds by repeatedly
Sparse Stiefel optimisation. There are numerous approaches for addressing general optimisation problems over the Stiefel manifold, including generic manifold optimisation techniques (e.g. B (4), or Stiefel-specific approaches (e.g. [0 [46)). In the following we will focus on works that consider sparse optimisation over the Stiefel manifold that most relevant
approach. In order regularisers can be utilised, for example via the
to promote sparse solutions, minimisation of the (non-convex) l,- 'norm' for p 1, or the (convex) l nn-nrrm [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 non-convex smooth and a convex non-smooth function is optimised via proximal gradient method. In (27]. the authors consider the optimisation of sum of non-smooth weakly convex functions over the Stiefel manifold using Riemannian subgradient-wye terms are considered [40] 29] 38 [5].
methods. Yet, in practice often
differentiable surrogates of non-smooth pparsity-promotinn Instead of minimising Ep- "-"'norms' with S p S 1, on the Stiefel manifold one may instead choose the maximisation of bp -o-mms with p > 2, such as the fy-norm [7), the fn-norm [49]. Further motivation for sparsity promoting higher-ordee 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 I8]
including sparse optimisation via -nnor 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-invariaa of the objective f(U)= f(UQ) for any Q E (dd) in Problem [ finding the dominant subspace W could also be posed as an optimisation problem over the Grassmannian only interested in identifying this subspace (which can for example be
manifold. However, Iterations algorithm (21), or via Grassmann Grassmann-Ralleigh Quotient Iterations [4]), done via the Orthogonal but we want to find a specific choice of coordinates for which the representation of the subspace is sparse.
3 Proposed Method for Sparse Quadratic Optimisation over the Stiefel
Input: W E Rm Uo E Rmdd
Output: U*
Initialise: f. 1-+00
repeat
Ue41Be+1 +- WU,Z(U:) // unique QR-decomposition
until convergence
U --UU++
Algorithm 1: Overview of our proposed algorithm. In terms of convergence properties our algorithm is equivalent the Orthogonal Iteration algorithm and thus produces a U in the dominant subspace of W. However, our modification introduces the matrix Z(U)) with the purpose of promoting sparsity of the solution (see Sec. 3.2 how choose (U.)).
3

---

In this section we introduce our Algorithm 1 that can be seen a modification of the popular Orthogonal Iteration algorithm. The main difference is that we iteratively weigh the matrix U, e St(m, d) by matrix ((U)) Rdxd, where the matrix function Z maps onto the set of full rank matrices. The purpose of the matrix Z(U+) (see Sec. 3.2 for our specific choice) 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 global maximiser of Problem (I). which, in some relaxed sense, is also sparse. However, we cannot in general augment the objective function in Problem ( with an explicit sparsity regulariser and expect that the respective solution is still maximiser of the original objective. Instead we steer the solution towards being more sparse by weighing our matrix of interest U, with the matrix (U,)).
3.1 Convergence
We start by ensuring that, under the assumption that (U,),) full rank, the sequence {Ut} generated by Algorithm [ 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: Ve+1 Re+1 <- WV,XD
Algorithm B: Ve+1Be+1 - WV.
where X₁ E iadxd 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 R&+1 and R441 are upper triangular with positive diagonal. Now, if ViQu+ E St(m, d), where Qi E (d), then Vy is equal to Vy up to rotation from the right for all t' d'>t t.
PetiQens Rest where Qc+1 and Re+1 are the two matrices in the QR-factorisation of the matrix ( This means that V4+1 = V4+1Q++1 This is due to the uniqueness of the QR-factoris,tion, so that also Re+1 = R Now the results follows readily by using induction.
Proof. It holds that Ro+1 WVXX, = WV.Q.X, = VRe
Thus, if (U,) is full rank for each L, we can identify our Algorithm 1 with Algorithm A in Lemma 3 and conclude that the sequence (U,) converges the dominant subspace of W as (-> oc. This because Algorithm B in Lemma 3 corresponds to the Orthogonal Iteration algorithm, which is known to converge the dominant subspace of W (2I).
We continue by investigating the behaviour at the limit where the columns of U2 are in the dominant subspace of W. Let us consider the update assume that the d largest eigenvalues of W are equal, which is for example the case for synchronisation problems assuming cycle consistency (cf. Sec.(4).
< WU, of Algorithm [ and further
Lemma 4. Assume that the d largest eigenvalues of the p.s.d. matrix W are all equal to ond and strictly larger than the other eigenvalues, and assume that the columns of U, E St(m, d) are contained in the dominant subspace of W. Provided ZUU) is full rank, it holds that
WU,Z(U:) = U,Z(U,) = U2+ R++1 U.Qu+1
where Ut+1 Ret1 is the unique QR-factorisation of U,Z(U:) and Q++1 Ret1 is the unique QR- factorisation of 7UU.).
Proof. Since all the d largest eigenvalues of W are equal to one and the columns of U, are in the dominant subspace it holds that WU, U,U,U, Us. which explains the first equality. While the second equality follows by definition, the third equality remains to be proven. Let Q++1 R441 be the unique QR-decomposition of ZU,). We want to show that U U,Q041 and R,+1 = R From the QR-decomposition of (U)) it follows that U.Z(U:) UtQc+1R44 while from the
If they are all equal to value A * 1, we can w.l.o.g. consider X W in place of W, since the optimiser of Problem ( is invariant to scaling W.

---

decompossiiionn of (UZ(U;) it follows that (U)( = R441. The uniqueness of the QR- decomposition implies U2Q++1 = -U++1 and Ro+1 R++1-
L
Hence, under the assumptions Lemma 4 the columns of Up span the dominant subspace W for any l' 2 I. The update in Algorithm reduces to the form
Q4+1 Z(U,),
(2)
UL41 t UQQ++1,
(3)
i.e. updating U, to obtain U441 simplifies to multiplication of U/ with the Q-matrix from the QR- factorisation of (U)).
3.2 Choosing Z(U,) to Promote Sparsity
We now turn our attention to the choice of (UU) in Algorithm [ We know that, as long as Z(U)) is full rank in each iteration the columns of Ug converge to the dominant subspace of W. In addition to our objective function f in Problem ( we now introduce the secondary objective
9(U) Σ ( (U. ()),
(4)
where for p being positive integer the notation U.A means raise each element in U to the power For U E St(m, d) and p larger than 2, the maximisation of promotes sparsity in relaxed sense. To be specific, if p is odd, few larger elements (with value close to lead 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 11 lead to larger value of compared to many smaller elements close to O.
Let us consider the maximisation of
q(U)) (r(h(U,
with respect to Q, where h(U,Q) means that we want to rotate U by Q in such a way that the secondary objective 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 in ④
Q 02, and U E St d). This
The extrinsic gradient with respect to Q (for Q relaxed to be in Rdxd) of g(UQ) at Q I is p-h(U, I). Thus the (manifold) gradient of g(U Q) at Q = is given by the projection of the extrinsic gradient at Q I onto the tangent space, which reads
Vog(UQ))e==== = -
(6)
For small step size a, we can perform a first-order approximation of gradient ascent if we choose Z = I oxp(h (U, - 1tt (U, 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 retraction 3 and serves as a first-order approximation of the exponential map. With that, the updates become
Q4t1Re41 + I+ (HHaU hT (U.,
U441 +- U.Q41
(8)
This choice of the matrix Z(U:) ensures that it is full rank, since adding the identity matrix and a skem-symmetric matrix results in a full rank matrix. Hence, for sufficiently small step size a, iteratively move in an ascent direction on Od by utilising the QR-retraction.
Application to Permutation Synchronisation
4.1 Permutation Synchronisation in Nutshell
Permutation synchronisation is procedure to improve matchings between multiple objects 23 and related concepts have been utilised to address diverse tasks, such as multi-alignment 7 5
5

---

251. suatt-shape matching [22 24 20], multi-image matching [5 42 9 13 12], or multi -graph matching 48
among many others. Permutation synchronisation refers 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 typical application where the points are feature descriptors or key-points, and the objects are images, as shown in Fig. 2
Let & denote the number of objects, where each object i contains m, points. For 1, being p- dimensional vector of all ones, and vector inequalities being understood in an element-wise sense, let P.y E P Fmam: == {X {0, 1,m.xxyy X1m; S 1m; 1 X S 1 } be the partial permutation matrix that represents the correspondence between the m4 points in object i and the my points in object j. In the case of bijective permutations, the set of pairwise permutations P = [Pi])s-== is
said to be cycle-consistent if for all i, j, it holds that PiPPs = Paj.
We define the set of partial permutation matrices with full row-rank as Pm.d:= := E Fm.d' x11 = ]. 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-univers matchings U == (P; E F Fm.d).-1 such that for all i.j we can write P.j = P.P.T (see (23 34 for details). The object-to- universe characterisation of cycle consistency also valid for the case of non-bijective (i.e. partial) permutations (see [42 g for details).
Given the noisy (i.e. not cycle-consistent) set of pairwise permutations P =(Ps) ig permutation synchronisation can be phrased as the optimisation problem
Σ 1.
[urPPTPPPPP))
arg argmax tr(P"WP), PEU
arg max (P.EPmmdl
where for m := E, m. we define the set U = P mid X... x F C Rmxdd the (mxd) dimeensimaal block matrix P=PP _PT..... ..... PT and the block matrix W = [Pojli E Rm in order to allow for a compact matrix representation of the objective.
4.2 Proposed Permutation Synchronisation Approach
The core idea of existing spectral permutation synchronisation approaches [34 39 31] is to replace the feasible set in Problem (9) with the Stiefel manifold St(m, d). so that we obtain an instance of Problem 11. We utilise Algorithm| | to obtain (globally optimal) solution U E 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(U'T = Σ Σj= Uß With that, in addition to U being an orthogonal matrix, it contains few large elements that are close to I and many smaller elements that are close to As such, we can readily project the matrix U onto the set U in terms of Euclidean projection. The Euclidean projection is given by projecting each of the h blocks of U = [UT ................... ... U individually onto the set of partial permutations, i.e. P, projF mo. (U?)
arg max af, Uf)) see e.g. 8 which amounts to a (partial) linear assignment problem (32 that we solve based on the efficient implementation in [00]
4.3 Experimental Results
We experimentally compare our proposed approach with various methods for permutation synchroni- sation 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 MATCHAL L.S do not guarantee cycle consistency (but instead aim to improve the initial matchings), in 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-r where p and r denote the precision and recall, respectively. All experiments are run on a Macbook Pro (2.8 GHz quad core 17, 16 GB RAM), where for = 10-5 we use ((U,) (U,)//(U+++ 2 e as convergence criterion in Algorithm| and a step size of a; IIh(U,,I) h (U,, I)
p-r
in

---

MatchEig MatchALS
Spectral
NrifSync
Ours
objective
runtime [s]
fscore
0.9
25
10
0.8
20
40
60
80
100
20
60
80
100
20
40
60
80
100
k
k
Figure Quantitative results on the CMU house sequence in terms of the fscore (t), objective value of Problem ( (1), and runtime (1). The individual instances of permutation synchronisation problems vary along the horizontal axis. Methods that do not guarantee cycle consistency are shown as dashed lines.
Real data. In this experiment we use the CMU house image sequence [ comprising 111 frames within the experimental protocol of [34] We generate a sequence of permutation synchronisation problem instances with gradually increasing number of objects & by extracting respective pairwise matchings for & objects. To this end, we vary & 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. [ 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 6² 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 sparser and less negative solution that 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 & that are to be matched, the observation rate and the error rate a (see [8] for details). One of these parameters varies in each experimental setting, while the others are kept fixed. Each individual experiment repeated 5 times with different random seeds. In Fig. 3 we compare the performance of MATCHEIG B1], MATCHALS [5I], SPECTRAL 34,, NMESYNC [ and OURS. The first row shows the ficcore and the second row the respective runtimes. Note that we did not run MATCHALS on the larger instances since the runtime 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 spectral relaxation (MATCHÉIG and SPECTRAL).
Discussion and Limitations
Our algorithm has several favourable properties: it finds a globally optimal solution to Problem [. 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 ( and thus converges to im (Va), the amount of sparsity that we can achieve limited by the sparsest orthogonal basis that spans (V2)... 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 is generally non-convex, we cannot guarantee that we attain its global optimum.
Since permutation synchronisation was our primary motivation is the main focus of this paper. In the case of permutations, the assumption that the eigenvalues corresponding to the Va2 are equal (as we explain Sec. 3.2p 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
7

---

INPUT
MATCHÉIG
MATCHALS
SPECTRAL.
NMFSYNC
OURS
Figure 2: Comparison of matchings between the first and last image the CMU house sequence obtained by several methods. The colour of the dots indicates the ground truth correspondence, and the lines show the obtained matchings (green: correct, red: wrong). Overall, our approach obtains the best matchings, see also Fig. 1
MatchEig
MatchALS
Spectral
NmfSync
p-
50, P
50,
50,
Iscore
iscore
fscore
0.8
0.8
0.5
0.6
runtinie [8]
runtime [s]
runtime [8]
500
runtime
100
50
20
80
0.3
0.5
0.8
0.4
0.6
0.8
200
300
400
500
Figure 3: Experimental comparison of permutation synchronisation methods on synthetic data. Each column shows different varying parameter. The first row shows the fscore (1) and the second row the runtime (L). Methods that do not guarantee cycle consistency are shown with dashed lines. Our method guarantees cycle consistency and leads to higher fscores in most settings.
not be valid will require further theoretical analysis and possibly a different strategy for choosing the matrix Z(U)). Studying the universality of the proposed method, as well as analysing different

---

convergence criteria and different matrix functions (()) 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 matrix. The key component in this procedure is the OR-factorisation, which is differentiable almost everywhere and is therefore well-suited to be applied within differentiable programming context (e.g. for end-to-end training of neural networks).
Conclusion
propose an efficient algorithm to find a matrix that forms an orthogonal basis of the dominant subspace of given matrix W while additionally promoting sparsity and non-negativity. Our procedure inherits favourable properties from the Orthogonal Iteration algorithm, namely is simple to implement, converges almost everywhere, and is computationally efficient. Morcover, method is designed to generate sparser solutions compared to the Orthogonal Iteration algorithm. This is achieved by rotating matrix of interest in each iteration by 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 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 not see any direct ethical concerns or negative societal impacts related to our work.
Overall, there numerous opportunities to the proposed method for various types of multi- matching problems over networks and graphs. The problem addressed is highly relevant within the fields of machine learning (e.g. for data canonicalisation faciliate efficient learning with non- Euclidean 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 account for orthogonality while promoting sparsity and non-negativity is an important contribution to the synchronisation community on the one hand. Moreover, the other hand, it has the potential have impact more general non- convex optimisation problems, particularly in the context of spectral relaxations of difficult and large combinatorial problems.
Acknowledgement
JT was supported by the Swedish Research Council (2019-04769). We thank Marvin Eisenberger for helpful feedback.
References
[1] CMU/VASC image database. http://www.. cs. http://wwwccccmee....... http::/p/wwoooooo.ccccc.////////
html/mmotii html/motion//
[2] P-A. Absil, C. G. Baker, and K. A. Gallivan. Trust-region methods riemannian manifolds. Foundations of Computational Mathematics, (33:3303-300, 2007.
[3] P.A. Absil, R. Mahony, and R. Sepulchre. Optimization algorithms on matrix manifolds. Princeton University Press, 2009.

---

[4] P.-.. Absil, R. Mahony, R. Sepulchre, and P. Van Dooren, A grasmmnnnraylleghh quotient iteration for computing invariant subspaces. SIAM review, 44(1):57-73, 2002.
[5] F. Arrigoni, B. Rossi, and A. Fusiello. Spectral synchronization of multiple views in se (3). SIAM Journal Imaging Sciences, 9(4):1963-199)
2016.
[6] F. Bernard, C. Theobalt, and M. Moeller. DS Tighter Lifting-Free Convex Relaxations for Quadratic
Matching Problems. In CVPR, 2018.
[7] F. Bernard, J. Thunberg, P. Gemmar, F. Hertel, A. Husch, and J. Goncalves. solution for multi-alignment
by transformation synchromisation. In CVPR, 2015.
[8] F. Bernard, Thunberg, J. Goncalves, and C. Theobalt. Synchronisation partial multi- matchings via non noganneatvve factorisations. Pattern Recognition, 92:146-1555 2019.
[9] F. Bernard, Thunberg, Swoboda, and Theobalt. Hippi: Higher-order projected power iterations for scalable multi-matching. CVPR, 2019.
[10] F. Bernard, N. Vlassis, P. Gemmar, A. Husch, J. Thunberg. J. Goncalves, and F. Hertel. Fast correspondences statistical shape models of brain structures. In Medical Imaging 2016: Image Processing, volume 9784, page 97840R. International Society for Optics and Photonics, 2016.
[11] D. Bertsekas. Network Optimization: Continuous and Discrete Models. Athena Scientific, 1998.
[12] T. Birdal, V. Golyanik, C. Theobalt, and L. J. Guibas. Quantum permutation syachronization. In CVPR,
2021.
[13] T. Birdal and U. Simsekli. Probabilistic permutation synchronization using the riemannian structure of the birkhoff polytope. In CVPR, 2019.
[14] N. Boumal, B. Mishra, P.-.. Absil, and R. Sepulchre. Manopt, a Matlab toolbox for optimization on manifolds. Journal of Machine Leurning Research, 15542): 2014.
[15] A. Breloy, S. Kumar, Sun, and D. P. Palomar. Maajorizationnminimizaa application to robust sparse pea. IEEF Transactions on Signal Processing, 69:1507-1520 2021. Zhang. Proximal gradient method for nonsmooth optimization over the stiefel manifold. SIAM Journal on Optimization, 30(1):210-25 2020.
on the stiefel manifold with
[16] S. Chen, S. Ma, Man- Cho So,
[17] T. Cour, P. Srinivasan, and J. Shi. Balanced graph matching. NIPS, 2006.
[18] L. De Lathauwer, L. Hoegaerts, and J. Vandewalle. A grassmant-rayleighh quotient iteration for dimen- sionality reduction in ica. In International Conference on Independent Component Analysis and Signal Separation, pages 335 342. Springer, 2004.
[19] A. Edelman, T. A. Arias, and S. T. Smith. The geometry of algorithms with orthogonality constraints.
SIAM journal on Matrix Analysis and Applications, 20(2):303-333, 1998.
[20] M. Gao, Z. Lähner, J. Thunberg, D. Cremers, and F. Bernard. Isometric multi-shape matching. In CVPR,
2021.
[21] G. H. Golub and C. F. Van Loan. Matrix computations, volume 4. JHU press, 2013.
[22] Q. Huang, Z. Liang. Wang, S. Zuo, and C. Bajaj. Tensor maps for synchronizing heterogeneous shape
collections. ACM Transactions on Gruphics (TOG), 38(4):1-18, 2019.
[23] Q.-X. Huang and I.. Guibas. Consistent shape maps via semidefinite programming. In Symposium on Geometry Processing, 2013.
[24] R. Huang. J. Ren, P. Wonka, and M. Ovsjanikov. Consistent zoomout: Efficient spectral map synchroniza-
tion. Computer Graphics Forum, volume 39, pages 265-278 Wiley Online Library, 2020.
[25] X. Huang, 7. Liang, X. 7hou, Y. Xie, I.. Guibas, and Q. Huang. Learning transformation synchronization.
In CVPR, 2019.
[26] M. Leordeanu and M. Hebert. A Spectral Technique for Correspondence Problems Using Pairwise Constraints. In ICCV, 2005.
[27] X. Li, S. Chen, Z. Deng, Q. Qu, Z. Zhu, and A. Man-Cho So. Weakly convex optimization over stiefel
manifold using riemannian subgradient-type methods. SIAM Journal on Optimization, 3053)1100-1164, 2021.
[28] Y. Li and Y. Bresler. Global geometry of multichannel sparse blind deconvolution on the sphere. In NewrIPS, pages 1132- 1143,114,,10 2018.
[29] C. Lu, S. Yan, and Z. Lin. Convex sparse spectral clustering: Single-view to multi-view. IEEE Transactions Image Processing, 2833)2883,,,,,,, 2016.
[30] H. Manton. Optimization algorithms exploiting unitary constraints. IEEE Transactions on Signal
Processing, 50(3):635-655 2002.
[31] E. Maset, Arrigoni, and A. Fusiello. Practical and Efficient Multi- Vietw Matching. In ICCV. 2017.
[32] J. Munkres. Algorithmss for the Assignment and Transportation Problems. Journal of the Society for Industrial and Applied Mathematics, 5(1):32-38, Mar. 1957.
[33] A. Ng, M. Jordan, and Y. Weiss. On spectral clustering: Analysis and an algorithm. NIPS, 14:849-856.
2001.
[34] D. Pachauri, R. Kondor, and V. Singh. Solving the multi-way matching problem by permutation synchro- nization. In NIPS, 2013.
10

---

[35] P. M. Pardalos, F. Readl, and H. Wolkowicz. The Quadratic Assignment Problem A Survey and Recent Developments. DIMACS Series Discrete Mathematics, 1993.
[36] Q. Qu. J. Sun, and J. Wright. Finding a sparse vector in a subspace: Linear sparsity using alternating directions. IEEE Trunsactions on Information Theory, 62( 62(10))5855-5800 2016.
[37] Q. Qu. Zhai. X. Li, Zhang, and Z. Zhu, Geometric analysis of nonconvex optimization landscapes for overcomplete learning. In International Conference on Learning Representations, 2019.
[38] Q. Qu, Zhu, X. Li, M. C. Tsakiris, Wright, and R. Vidal. Finding the sparsest vectors in subspace: Theory, algorithms, and applications. arXiv preprint arXiv:2001,06%3 06970, 2020.
[39] Y. Shen, Q. Huang, N. Srcbro, and S. Sanghavi. Normalized spectral map synchronization. NIPS,
29:49255-4333, 2016.
[40] J. Song, P. Babu, and D. P. Palomar. Sparse generalized eigenvalue problem via smooth optimization.
IEEE Transactions on Signal Processing, 63(7):1627-1442 2015.
[41] P. Swohoda, D. Kainmüller, A. Mokarian, C. Theobalt, and Bernard. convex relaxation for multi-graph
matching. In CVPR, 2019.
[42] R. Tron, X. Zhou, C. Esteves, and K. Daniilidis. Fast multi-image matching via density-based clustering.
In CVPR, 2017.
[43] S. Umeyama. An eigendecomposstioon approach weighted graph matching problems. IEEE transactions on pattern analysis and machine intelligence, 10(5):695-703, 1988.
[44] U. Von Luxburg. A tutorial on spectral clustering. Statistics and computing, 17(4):395-416, 2007.
[45] Q. Wang, Gao, and H. Li. Grassmannian manifold optimization assisted sparse spectral clustering, In Proceedings of the IEFE Conference on Computer Vision and Pattern Recognition, pages 5258-5266,
2017.
[46] Z. Wen and W. Yin. A feasible method for optimization with orthogonality constraints. Mathematical Programming. 14377443,,,433, 2013.
[47] Y. Xue, Y. Shen, V. Lau, J. Zhang. and K. B. Letaief. Blind data detection in massive mimo via fo-morm maximization over the stiefel manifold. IEEE Transactions Wireless Communications, 2020.
[48] J. Yan, X.-. Yin, W. Lin, C. Deng, H. Zha, and X. Yang. A short survey of recent advances in graph matching. In Proceedings of the 2016 ACM on International Conference Multimedia Retrieval, pages
167-174, 2016.
[49] Y. Zhai, Z. Yang, Z. Liao, J. Wright, and Y. Ma. Complete dictionary learning via 14-norm maximization
over the orthogonal group. Journal of Machine Learning Research, 21(165): 2020.
[50] Y. Zhang. H.-W. Kuo, and J. Wright. Structured local optima in sparse blind deconvolution. IEEE Transactions on Information Theory, 66(1):4194452 2019.
[51] X. Zhou, M. Zhu, and K. Danillidis, Multi-image matching via fast alternating minimization. In ICCV, 2015.
[52] Z. Zhu, T. Ding, D. Robinson, M. Tsakiris, and R. Vidal. A lincarly convergent method for non-smooth non-convex optimization the grassmannian with applications robust subspace and dictionary learning. NeurlPS, 32:9442-9452, 2019.


Reviews:
[{"review_id": "TUOncEBifgA", "reviewer": "Reviewer_fYRj", "summary": "They address sparse quadratic optimization over the Stiefel manifold with application to permutation synchronization.", "questions": "", "limitations": "", "rating": 7, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Existing solvers, e.g. based on eigenvalue decomposition, are unable to account for sparsity while at the same time maintaining global optimality guarantees. They 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. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "10"}, "extra_text": {}}, {"review_id": "PAqqFMUXp5-", "reviewer": "Reviewer_C7cB", "summary": "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. ", "questions": "", "limitations": "", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The idea of forcing sparsity on U towards recovering a solution of the original problem (the one whose relaxation is [*]) is interesting and important. Moreover, the idea of modifying the power iteration method with the sparsity promoting step of (4) is interesting as well. \n\nOn the other hand, a weakness in the paper is that the authors do not treat the attribute \"sparse\" precisely.  In page 4, line 136, when they introduce their algorithm, they say that they are interested in an optimal solution of [*] which is sparse \"in some relaxed sense\". Only in the discussion section at page 7 they write \"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.\" I feel that this should have been clarified very early on. \n\nAnother weakness is that the authors do not provide any experimental analysis of Algorithm 1 alone.  \nI run the code of the authors on a very simple case, where W=diag(1,1,1,0.5,0.5,0.5), with d=3 and various values of p ranging from 2 to 20. Ideally, the algorithm should be converging to the matrix E=[e1, e2, e3] where e1,e2,e3 are the first three standard basis vectors of R^6. The algorithm converges very quickly (in two iterations or so) to a matrix that is far from E (in fact most elements in the matrix seem to have the same order of magnitude), but if one only selects the largest in magnitude element per column and sets the rest to zero, then one does get E M, where M is a diagonal matrix (this is more clear for larger values of p). More precisely, if one runs the linear assignment post-processing step that the authors propose, then we get back E for p=10, but a wrong result for p=2,3. This shows that both 1) the linear assignment step  and 2) the value of p, are crucial. I also did another experiment at the other extreme, where I chose W to be a generic 6x6 positive semi-definite matrix. With d=5 one expects 4+3+2+1 zeros (this is the highest sparsity of an orthonormal basis of a generic 5-dimensional subspace of R^6). With p=10, while again most elements have the same order of magnitude, there is indeed a clear cluster of 10 smaller elements. \n\nIn conclusion, Algorithm 1 is an interesting heuristic, but more investigation is needed.\nSince it is so simple, perhaps it is possible to actually provide a precise mathematical analysis. For this, the authors should define their approximate notion of sparsity rigorously and prove that if the d-th principal eigenspace of W admits an orthonormal basis of sparsity level \\alpha, then for p large enough, Algorithm 1 converges to an orthonormal basis with approximate sparsity \\alpha. \nI believe this would make a strong paper. \n\nThree more suggestions:\n1) I think it would be more interesting to apply permutation synchronization to a less explored and modern application, such as record linkage, instead of feature matching in computer vision, where there is a vast literature. \n2) Lemmas 1 and 2 are textbook level, there is no need to include the proofs (in fact the given proof of Lemma 2 is not complete). \n3) The notation U* \\in im(bar(U)) is not precise. \n\n**** Post-Rebuttal ****\n\nI thank the authors for their responses. I would like to keep my rating as is, in anticipation of a stronger version of this work. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "No issues as far as I can see.", "review_text": "", "extra_scores": {"ethical_concerns": "No issues as far as I can see.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "14"}, "extra_text": {}}, {"review_id": "NM3LpGptADM", "reviewer": "Reviewer_9Pkc", "summary": "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. ", "questions": "", "limitations": "", "rating": 6, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "1. The proposed algorithm has strong intuitions. And the presentation of this work is easy to follow.   To promote the sparsity, the objective function $g(UQ)$ is introduced. The authors point out that $g(UQ)$ is a relaxation of the sparsity and non-negativity.    Since $g(UQ)$ is nonconvex, the step (7) may converge to the local minimum points of (5). The related discussions are given in Appendix. It's better to mention this in the main paper. \n\n2. There is also a post-processiong step after algorithm 1  as stated in Section 4.2. It would be better to put them together in the Algorithm environment. \n\n3. The reviewer is curious about the numerical results of the case $p=4$ and minimizing the following function $\\sum_{ij} |(UQ)_{ij}|$. For the latter nonsmooth function, you may try  the Riemannian subgradient method [Li. et al 2021] with QR retraction.\n\nMinor:\n1. The proofs of Lemma 1 and 2 are trivial. I suggest putting them in the Appendix. \n2. Line 138: 'solution to (1)'\n\nReference:\n[Li. et al 2021] Xiao Li, Shixiang Chen, Zengde Deng, Qing Qu, Zhihui Zhu, Anthony Man-Cho So. Weakly Convex Optimization over Stiefel Manifold Using Riemannian Subgradient-Type Methods. SIAM Journal on Optimization (2021) 31(3):1605-1634.\n\n-----------------------------------------------------------------------\nPost rebuttal:\nThe authors sucessfully address my concerns.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Limitations:\nThe limitations have been discussed by the authors.\n\nSocietal Impact:\nNot apply.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "1u2Ghbv_qYF", "reviewer": "Reviewer_ZQ5Z", "summary": "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.\n\n", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The method is a clever modification of a very well-known algorithm. The novelty is the product with a matrix $Z(U_t)$, which doesn't affect the convergence of the original method, and at the same time it allows to incorporate a second objective function, which is chosen as the $l_p^p$ norm of the vectorized matrix.\nThe existing literature is correctly described, and the difference with respect to existing approaches is discussed.\n\n\nThe presented theoretical results seem to be correct. The proofs of the first part are simple, using basically linear algebra tools, but this is sufficient to prove the convergence and the bahaviour of the method in general.\n\nThe paper is correctly written, and it provides sufficient elements to reproduce the results. I wouldn't say is an \"easy to follow\" paper. I had to go back and forth several times in order to understand some parts. At one point I kept reading even though I was not completely understanding the goal of some points, with the hope to understand it later, which I think I did. I think this is in part due to the method itself, and not fault of the authors.\n\nI think that the contributions are significant, in the sense that the problem is relevant, this method seems to perform well for some difficult tasks, and the modifications are simple but clever. I'm not sure if those ideas could be used for other methods, not based on the QR methodology. Maybe they can.\n\nSome comments and minor comments:\nIn the experimental section, I would have liked to see an experiment showing the effect of the sparsity promoting addition. Maybe comparing the sparsity of resulting matrix from the original Orthogonal Iteration algorithm (which will be dense most likely), with sparsity resulting from your method.\n\nIn lines 112-113, the sentence might be read as if $l_p$ norms with $p<1$ are smooth (I know the observation is about the convexity).\n\nPost rebuttal: The authors addressed my concerns\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The limitations are correctly stated in Section 5, but specially the one regarding the equal eigenvalues is easy to overlook in the rest of the text.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The proposed method for sparse quadratic optimization over the Stiefel manifold and its analysis are interesting. However, the technical presentation and theoretical results have to be strengthened.
Decision: Accept (Poster)