--- title: "Non-smooth Optimization over Stiefel Manifolds with Applications to Dimensionality Reduction and Graph Clustering" authors: - "Fariba Zohrizadeh" - "Mohsen Kheirandishfard" - "Farhad Kamangar" - "Ramtin Madani" doi: "10.24963/ijcai.2019/183" url: "https://ijcai.org/proceedings/2019/183" pdf: "https://ijcai.org/proceedings/2019/0183.pdf" venue: "Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence" year: 2019 pages: "1319-1326" track: "Main track" keywords: - "Machine Learning: Explainable Machine Learning" - "Heuristic Search and Game Playing: Combinatorial Search and Optimisation" --- # Non-smooth Optimization over Stiefel Manifolds with Applications to Dimensionality Reduction and Graph Clustering Fariba Zohrizadeh, Mohsen Kheirandishfard, Farhad Kamangar, Ramtin Madani *Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence*, pages 1319-1326. DOI: [10.24963/ijcai.2019/183](https://doi.org/10.24963/ijcai.2019/183) ## Abstract This paper is concerned with the class of non-convex optimization problems with orthogonality constraints. We develop computationally efficient relaxations that transform non-convex orthogonality constrained problems into polynomial-time solvable surrogates. A novel penalization technique is used to enforce feasibility and derive certain conditions under which the constraints of the original non-convex problem are guaranteed to be satisfied. Moreover, we extend our approach to a feasibility-preserving sequential scheme that solves penalized relaxation to obtain near-globally optimal points. Experimental results on synthetic and real datasets demonstrate the effectiveness of the proposed approach on two practical applications in machine learning. ## Full text ## 1 Introduction Consider the following optimization problem where g0 : Rn×m → R is a convex piecewise linear function and ¯ fk : Rn×m → R is an arbitrary quadratic function of the form ¯ fk(P)≜〈Mk,PP⊤〉+〈Nk,P〉+qk, for every k ∈ {0, 1, . . . , p}, and {Mk ∈ Sn}pk=0 , {Nk ∈ Rn×m}pk=0 and {qk ∈ R}pk=0 are given. With no loss of generality, we assume that q0 = 0 and write g0 in the form of g0(P) = ‖α(P) + b‖1 , where b ∈ Rw is a given vector, α : Rn×m → Rw is a linear matrix function defined as α(Y ) ≜ ∑wi=1 〈Ai, Y 〉ei, the matrices {Ai ∈ Rn×m}wi=1 are given, and {ei ∈ Rw }wi=1 represent the standard basis for Rw . The formulation (1a) - (1c) encompasses a broad class of computationally-hard optimization problems with a variety of practical applications in discriminative dimensionality reduction [Bian and Tao, 2011], graph matching [Jiang et al., 2017], feature selection [Tang and Liu, 2012; Yang et al., 2011], compressed modes [Ozolin¸ˇs et al., 2013; Chen et al., 2016], among other areas of machine learning. The majority of methods in the literature are focused on a special case of (1a) - (1c) that involves the minimization of a convex and smooth objective function over non-convex sets of the form Sn,m ≜ {P ∈ Rn×m | P ⊤P = Im}, known as the Stiefel manifolds. There are various iterative local search algorithms which preserve the structure of Stiefel manifolds via geodesics steps [Edelman et al., 1998; Abrudan et al., 2008] or retractions [Absil et al., 2009; Wen and Yin, 2013]. Although these algorithms exhibit satisfactory performance in dealing with orthogonality constraints, they mostly restrict the objective function to the class of smooth functions and are not compatible with additional constraints [Gao et al., 2018]. To overcome these limitations, general algorithms are proposed that work with either smooth or non-smooth objective functions [Bian and Tao, 2011; Chen et al., 2016]. The paper [Bian and Tao, 2011] uses a family of semidefinite programming (SDP) problems to generate a converging sequence of points Stiefel manifolds. The paper [Chen et al., 2016] introduces an inner-outer iteration scheme for solving ℓ1-regularized optimization problems with orthogonality constraints based on the augmented Lagrangian method from [Fortin and Glowinski, 2000] and the proximal alternating minimization technique from [Attouch et al., 2013]. Moreover, a series of splitting techniques are proposed in [Ozolin¸ˇs et al., 2013; Lai and Osher, 2014; Kovnatsky et al., 2016] that can efficiently handle nonsmooth objective functions. They partition the problem into multiple sub-problems with analytical solutions and employ Bregman iterations [Yin et al., 2008] or its variants [Boyd et al., 2011] to obtain optimal solutions for orthogonalityconstrained problems. In the more recent paper [Zhu et al., 2017], an extended proximal alternating linearized minimization method is introduced to minimize convex functions subject to linear constraints and generalized orthogonality constraints. The success of related sequential frameworks and penalized relaxations for non-convex optimization is demonstrated in [Ibaraki and Tomizuka, 2001; Kheirandishfard et al., 2018b; Kheirandishfard et al., 2018a]. In [Ibaraki and Tomizuka, 2001], a sequential framework is introduced for solving bilinear matrix inequalities without theoretical guarantees. In [Kheirandishfard et al., 2018b; Kheirandishfard et al., 2018a], this approach is further investigated and theoretical results are offered through the notion of generalized Mangasarian-Fromovitz regularity condition. Another sequential SDP-based algorithm for pattern recognition is introduced in [Bian and Tao, 2011] that is not feasibility preserving. ## 1.1 Contributions Differentiated from the existing literature, we propose a computational approach with theoretical analysis for solving problems of the form (1a) - (1c), that guarantees the recovery of feasible points. The proposed approach generalizes the existing literature by including additional quadratic inequality constraints. The core of our approach is based on a novel and computationally efficient convex relaxation which transforms the non-convex problem (1a) - (1c) into a convex quadratically-constrained quadratic program (QCQP). To ensure that the solution of the relaxed problem is feasible for (1a) - (1c), we incorporate a penalty term into the objective function and derive certain conditions that guarantee the recovery of feasible points. Moreover, under certain conditions, we prove that by starting from any arbitrary initial point on a Stiefel manifold (not necessarily feasible), a sequence of penalized relaxations can be solved to find a feasible and near-optimal point. Unlike the existing algorithms, if mild assumptions are satisfied, the proposed sequential scheme is feasibility-preserving and improves the objective monotonically at every step. To corroborate the effectiveness of our method, we perform experiments on two practical applications with both synthetic and real datasets. The experimental results demonstrate that the proposed approach exhibits comparable results for both applications. ## 1.2 Notation Throughout this paper, the scalars, vectors, and matrices are shown by italic, bold lower-case and bold upper-case letters, respectively. The symbols Rn , Rn×m, Sn, and S+n denote the set of real n-dimensional vectors, real n × m matrices, real symmetric n × n matrices, and real positive semidefinite matrices, respectively. The symbols tr{.} and (.)⊤ are indicative of the trace and transpose operators, respectively. Given a vector a and a matrix A, the symbols ai and Aij, respectively, refer to the ith element of a and the (i, j)th element of A. The notation A ⪰ 0 states that A is symmetric positive semidefinite. Given matrices A and B of the same size, 〈A,B〉 ≜ tr{A⊤B} and A ◦ B, respectively, denote the Frobenius inner-product and the Hadamard product of A and B. The operator diag(.) gets a vector and forms a diagonal matrix with its input on the diagonal elements. The notation ‖.‖p refers to either matrix norm or vector norm depending on the context, ‖.‖F shows the Frobenius norm, and |.| indicates the absolute value or the cardinality of a set depending on the context. The symbol Im denotes the identity matrix of size m and the letter K is used as a shorthand for the set {1, . . . , p}. The symbol Sn,m as the set of real n × m matrices with orthonormal columns, i.e., Sn,m ≜ {P ∈ Rn×m | P ⊤P = Im}. The projection operator projSn,m : Rn×m → Sn,m is defined as projSn,m H = argmin{‖P -H‖F | P ∈ Sn,m}. ## 2 ProblemFormulation Optimization problems of the form (1a) - (1c) can be computationally challenging due to the non-convexities of the objective function and constraints. In order to derive convex relaxations, we first lift the problem into a higher dimensional space by introducing an auxiliary variable X ∈ Sn, accounting for the quadratic term PP ⊤. For every k ∈ {0} ∪ K, define fk : Rn×m × Sn → R as: Using the auxiliary variable X, the optimization problem (1a) - (1c) can be equivalently reformulated as with a convex objective function and convex linear inequality constraints (3b). The above formulation is still not convex due to the presence of the constraints (3c) and (3d) that capture all non-convexities of the problem. ## 2.1 Convex Relaxation In order to convexify the lifted problem (3a) - (3d), we relax the constraints (3c) and (3d) to where C ⊆ Sm and D ⊆ Sn are convex cones to be defined. In this work, we consider the common-practice semidefinite programming (SDP) relaxation and introduce a novel convex relaxation that transforms the problem (3a) - (3d) into a convex quadratically-constrained quadratic program (QCQP). ## Semidefinite Programming Relaxation This relaxation provides a powerful method for tackling nonconvex polynomial optimization problems [Boyd and Vandenberghe, 2004]. The SDP relaxation of the problem (3a) - (3d) can be derived by having C = S+m and D = S+n . Despite the effectiveness of this relaxation in providing high-quality solutions, its applicability is limited to the problems of moderate size due to the computational cost of imposing highdimensional conic constraints. ## Convex Quadratic Relaxation We propose a computationally efficient convex relaxation as an alternative to the SDP relaxation. In order to formulate the proposed relaxation for the problem (3a) - (3d), we need to have C = Vm and D = Vn, where for every positive integer o, set Vo ⊆ So is defined as follows Remark 1. It can be easily observed that if (C,D)=(Vm,Vn), the constraints (3c) and (3d) boil down to convex quadratic inequalities. Hence, the proposed relaxation reduces (3a) - (3d) to a convex QCQP. Notice that either of the aforementioned relaxations may fail to produce a feasible point for (1a) - (1c), because in general, an optimal solution to a convex relaxation does not necessarily satisfy the constraints (3c) and (3d). In what follows, we propose a penalization technique that guarantees the recovery of feasible points for (1a) - (1c) under certain conditions. ## 3 Penalization In this section, we show that by including a penalty term in the objective, one can obtain feasible points for the nonconvex problem (3a) - (3d). Given an arbitrary initial point ˇP ∈ Sn,m, that is not necessarily feasible, we transform the problem (3a) - (3d) into the following convex relaxation with revised objective function: where (C, D) ∈ {(S+m, S+n ), (Vm, Vn)}, and the fixed parameter µ > 0 sets a trade-off between the original objective function and the linear penalty term 〈P, ˇP〉. Remark 2. If an optimal solution ( ∗P , ∗X) of the problem (5a) - (5e) satisfies the constraints (3c) and (3d), then ∗P is feasible for (1a) - (1c). In the remainder of this section, certain conditions are introduced to guarantee that the penalized relaxation (5a) - (5e) produces feasible points for the non-convex problem (3a) - (3d). Definition 1. Define feasibility distance dF :Rn×m →R as where F denotes the feasible set of the problem (1a) - (1c). Definition 2. Define the singularity function s :Sn,m→R as: where ZP ≜ {D ∈ Rn×m | P⊤D = 0 ∧ ‖D‖F ≤ 1}. A point P ∈ Sn,m is said to satisfy the Mangasarian-Fromovitz constraint qualification (MFCQ) condition if it is feasible for the problem (1a) - (1c) and s(P ) > 0. Theorem 1. Define the constants and let ˇP ∈F be a feasible point for the problem (1a) - (1c) that satisfies the MFCQ condition. If then the penalized relaxation (5a) - (5e) has a unique optimal solution ( ∗P , ∗X), that satisfies (3c) and (3d). Moreover, ∗P is feasible for (1a) - (1c) and ¯ f0( ∗P )+g0( ∗P) ≤ ¯ f0( ˇP )+g0( ˇP ). Remark 3. For every point P ∈ Sm,n, it is straightforward to calculate s(P) by solving the following convex problem: Notice that β can be simply lower- and upper-bounded by any arbitrary member of the set Sm,n and the constant Ψ, respectively. This certifies the existence of a bounded µ that satisfies (9). In practice, there is no need to compute s(P) for fine-tuning parameter µ, since (9) offers a conservative sufficient condition and there mostly exist smaller µ that satisfies Theorem 1. In Section 4, we assess the sensitivity of our approach with respect to different choices of µ. Theorem 1 is concerned with the case where the initial point ˇP is feasible for the original problem (1a) - (1c). However, finding a feasible starting point can be difficult due to the presence of the non-convex quadratic inequality constraints (1b). The next theorem states that even if ˇP is not feasible, the proposed penalized relaxation can still result in a feasible point for the non-convex problem (1a) - (1c). Theorem 2. Consider an initial ˇP ∈ Sn,m that satisfies where κ is defined in (8c). If µ is sufficiently large, then the penalized convex relaxation (5a) - (5e) has a unique optimal solution ( ∗P , ∗X) that satisfies (3c) and (3d). Moreover, ∗P is feasible for (1a) - (1c). Proof. See Section 5 for the proof. ## 3.1 Sequential Penalized Relaxation Motivated by Theorems 1 and 2, this section presents a sequential approach that solves a sequence of penalized relaxations of the form (5a) - (5e) to infer high-quality feasible points for the non-convex problem (1a) - (1c). The proposed scheme starts from an initial point ˇP on the Stiefel manifold. In each round, the solution of the penalized relaxation (5a) - (5e) is projected onto the Stiefel manifold and then the projected point is employed as an initialization for the next round. Once a feasible point for (1a) - (1c) is obtained, according to Theorem 1, the proposed scheme preserves feasibility and generates a sequence of points whose objective values monotonically improves. The details of the sequential scheme are delineated in Algorithm 1. Figure 1: Two dimensional data representation on a training set from the synthetic data set. Left: MMDA [Bian and Tao, 2011], middle: SPRS, right: SPR-Q. The results show that the SPR-S and SPR-Q algorithms have provided more discriminative 2D representations compared to the MMDA method. ![Image](https://ijcai.org/proceedings/2019/fulltext/183/image_000000_60b397ab483ad37f49e7493e166b78ec67c6627492e149f16ad6bd7cc34f2d4a.png) ## Algorithm 1 Sequential Penalized Relaxation Input: ˇP ∈Sn,m, a fixed parameter µ > 0, and k = 0, - 1: repeat - 2: k ← k +1 - 3: P k ← solve (5a) - (5e) with the penalty µ〈P, ˇP〉 - 4: ˇP ← projSn,m P k - 5: until stopping criteria is met Output: P k ## 4 Experimental Results In this section, we conduct numerical experiments on real and synthetic datasets to verify the effectiveness of the proposed sequential approach, termed SPR, in solving non-convex optimization problems with orthogonality constraints. In Sections 4.1 and 4.2, we apply SPR on two practical problems involving orthogonality constraints. We use SPR-S and SPRQ to refer to the combination of Algorithm 1 with the SDP relaxation and the proposed convex quadratic relaxation, respectively. To solve the penalized relaxations in each round of the algorithm, we use MOSEK version 7.0 [Mosek, 2015]. Through the experiments, we leverage the inherent sparsity patterns of the problems to reduce the computational cost of solving large-scale semidefinite programs. This enables us to break down large-scale conic constraints into a set of smaller ones [Nakata et al., 2003]. Since finding a feasible point for (1a) - (1c) can be computationally demanding, we initialize Algorithm 1 with an arbitrary starting point on the Stiefel manifold and aim to improve the quality of the point. If the algorithm can recover a feasible point for (1a) - (1c), according to Theorem 1, it can generate a sequence of feasible points whose objective values monotonically improve. To measure the level of infeasibility, define tr{ ¯X - ¯P ¯P⊤} as the feasibility violation of an arbitrary feasible point ( ¯P, ¯X) of the problem (5a) - (5e). We terminate the sequential algorithm once the feasibility violation and objective value improvement are less than 10-5 or if the round number exceeds 100. Notice that the Nesterov acceleration technique can be employed to improve the convergence behaviour of the SPR algorithm. However, in this case, the algorithm may fail to preserve the monotonically decreasing order of the objective values even if the initial point is feasible. We apply the sequential algorithm on two fundamental machine learning problems of discriminative dimensionality reduction and graph clustering. Notice that each of these problems are well-studied in the literature and several approaches have been developed to efficiently target these applications. Therefore, it is not the intent of this work to compete with these state-of-the-art problem-specific approaches, but rather to demonstrate the potential of Algorithm 1 in solving the problems of form (1a) - (1c) that widely arise in different areas of machine learning. ## 4.1 Experiment I: Discriminative Dimensionality Reduction Given a collection of high-dimensional data points from c different classes, the problem of discriminative dimensionality reduction aims to learn a low-dimensional subspace on which the projection of different classes are well-separated. To this end, [Bian and Tao, 2011] proposed a max-min distance analysis (MMDA) that maximizes the minimum distance between all class pairs. This problem can be cast as a non-convex and non-smooth optimization problem of form where each Aij ∈ Sn is a given weighted distance matrix between the ith and j th classes. In this experiment, we evaluate the performance of the SPR algorithm for solving the problems of form (12a) - (12b). Closely related to our work, [Bian and Tao, 2011] uses a sequence of local SDP relaxations to find the solution of problem (12a) - (12b). We benchmark the SPR method against the MMDA on both real and synthetic datasets. To ensure the comparison is fair, both methods use the same initial point and the same distance matrices Aij which are computed based on [Bian and Tao, 2011]. Other parameter settings of the MMDA are set to their default values. Following [Bian and Tao, 2011], we conduct 100 independent experiments on 10-dimensional synthetic data from seven classes. For each class i, a mean vector ηi ∈ R10 is sampled from 10-dimensional zero mean Gaussian distribution with co-variance matrix 2I10 and then a pair of training and testing sets, each with 100 members, is generated based on the Gaussian distribution N(ηi, I10). Figure 2: Performance of SPR comparing to MMDA [Bian and Tao, 2011] on synthetic dataset (left), YALE dataset (right) [Georghiades et al., 2001]. Best viewed in color. ![Image](https://ijcai.org/proceedings/2019/fulltext/183/image_000001_6fd07cf49eca0560e0a73525fc14854755fbc468a6d812a83f631748149632ff.png) To compare the classification error rate, we project each test set into subspaces with varying dimensions, learned on its corresponding training set. The projected instances are then classified using the nearest mean classifier. Figure 2 (left) shows the average classification error rate with respect to the reduced dimensionality on the synthetic datasets. To run the experiment on the synthetic datasets, we set µ to 100 and 200 for SPR-S and SPR-Q, respectively. Moreover, we conduct this experiment on the YALE dataset consisting of 165 frontal face images of 15 individuals under different illumination and lightening conditions [Belhumeur et al., 1997]. Each image is of size 32×32 pixels. The results of this experiment are illustrated in Figure 2 (right). According to Figure 2, SPR-S and SPR-Q perform on par or better than the MMDA algorithm on both real and synthetic datasets in the problem of discriminative dimensionality reduction. In the experiment on the YALE dataset, we set µ to 5000 and 10000 for SPR-S and SPR-Q, respectively. To qualitatively compare the methods, Figure 1 visualizes the results of projecting a randomly chosen training set from the synthetic dataset on the 2D space. Observe that comparing to the MMDA method, the SPR-based algorithms learn more discriminative 2D representations that are suitable for classification tasks. To assess the sensitivity of the SPR algorithm with respect to the parameter µ, we perform the discriminative dimensionality reduction experiment with m = 2 on YALE dataset and report the results in Figure 3 for various choices of µ. Observe that the final solution obtained by the proposed algorithm is not very sensitive to the choice of µ. According to the figure, the SPR-S requires smaller values of µ to recover feasible points, e.g. µ = 5000, while SPR-Q fails to find feasible points for such choice of µ. Moreover, it can be seen that if µ exceeds a certain threshold, both SPR-S and SPR-Q provide the same sequence of feasible points. ## 4.2 Experiment II: Graph Clustering Given a weighted graph G with n vertices, the graph clustering problem aims to partition G into a set of sub-graphs such that the vertices within each one are more densely connected to each other than those belonging to different sub-graphs. Inspired by the well-known spectral clustering technique [Ng et al., 2002], this experiment incorporates a set of non-negative constraints to formulate the graph clustering problem as the following optimization [Han et al., 2017]: Figure 3: Sensitivity analysis of SPR-S (left) and SPR-Q (right) with respect to different choices of parameter µ for the discriminative dimensionality reduction problem, where m=2. This experiment is performed on the YALE dataset. Best viewed in color. ![Image](https://ijcai.org/proceedings/2019/fulltext/183/image_000002_e2797a685d4caca75a07db825feb8c197d323b38d273de346642abd3f8023f01.png) where L denotes the Laplacian matrix of the weighted graph G and ≥ is the element-wise inequality operator. Comparing to the spectral clustering, formulation (13a) - (13c) offers a more interpretable clustering framework which requires no further post-processing steps to identify the cluster members. Given ∗P, the optimal solution of the above problem, each vertex i is assigned to a cluster with label argmaxj ∗ Pij. [Han et al., 2017] proposed a fast and scalable heuristic, denoted by ONGR, to solve large-scale instances of the form (13a) - (13c). Due to the fact that this problem is a special case of (1a) - (1c), we apply the SPR algorithm to find the solution of (13a) - (13c) and use the same procedure as [Han et al., 2017] to create the Laplacian matrix L. To make a fair comparison between the ONGR and SPR, we use the same initialization for both methods. Table 1 reports the clustering performance of the SPR against [Han et al., 2017] on wellknown datasets from the UCI machine learning repository [Dua and Graff, 2017] and shape sets [Gionis et al., 2007; Jain and Law, 2005]. For each dataset, n, Dim, and m refer to the number of sample points, dimension of each point, and the number of classes, respectively. The scores for each method is computed by averaging over 30 independent runs for each dataset. As the results indicate, SPR-S and SPR-Q exhibit better performance compared to [Han et al., 2017] on most of the datasets. Through this experiment, we set µ= 1000 in the SPR algorithm and use the default parameter settings for the ONGR algorithm. ## 5 Proofs This section presents the proof of Theorems 1 and 2. Before proceeding with the proofs, we provide some prerequisite lemmas. We defer the proofs of the lemmas to the full version of the paper which would be posted on the authors webpage. Using the well-known epigraph technique [Boyd and Vandenberghe, 2004], the non-smooth term g0(P) in (3a) can be removed by adding a pair of linear constraints and incorporating an additional term into the objective function. This reformulation of (3a) - (3d) leads to the following penalized non-convex problem: Table 1: Clustering performance (%) on the UCI datasets [Dua and Graff, 2017] and shape sets [Gionis et al., 2007; Chang and Yeung, 2008; Jain and Law, 2005]. | Dataset | n | Dim. | m | ONGR | SPR-S | SPR-Q | |-------------|-----|--------|-----|--------|---------|---------| | Iris | 150 | 4 | 3 | 79.84 | 86.71 | 81.23 | | Spiral | 312 | 2 | 3 | 87.44 | 95.76 | 94.15 | | Jain | 373 | 2 | 2 | 88.42 | 92.33 | 90.26 | | Compound | 399 | 2 | 6 | 74.57 | 74.25 | 76.48 | | R15 | 600 | 2 | 15 | 86.07 | 85.36 | 86.94 | | Aggregation | 788 | 2 | 7 | 87.84 | 86.39 | 84.66 | with ¯γ ∈ Rw , ¯γ ∈ Rw , λ ∈ R|K| , and Ω ∈ Sm as the dual variables associated with the constraints (14b), (14c), (14d), and (14e), respectively. Observe that the problems (14a) - (14e) and (1a) - (1c) are equivalent, if µ=0. In what follows, we show that under certain conditions, the optimal solution of (14a) - (14e) can be obtained in polynomial time via convex relaxation. The next lemma guarantees the existence of Lagrange multipliers corresponding optimal solutions of the problem (14a) - (14e). Lemma 1. Consider an arbitrary ˇP ∈ Sn,m that satisfies If the following inequality holds true, then for every primal optimal pair ( ∗P, ∗t) of (14a) - (14e), there exists Lagrange multipliers ( ∗¯γ, ∗ ¯γ, ∗λ, ∗Ω) ∈ Rw ×Rw× R|K| ×Sm that satisfy the following Karush-Kuhn-Tucker (KKT) conditions where L(P, t, ¯γ, ¯γ, λ,Ω) represents the Lagrangian function of (14a) - (14e), defined as and, β and κ are defined in (8a) and (8c). The next lemma provides an upper bound on the Lagrange multipliers of the problem (14a) - (14e), that will be used to show that this problem can be relaxed to (5a) - (5e) with no effect on the solution. Lemma 2. Consider an arbitrary ˇP ∈ Sn,m that satisfies (15) and let µ satisfy (16). For every solution ( ∗P , ∗t) of (14a) - (14e), there exist Lagrange multipliers ( ∗¯γ, ∗ ¯γ, ∗λ, ∗Ω) ∈ Rw × Rw×Rp ×Sm that satisfy the KKT conditions (17a) - (17f) as well as the inequalities: where constant κ2 is given by and β, ψ and κ are defined in (8a) - (8c). Using Lemma 2, the next lemma offers conditions to guarantee that penalized relaxations give feasible points for (14a) - (14e). Lemma 3. Consider an initial point ˇP ∈ Sn,m and µ > 0. Let ( ∗P, ∗t) be a primal optimal solution of (14a) - (14e) with the corresponding Lagrange multipliers ( ∗¯γ, ∗ ¯γ, ∗λ, ∗Ω) that satisfy the KKT conditions (17a) - (17e). Define If the following inequalities hold true then the pair ( ∗P , ∗P ∗P⊤) is the unique primal solution of the penalized convex relaxation (5a) - (5d), where κ and κ2 are defined in (8c) and (19), respectively. Proof of Theorem 1. Due to the main assumption, it is straightforward to verify the following three inequalities: Consider an arbitrary optimal solution ( ∗P , ∗t) of (14a) - (14e). The point ∗P is consequently feasible for (1a) - (1c). Therefore dF( ˇP) = 0 and the inequalities (15) and (16) are satisfied. According to Lemma 2, there exist Lagrange multipliers ( ∗¯γ, ∗ ¯γ, ∗λ, ∗Ω) ∈ Rw ×Rw×Rp ×Sm corresponding to ( ∗P , ∗t) that satisfy the KKT conditions (17a) - (17f) as well as the inequalities (18a) and (18b). Based on Lemma 3 and since dF( ˇP) = 0, in order to prove the theorem, it suffices to show that: - (23a) is the direct consequence of (22a): - (23b) is the direct consequence of (18a), (22b), and (22c): - (23c) can be concluded from (18b), (22a), and (22c): Hence, according to Lemma 3, the point ( ∗P , ∗P ∗P⊤) is the unique optimal solution for the penalized relaxation (5a) - (5e), for which the relaxed constraints (3c) and (3d) are satisfied. Finally, due to the feasibility of pair ( ˇP, ˇP ˇP⊤), we have: and the proof is completed. Proof of Theorem 2. Consider an arbitrary optimal solution ( ∗P, ∗t) of (14a) - (14e). Due to the main assumption, (15) is satisfied and if µ is large, then (16) is satisfied as well. Moreover, according to Lemma 2, there exist Lagrange multipliers ( ∗¯γ, ∗ ¯γ, ∗λ, ∗Ω) ∈ Rw ×Rw×Rp ×Sm corresponding to ( ∗P , ∗t) that satisfy the KKT conditions (17a) - (17f) as well as the inequalities (18a) and (18b). According to Lemma 3, the proof follows directly from the fact that and therefore, if µ is sufficiently large, the inequalities (18a) and (18b) conclude (21a) - (21c). As a result, if µ is large, ( ∗P , ∗P ∗P⊤) is the unique primal solution of the penalized convex relaxation (5a) - (5d). ## 6 Conclusions This work introduces convex relaxations for solving a broad class of non-convex and non-smooth optimization problems involving orthogonality constraints. The proposed approach relies on solving a sequence of penalized convex relaxations to find feasible and near-globally optimal points for a given non-convex orthogonality-constrained problem. Experimental results on two fundamental problems in machine learning demonstrate the potential and effectiveness of the proposed approach in solving practical problems. ## Acknowledgements This work is in part supproted by a University of Texas System STARs award. ## References 1. [Abrudan et al., 2008] Traian E Abrudan, Jan Eriksson, and Visa Koivunen. Steepest descent algorithms for optimiza- tion under unitary matrix constraint. IEEE T SIGNAL PROCES, 56(3):1134-1147, 2008. 2. [Absil et al., 2009] P-A Absil, Robert Mahony, and Rodolphe Sepulchre. Optimization algorithms on matrix manifolds. Princeton University Press, 2009. 3. [Attouch et al., 2013] Hedy Attouch, Jérôme Bolte, and Be- nar Fux Svaiter. Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods. MATH PROGRAM, 137(1-2):91-129, 2013. 4. [Belhumeur et al., 1997] Peter N Belhumeur, João P Hes- panha, and David J Kriegman. Eigenfaces vs. fisherfaces: Recognition using class specific linear projection. IEEE TPAMI, (7):711-720, 1997. 5. [Bian and Tao, 2011] Wei Bian and Dacheng Tao. Max- min distance analysis by using sequential SDP relax- ation for dimension reduction. IEEE T PATTERN ANAL, 33(5):1037-1050, 2011. 6. [Boyd and Vandenberghe, 2004] Stephen Boyd and Lieven Vandenberghe. Convex optimization. Cambridge univer- sity press, 2004. 7. [Boyd et al., 2011] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein. Distributed opti- mization and statistical learning via the alternating direc- tion method of multipliers. Foundations and Trends R in Machine Learning, 3(1):1-122, 2011. 8. [Chang and Yeung, 2008] Hong Chang and Dit-Yan Ye- ung. Robust path-based spectral clustering. PATTERN RECOGN, 41(1):191-203, 2008. 9. [Chen et al., 2016] Weiqiang Chen, Hui Ji, and Yanfei You. An augmented Lagrangian method for ℓ 1 -regularized op- timization problems with orthogonality constraints. SIAM J SCI COMPUT, 38(4):B570-B592, 2016. 10. [Dua and Graff, 2017] Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. 11. [Edelman et al., 1998] Alan Edelman, Tomás A Arias, and Steven T Smith. The geometry of algorithms with orthog- onality constraints. SIAM J MATH ANAL, 20(2):303-353, 1998. 12. [Fortin and Glowinski, 2000] Michel Fortin and Roland Glowinski. Augmented Lagrangian methods: applica- tions to the numerical solution of boundary-value prob- lems, volume 15. Elsevier, 2000. 13. [Gao et al., 2018] Bin Gao, Xin Liu, Xiaojun Chen, and Ya- xiang Yuan. A new first-order algorithmic framework for optimization problems with orthogonality constraints. SIAM J OPTIMIZ, 28(1):302-332, 2018. 14. [Georghiades et al., 2001] Athinodoros S. Georghiades, Pe- ter N. Belhumeur, and David J. Kriegman. From few to many: Illumination cone models for face recognition un- der variable lighting and pose. IEEE T PATTERN ANAL, 23(6):643-660, 2001. 15. [Gionis et al., 2007] Aristides Gionis, Heikki Mannila, and Panayiotis Tsaparas. Clustering aggregation. ACM T KNOWL DISCOV D, 1(1):4, 2007. 16. [Han et al., 2017] Junwei Han, Kai Xiong, and Feiping Nie. Orthogonal and nonnegative graph reconstruction for large scale clustering. In IJCAI, 2017. [Ibaraki and Tomizuka, 2001] Soichi Ibaraki and Masayoshi Tomizuka. Rank minimization approach for solving BMI problems with random search. In ACC, volume 3, pages 1870-1875. IEEE, 2001. 17. [Jain and Law, 2005] Anil K Jain and Martin HC Law. Data clustering: A user's dilemma. In ICPRAI, 2005. 18. [Jiang et al., 2017] Bo Jiang, Jin Tang, Chris HQ Ding, and Bin Luo. Nonnegative orthogonal graph matching. In AAAI, 2017. 19. [Kheirandishfard et al., 2018a] Mohsen Kheirandishfard, Fariba Zohrizadeh, Muhammad Adil, and Ramtin Madani. Convex relaxation of bilinear matrix inequalities part II: Applications to optimal control synthesis. In CDC, 2018. 20. [Kheirandishfard et al., 2018b] Mohsen Kheirandishfard, Fariba Zohrizadeh, and Ramtin Madani. Convex relax- ation of bilinear matrix inequalities part I: Theoretical results. In CDC, 2018. [Kovnatsky et al., 2016] Artiom Kovnatsky, Klaus Glashoff, and Michael M Bronstein. MADMM: a generic algorithm for non-smooth optimization on manifolds. In ECCV, 2016. 21. [Lai and Osher, 2014] Rongjie Lai and Stanley Osher. A splitting method for orthogonality constrained problems. SIAM J SCI COMPUT, 58(2):431-449, 2014. 22. [Mosek, 2015] ApS Mosek. The MOSEK optimization tool- box for MATLAB manual. Version 7.1, 2015. 23. [Nakata et al., 2003] Kazuhide Nakata, Katsuki Fujisawa, Mituhiro Fukuda, Masakazu Kojima, and Kazuo Murota. Exploiting sparsity in semidefinite programming via ma- trix completion II: Implementation and numerical results. Mathematical Programming, 95(2):303-327, 2003. 24. [Ng et al., 2002] Andrew Y Ng, Michael I Jordan, and Yair Weiss. On spectral clustering: Analysis and an algorithm. In NIPS, pages 849-856, 2002. 25. [Ozolin ¸š et al., 2013] Vidvuds Ozolin ¸š, Rongjie Lai, Rus- sel Caflisch, and Stanley Osher. Compressed modes for variational problems in mathematics and physics. PNAS, 110(46):18368-18373, 2013. 26. [Tang and Liu, 2012] Jiliang Tang and Huan Liu. Unsuper- vised feature selection for linked social media data. In SIGKDD, 2012. 27. [Wen and Yin, 2013] Zaiwen Wen and Wotao Yin. A fea- sible method for optimization with orthogonality con- straints. MATH PROGRAM, 142(1-2):397-434, 2013. 28. [Yang et al., 2011] Yi Yang, Heng Tao Shen, Zhigang Ma, Zi Huang, and Xiaofang Zhou. ℓ 2,1 -norm regularized dis- criminative feature selection for unsupervised learning. In IJCAI, 2011. 29. [Yin et al., 2008] Wotao Yin, Stanley Osher, Donald Gold- farb, and Jerome Darbon. Bregman iterative algorithms for ℓ 1 -minimization with applications to compressed sens- ing. SIAM J IMAGING SCI, 1(1):143-168, 2008. 30. [Zhu et al., 2017] Hong Zhu, Xiaowei Zhang, Delin Chu, and Li-Zhi Liao. Nonconvex and nonsmooth optimiza- tion with generalized orthogonality constraints: An ap- proximate augmented Lagrangian method. SIAM J SCI COMPUT, pages 1-42, 2017.