A Double Tracking Method for Optimization with Decentralized Generalized Orthogonality Constraints

Decentralized optimization has experienced significant advancements in recent decades, particularly in the Euclidean space. Various algorithms have been proposed to tackle different types of problems, such as gradient-based algorithms [24, 41, 27, 31], primal-dual frameworks [29, 22, 8, 18], and second-order methods [5, 43, 13]. In general, these algorithms are only capable of handling scenarios where variables are restricted in a convex subset of the Euclidean space. Consequently, these algorithms cannot be directly applied to solve (1.1b), where the feasible region is typically non-convex. Interested readers can refer to some recent surveys [25, 9] for more comprehensive information.

In order to solve decentralized optimization problems on manifolds, many algorithms have adopted geometric tools derived from Riemannian optimization, including tangent spaces and retraction operators. For instance, the algorithms delineated in [11, 14] extend the Riemannian gradient descent method [1] to the decentralized setting, which can be combined with gradient tracking techniques [32] to achieve the exact convergence. Building on this foundation, Chen et al. [10] further introduces the decentralized Riemannian conjugate gradient method. Additionally, there are several algorithms devised to address nonsmooth optimization problems, such as subgradient algorithm [34] and proximal gradient algorithm [38]. It is crucial to underscore that the computation of tangent spaces and retraction operators requires complete information about the matrix $M$, which is unattainable in a decentralized environment. This inherent limitation hinders the application of the aforementioned algorithms to the optimization problems with decentralized generalized orthogonality constraints.

There is another class of methodologies [35, 37, 36, 33] that focuses on employing infeasible approaches, such as augmented Lagrangian methods, to tackle nonconvex manifold constraints. These algorithms do not require each iterate to strictly adhere to manifold constraints, thereby allowing the pursuit of global consensus directly in the Euclidean space. Compared to Riemannian optimization methods, this type of algorithm requires only a single round of communication per iteration to guarantee their global convergence. However, it is noteworthy that these algorithms are tailored for Stiefel manifolds and cannot handle scenarios where the constraints themselves exhibit a distributed structure. Although we can draw inspiration from these algorithms to tackle the problem (1.1), the resulting penalty model remains challenging to solve. On the one hand, the penalty function usually forfeits the distributed structure inherent in the constraint, which is no longer separable across the agents. On the other hand, without full knowledge of the matrix $M$, each agent can not independently compute its own local gradients. Consequently, it is impossible to straightforwardly extend existing algorithms to solve the problem (1.1).

In decentralized optimization, current researches mainly focus on scenarios where the objective function exhibits a distributed structure. However, in problem (1.1), the generalized orthogonal constraint also exhibits a similar structure, which triggers off an enormous difficulty in solving it.

To develop a decentralized algorithm for the problem (1.1), we employ the constraint dissolving operator introduced in [40] to construct an exact penalty model, which can not be solved by existing algorithms. To efficiently minimize our proposed penalty model, we devise an approximate direction for the gradient of the penalty function, which is composed of separable components, including the gradient of the objective function and the Jacobian of the constraint mapping. Then, we propose to track these two components simultaneously across the network to assemble them in the approximate direction. This double-tracking strategy is quite efficient to reach a consensus on the generalized Stiefel manifolds.

Based on the aforementioned techniques, we develop a novel constraint dissolving algorithm with double tracking (CDADT). To the best of our knowledge, this is the first algorithm capable of solving optimization problems with decentralized generalized orthogonality constraints. Under rather mild conditions, we establish the global convergence of CDADT and provide its iteration complexity. Preliminary numerical results demonstrate the great potential of CDADT.

The rest of this paper is organized as follows. Section 2 draws into some preliminaries related to the topic of this paper. In Section 3, we develop a constraint dissolving algorithm with double tracking to solve the problem (1.1). The convergence properties of the proposed algorithm are investigated in Section 4. Numerical results are presented in Section 5 to evaluate the performance of our algorithm. Finally, this paper concludes with concluding remarks and key insights in Section 6.

The following notations are adopted throughout this paper. The Euclidean inner product of two matrices $Y_{1},Y_{2}$ with the same size is defined as $\left\langle Y_{1},Y_{2}\right\rangle=\mathrm{tr}(Y_{1}^{\top}Y_{2})$, where $\mathrm{tr}(B)$ stands for the trace of a square matrix $B$. The $p\times p$ identity matrix is represented by $I_{p}\in\mathbb{R}^{p\times p}$. We define the symmetric part of a square matrix $B$ as $\mathrm{sym}(B):=(B+B^{\top})/2$. The Frobenius norm and 2-norm of a given matrix $C$ are denoted by $\left\|C\right\|_{\mathrm{F}}$ and $\left\|C\right\|_{2}$, respectively. The $(i,j)$-th entry of a matrix $C$ is represented by $C(i,j)$. The notations $\mathbf{1}_{d}\in\mathbb{R}^{d}$ and $\mathbf{0}_{d}\in\mathbb{R}^{d}$ stand for the $d$-dimensional vector of all ones and all zeros, respectively. The notations $\sigma_{\min}(X)$ and $\sigma_{\max}(X)$ represent the smallest and largest singular value of a matrix $X$, respectively. The Kronecker product is denoted by $\otimes$. For any $d\in\mathbb{N}_{+}$, we denote $[d]:=\{1,2,\dotsc,d\}$. We define the distance between a point $X$ and a set $\mathcal{C}$ by $\mathrm{dist}(X,\mathcal{C}):=\inf\{\left\|Y-X\right\|_{\mathrm{F}}\mid Y\in\mathcal{C}\}$. Given a differentiable function $g(X):\mathbb{R}^{n\times p}\to\mathbb{R}$, the Euclidean gradient of $g$ with respect to $X$ is represented by $\nabla g(X)$. Further notations will be introduced wherever they occur.

We assume that the $d$ agents are connected by a communication network. And they can only exchange information with their immediate neighbors. The network $\mathtt{G}=(\mathtt{V},\mathtt{E})$ captures the communication links diffusing information among the agents. Here, $\mathtt{V}=[d]$ is composed of all the agents and $\mathtt{E}=\{(i,j)\mid i\text{~{}and~{}}j\text{~{}are connected}\}$ represents the set of communication links. Throughout this paper, we make the following assumptions on the network.

The conditions in Assumption 2 follow from standard assumptions in the literature [39, 25], which is dictated by the underlying network topology. According to the Perron-Frobenius Theorem [26], we find that the eigenvalues of $W$ fall within the range $(-1,1]$, and hence,

The parameter $\lambda$ serves as a key indicator of the network connectivity and is instrumental in the analysis of decentralized methods. Generally speaking, the closer $\lambda$ approaches $1$, the poorer the network connectivity becomes.

In this subsection, we delve into the first-order stationarity condition of the problem (1.1). Towards this end, we introduce some geometric concepts of Riemannian manifolds. For each point $X\in\mathcal{S}_{M}^{n,p}$, the tangent space to $\mathcal{S}_{M}^{n,p}$ at $X$ is referred to as $\mathcal{T}_{X}:=\{D\in\mathbb{R}^{n\times p}\mid D^{\top}MX+X^{\top}MD=0\}$. In this paper, we consider the Riemannian metric $\left\langle\cdot,\cdot\right\rangle_{M}$ on $\mathcal{T}_{X}$ that is induced from the inner product, i.e., $\left\langle V_{1},V_{2}\right\rangle_{M}=\left\langle V_{1},MV_{2}\right\rangle=\mathrm{tr}(V_{1}^{\top}MV_{2})$. The corresponding Riemannian gradient of a smooth function $f$ is given by

which is nothing but the projection of $\nabla f(X)$ onto $\mathcal{T}_{X}$ under the metric $\left\langle\cdot,\cdot\right\rangle_{M}$. Finally, the first-order stationarity condition of the problem (1.1) can be stated as follows.

Since this paper focuses on infeasible algorithms for the problem (1.1), we introduce the following definition of $\epsilon$-stationary point.

A part of the convergence results developed in this paper falls in the scope of a general class of functions that satisfy the Kurdyka-Łojasiewicz (KŁ) property [23, 20]. Below, we introduce the basic elements to be used in the subsequent theoretical analysis.

For any $\tau>0$, we denote by $\Phi_{\tau}$ the class of all concave and continuous functions $\phi:[0,\tau)\to\mathbb{R}_{+}$ which satisfy the following conditions,

1. (i)

   $\phi(0)=0$;

2. (ii)

   $\phi$ is continuously differentiable on $(0,\tau)$ and continuous at $0$;

3. (iii)

   $\phi^{\prime}(t)>0$ for any $t\in(0,\tau)$.

Now we define the KŁ property.

KŁ functions are ubiquitous in many practical applications, which covers a wealth of nonconvex nonsmooth functions. For example, tame functions constitutes a wide class of KŁ functions, including semialgebraic and real subanalytic functions. We refer interested readers to [6, 2, 3, 4, 7] for more details.

The constraint dissolving operator proposed in [40] offers a powerful technique to handle manifold constraints, which can be leveraged to construct an unconstrained penalty model for Riemannian optimization problems. Specifically, a constraint dissolving operator of the generalized orthogonality constraint (1.1b) is given by

Then solving the problem (1.1) can be converted into the unconstrained minimization of the following penalty function,

where $\beta>0$ is a penalty parameter.

Xiao et al. [40] have proved that (1.1) and (2.2) share the same first-order stationary points, second-order stationary points, and local minimizers in a neighborhood of $\mathcal{S}_{M}^{n,p}$. However, it is intractable to solve the problem (2.2) under the decentralized setting. It is worth noting that, in the construction of the penalty function $h(X)$, the constraint (1.1b) is integrated into the original objective function and the quadratic penalty term, resulting in the loss of the separable structure. To the best of our knowledge, there are currently no algorithms equipped to tackle such problems. Therefore, in the next section, we will design a novel algorithm to solve the penalty model under the decentralized setting.

Under the conditions in Assumption 1, the penalty function $h(X)$ in (2.2) is continuously differentiable, the gradient of which takes the following form:

Therefore, each agent $i$ can not compute the local gradient $\nabla f_{i}(Z)\mid_{Z=\mathcal{A}(X)}$ individually since the evaluation of $\mathcal{A}(X)$ requires the accessibility of $\{M_{i}\}$ over all the agents. In light of the property that $\left\|\mathcal{A}(X)-X\right\|_{\mathrm{F}}=\mathcal{O}(\left\|X^{\top}MX-I_{p}\right\|_{\mathrm{F}})$ whenever $X$ is not far away from $\mathcal{S}_{M}^{n,p}$, we propose to approximate the local gradient $\nabla f_{i}(Z)\mid_{Z=\mathcal{A}(X)}$ by $\nabla f_{i}(Z)\mid_{Z=X}$. Hereafter, $\nabla f_{i}(Z)\mid_{Z=X}$ is denoted by $\nabla f_{i}(X)$ for simplicity. As a result, we can obtain the following approximation of $\nabla h(X)$:

where

and

The approximate direction $H(X)$ of $\nabla h(X)$ possesses the following desirable properties.

Lemma 3.1 reveals that the norm of $\mathrm{grad}\,f(\mathcal{P}(X))$ and the feasibility violation of $X$ are both controlled by the norm of $H(X)$ as long as $X\in\mathcal{R}$ and $\beta$ is sufficiently large. Therefore, as an approximation of $\nabla h(X)$, $H(X)$ can serve as a search direction to solve the problem (2.2).

In the construction of $H(X)$, each agent $i\in[d]$ is able to compute the local gradient $\nabla f_{i}(X)$ independently. However, the evaluation of $H(X)$ still fails to be distributed into $d$ agents since it is not separable. Nevertheless, we find that the components constituting $H(X)$ exhibit a separable structure as follows,

where

This observation inspires us to propose a double-tracking strategy. Specifically, we can first track $U(X)$ and $V(X)$ separately across the whole network by resorting to the dynamic average consensus [45] protocol. Then, these two components are collected together to assemble a global estimate of $H(X)$. It is worth mentioning that $V(X)$ is exactly the Jacobian of the constraint mapping.

In this subsection, we describe the proposed algorithm to solve the problem (2.2). Hereafter, the notation $X_{i}^{(k)}$ represents the $k$-th iterate of $X_{i}$. Our algorithm introduces three auxiliary local variables $U_{i}^{(k)}\in\mathbb{R}^{n\times p}$, $V_{i}^{(k)}\in\mathbb{R}^{n\times p}$, and $H_{i}^{(k)}\in\mathbb{R}^{n\times p}$ for each agent $i$ at the $k$-th iteration. Specifically, $U_{i}^{(k)}$ and $V_{i}^{(k)}$ track $\nabla f(X_{i}^{(k)})$ and $MX_{i}^{(k)}$ respectively, through the exchange of local information. In addition, $H_{i}^{(k)}$ aims at estimating the search direction based on the formulation (3.3). The key steps of our algorithm from the perspective of each agent are outlined below.