Non-smooth Optimization over Stiefel Manifolds with Applications to Dimensionality Reduction and Graph Clustering

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.

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).

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]:

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:

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.