On the Graph Fourier Transform for Directed Graphs

I Introduction

In this paper, we propose a novel alternative approach to build the GFT basis for the general case of directed graphs. Rather than starting from the decomposition of one of the graph matrix descriptors, either adjacency or Laplacian, we start identifying an objective function to be minimized and then we build an orthogonal matrix that minimizes that objective function. More specifically, we choose as objective function the graph cut size, as its minimization leads to identifying clusters. We consider the general case of directed graphs, which subsumes the undirected graphs as a particular case. The cut function is a set function and its minimization is NP-hard, however exploiting the sub-modularity property of the cut size, it has been shown that there exists a lossless convex relaxation of the cut size, named its Lovász extension , , whose minimization preserves the optimality of the solution of the original non-convex problem. Interestingly, the Lovász extension of the cut size gives rise to an alternative definition of total variation of a graph signal that captures the edges’ directivity. Furthermore, in the case of undirected graphs, the Lovász extension reduces to the $l_{1}$ norm total variation of a graph signal, which represents the discrete counterpart of the total variation of continuous-time signals, which plays a fundamental role in the continuous time Fourier Transform, see, e.g., ,. We define the GFT basis as the set of orthonormal vectors that minimize the Lovász extension of the cut size. Unfortunately, even though the objective function is convex, the resulting problem is non-convex, because of the orthogonality constraint imposed on the basis vectors. Thus, to find a (possibly local) solution of the problem in an efficient manner, we exploit two recently developed methods that are specifically tailored to handle non-convex orthogonality constraints, namely, the splitting orthogonality constraints (SOC) method , and the proximal alternating minimized augmented Lagrangian (PAMAL) method . SOC method is quite simple to implement and, even if no convergence proof has been provided yet, extensive numerical results validate the effectiveness and robustness of such a strategy. Conversely, PAMAL algorithm, which hybridizes the augmented Lagrangian method and the proximal minimization scheme, is known to guarantee convergence. Furthermore, any limit point of each sequence generated by PAMAL method satisfies the Karush-Kuhn Tucker conditions of the original non-convex problem . Finally, to prevent the resulting basis vectors to be excessively sparse vectors, we consider the minimization of a continuous relaxation of the balanced cut size. To solve the corresponding non-convex fractional problem, we adopt an efficient and convergent algorithm based on the explicit-implicit gradient method .

The paper is organized as follows. Sec. II introduces the graph signal variations as the continuous Lovász extension of the min-cut size. In Sec. III, we define the GFT as the set of optimal orthonormal vectors minimizing the graph signal variation, and in Sec. IV we illustrate the optimization methods used for solving the resulting non-convex problem. Therefore, in Sec. V we conceive the GFT as the solution of a balanced min cut problem, while Sec. VI illustrates some numerical examples validating the effectiveness of the proposed approaches. Finally, Sec. VII draws some conclusions.

II Min-cut size and its Lovász extension

One of the basic operations over graphs is clustering, i.e. the partition of the graph onto disjoint subgraphs, such that the vertices within each subgraph (cluster) are highly interconnected, whereas there are only a few links between different clusters. Finding a good partition can be formulated as the minimization of the cut size , whose definition is reported here below. Let us consider a subset of vertices $\mathcal{S}\subset\mathcal{V}$, and its complement set in $\mathcal{V}$ denoted by $\bar{\mathcal{S}}$. The edge boundary of $\mathcal{S}$ is defined as the set of edges with one end in $\mathcal{S}$ and the other end in $\bar{\mathcal{S}}$. The cut size between $\mathcal{S}$ and $\bar{\mathcal{S}}$ is defined as the sum of the weights over the boundary , i.e.

Note that f(\mbox{\boldmathx}) is piecewise affine w.r.t. $x$, and $F(\mathcal{S})=f(\mathbf{1}_{\mathcal{S}})$ for all $\mathcal{S}\subseteq\mathcal{V}$. An interesting class of set functions is given by the submodular set functions, whose definition follows next.

A fundamental property of a submodular set function is that its Lovász extension is a convex function. This is formally stated in the following proposition [24, p.23].

Let $F:2^{\mathcal{V}}\rightarrow\mathcal{R}$ be a submodular function and $f$ be its Lovász extension. Then, it holds

Moreover, the set of minimizers of f(\mbox{\mathbf{x}}) on $^{N}$ is the convex hull of the minimizers of f(\mbox{\mathbf{x}}) on $\{0,1\}^{N}$.

The cut size function in (1) is known for being submodular, see, e.g., , . More specifically, as shown in [24, p.54], the cut function is equal to the positive linear combination of the function $G_{ij}\,:\,\mathcal{S}\mapsto(\mathbf{1}_{\mathcal{S}})_{i}[1-(\mathbf{1}_{\mathcal{S}})_{j}]$, i.e.

The function $G_{ij}$ is the extension to $\mathcal{V}$ of a function $\widetilde{G}_{ij}$ defined only on the power set of $\{i,j\}$, where $\widetilde{G}_{ij}(\{i\})=1$ and all other values are zero, so that, from (2), its Lovász extension is $\widetilde{G}_{ij}(x_{i},x_{j})=[x_{i}-x_{j}]_{+}$ with $[y]_{+}:=\text{max}\{y,0\}$. Therefore the Lovász extension of the cut size function, in the general case of directed graphs, is given by:

We term this function the Graph Directed Variation (GDV), as it captures the edges’ directivity. For undirected graphs, imposing $a_{ij}=a_{ji}$, the Lovász extension of the cut size boils down to

Interestingly, this function, which we call Graph Absolute Variation (GAV), represents the discrete counterpart of the $l_{1}$ norm total variation, which plays a key role in the classical Fourier Transform of continuous time signals , . It is easy to show that the directed variation GDV satisfies the following properties:

\text{GDV}(\mbox{\mathbf{x}})=0, \forall\,\mbox{\mathbf{x}}=c\mathbf{1} with $c\geq 0$;

\text{GDV}(\alpha\,\mbox{\mathbf{x}})=\alpha\,\text{GDV}(\mbox{\mathbf{x}}), $\forall\,\alpha\geq 0$, i.e. it is positively homogeneous;

GDV is neither a proper norm nor a semi-norm, since, in this latter case, it should be absolutely homogeneous. However, it meets the desired property ii) ensuring that a constant graph signal has zero total variation.

III Graph Fourier Basis and Directed Total Variation

Alternative definitions of GFT have been proposed in the literature, depending on the different perspectives used to emphasize specific signal features. In case of undirected graphs, the GFT of a vector $s$ was defined as

where the columns of matrix $\mathbf{U}$ are the eigenvectors of the Laplacian matrix $\mathbf{L}$, i.e. \mbox{\mathbf{L}}=\mbox{\mathbf{U}}\mbox{\mathbf{\Lambda}}\mbox{\mathbf{U}}^{T}. This definition is basically rooted on the clustering properties of these eigenvectors, see, e.g., . In fact, by definition of eigenvector, the Fourier basis used in (5) can be thought as the solution of the following sequence of optimization problems:

where $\mathbf{V}$ comes from the Jordan decomposition of the nonsymmetric adjacency matrix $\mathbf{A}$, i.e. \mbox{\mathbf{A}}=\mbox{\mathbf{V}}\mbox{\mathbf{J}}\mbox{\mathbf{V}}^{-1}. To estimate variations of the graph Fourier basis and to identify an order among frequencies, the total variation of a vector was defined in as

where \mbox{\mathbf{A}}_{\text{norm}}\!:=\!\mbox{\mathbf{A}}/|\lambda_{\text{max}}(\mbox{\mathbf{A}})|. The previous definition leads to the elegant theory of algebraic signal processing over graphs . However, there are some critical issues associated to that definition that need to be further explored. First, the definition of total variation as given in (8) does not ensure that a constant graph signal has zero total variation, and this collides with the common meaning of total variation , , . Second, the columns of $\mathbf{V}$ are linearly independent complex generalized eigenvectors, but in general they are not orthogonal. This gives rise to a GFT that does not preserve inner products when passing from the observation to the transformed domain. Furthermore, the computation of the Jordan decomposition incurs into serious and intractable numerical instabilities when the graph size exceeds even moderate values and more stable matrix decomposition methods have to be adopted to tackle its instability issues . To overcome some of these criticalities, very recently the authors of proposed a shift operator based on the directed Laplacian of a graph. Using the Jordan decomposition, the graph Laplacian is decomposed as

To quantify oscillations in the graph harmonics and to order the frequencies, the total variation was defined in as

This definition of total variation ensures a zero value for constant graph signals. Furthermore, the eigenvalues with small absolute value correspond to low frequencies. Nevertheless, the GFT given by \mbox{\mathbf{F}}=\mbox{\mathbf{V}}^{-1}_{L} is still a non-unitary transform and its computation is affected by the numerical instabilities associated to the Jordan decomposition.

The constraints are used to find an orthonormal basis and to prevent the trivial null solution. Although the objective function is convex, problem $\mathcal{P}$ is non-convex due to the orthogonality constraint. In the next section, we present two alternative optimization strategies aimed at solving the non-convex, non-differentiable problem $\mathcal{P}$ in an efficient manner.

IV Optimization Algorithms

To avoid handling the non-convex orthogonality constraints directly, several methods have been proposed in the literature based on the solution of a sequence of unconstrained problems approaching the feasibility condition, such as the penalty methods , and the augmented Lagrangian based methods , . The penalty method is generally simple, but it suffers from slow-convergence and ill-conditioning. On the other hand, the standard augmented Lagrangian method solves a sequence of sub-problems that usually have no analytical solutions and the choice of the initial points, ensuring a fast convergence rate, is usually nontrivial. To cope with these issues, in this section we present two alternative iterative algorithms to solve the non-convex, non-smooth problem $\mathcal{P}$, hinging on some recently developed methods for solving non-differentiable problems with non-convex constraints ,. The first method, introduced in , called splitting orthogonality constraints (SOC) method, is based on the alternating method of multipliers (ADMM) , and the split Bregman method ,. The SOC method leads to some important benefits, as it is simple to implement and the resulting non-convex sub-problem with orthonormal constraint admits a closed form solution. Although no convergence proof of SOC method has been provided yet, numerical results validate its value and robustness. An alternative optimization method that tackles the non-convex minimization problem $\mathcal{P}$ and guarantees convergence is the PAMAL algorithm recently developed in . The algorithm combines the augmented Lagrangian method with proximal alternating minimization. A convergence proof was provided in . More specifically, this method has the so-called sub-sequence convergence property, i.e. there exists at least one convergent sub-sequence, and any limit point satisfies the Karush-Kuhn Tucker (KKT) conditions of the original nonconvex problem. Building on these algorithms, in the sequel we introduce two efficient optimization strategies that build the basis for the Graph Fourier Transform, as the solution of problem $\mathcal{P}$.

The SOC algorithm was developed in and tackles orthogonality constrained problems by iteratively solving a convex problem and a quadratic problem that admits a closed-form solution. More specifically, introducing an auxiliary variable \mbox{\mathbf{P}}=\mbox{\mathbf{X}} to split the orthogonality constraint, problem $\mathcal{P}$ is equivalent to

The first constraint is linear and, as discussed in , it can be solved using Bregman iteration. Therefore, by adding the Bregman penalty function , problem (12) is equivalent to the following simple two-step procedure:

where $\beta$ is a strictly positive constant. Similarly to ADMM and split Bregman iteration , the above problem can be solved by iteratively minimizing with respect to $\mathbf{X}$ and $\mathbf{P}$:

The interesting aspect of this formulation is that subproblem $\mathcal{P}_{k}$ is convex and the second constrained quadratic problem $\mathcal{Q}_{k}$ has a closed-form solution, as illustrated in the following proposition.

Define \mbox{\mathbf{Y}}^{k}=\mbox{\mathbf{X}}^{k}+\mbox{\mathbf{B}}^{k-1} and let

Combining (13) and Proposition 2, the main steps of the SOC method are summarized in Algorithm 1. It is important to remark that the choice of the coefficient $\beta$ strongly affects the convergence behavior of the algorithm: a large value of $\beta$ will force a stronger equality constraint, while a too small $\beta$ might not be able to guarantee the solution to satisfy the orthogonality constraint. Then, a proper tuning of the coefficient $\beta$ is important to ensure a fast convergence of the algorithm. Although, as remarked in , the convergence analysis of SOC algorithm is still an open problem, we will show next that the numerical results testify the validity and robustness of this method when applied to our case.

IV-B PAMAL method

Given these symbols, problem (12) is equivalent to the following one:

The basic idea to solve a problem in the form of $\mathcal{P}_{e}$ was proposed in , and combines the augmented Lagrangian method , with the alternating proximal minimization algorithm. The result is known as the PAM method , which deals with non-smooth, non-convex optimization. According to the augmented Lagrangian method, we add a penalty term to the objective function in order to associate a high cost to unfeasible points. In particular, the augmented Lagrangian function associated to the non-smooth problem $\mathcal{P}_{e}$, is

Compute the critical point (\mbox{\mathbf{X}}^{k},\mbox{\mathbf{P}}^{k}) of the function \mathcal{L}(\mbox{\mathbf{X}},\mbox{\mathbf{P}},\mbox{\mathbf{\Lambda}}^{k};\rho^{k}) by solving

Update the multiplier estimates \mbox{\mathbf{\Lambda}}^{k};

We will show next how to implement the previous steps, which are described in detail in Algorithm 2. Computation of the critical points (\mbox{\mathbf{X}}^{k},\mbox{\mathbf{P}}^{k}). The optimal solution (\mbox{\mathbf{X}}^{k},\mbox{\mathbf{P}}^{k}) of problem (15) is computed using an approximate algorithm, i.e. finding a subgradient point \mbox{\mathbf{\Theta}}^{k}\in\partial\mathcal{L}(\mbox{\mathbf{X}}^{k},\mbox{\mathbf{P}}^{k},\mbox{\mathbf{\Lambda}}^{k};\rho^{k}) satisfying, with a prescribed tolerance value $\epsilon^{k}$, the following inequality

with \mbox{\mathbf{P}}^{k}\in\mathcal{S}_{t}. To evaluate such point, we exploit a coordinate-descent method with proximal regularization based on the PAM method proposed in . More specifically, at the $k$-th outer iteration of the algorithm, we compute (\mbox{\mathbf{X}}^{k},\mbox{\mathbf{P}}^{k}) by iteratively solving, at each inner iteration $n$, the following proximal regularization of a two blocks Gauss-Seidel method:

where the proximal parameters $c_{i}^{k,n}$ can be arbitrarily chosen as long as they satisfy

The convergence claim for Algorithm 2 to a stationary solution of problem $\mathcal{P}_{e}$ is stated in the following theorem.

The proof follows similar arguments as in [26, Th. 3.1-3.5], and thus is omitted due to space limitation. ∎

Remark $1$. Note that both Algorithms $1$ and $3$ at each step of their loops have to compute the SVD of an $N\times N$ matrix. Therefore, at each iteration their computational cost is proportional to $\mathcal{O}(N^{3})$. So, clearly, there is a complexity issue that deserves further investigations to enable the application to large size graphs. In this paper, we have not investigated methods to reduce the complexity of the approach exploiting, for instance, the sparsity of the graphs under analysis. Also, we have not optimized the selection of the parameters involved in both SOC and PAMAL methods. However, even if complexity is an issue, the proposed approach is more numerically stable than the only method available today for the analysis of directed graphs, based on the Jordan decomposition.

Remark $2$. The two alternative methods proposed above to solve the non-convex problem $\mathcal{P}$ are robust to random initializations, as testified also by the numerical results presented in the sequel. In terms of implementation complexity, SOC algorithm is easier to code even though, to the best of our knowledge, a theoretical proof of its convergence is still lacking.

V Minimization of balanced total variation

The minimization of the total variation as in (III) is inspired by the min-cut problem. However, in some cases, this might favor the appearance of very sparse vectors or of very small clusters, possibly also isolated nodes. One way to prevent these undesired solutions passes through the introduction of the balanced cut , . A popular definition for the balanced cut of undirected graph is the Cheeger cut , which is given by:

Note that $\min(|\mathcal{S}|,|\bar{\mathcal{S}}|)$ attains its maximum when $|\mathcal{S}|=|\bar{\mathcal{S}}|=N/2$, so that, for a given value of $\text{cut}(\mathcal{S},\bar{\mathcal{S}})$, the minimum occurs when $\mathcal{S}$ and $\bar{\mathcal{S}}$ have approximately equal size. While the problem stated above is NP-hard, a tight continuous relaxation of the balanced cut problems has recently been shown to provide excellent clustering results . In , it was proved that the balanced Cheeger cut problem in (21) for undirected graphs admits the following exact continuous relaxation

where \text{m}(\mbox{\mathbf{x}}) stands for the median value of $\mathbf{x}$. Note that since it holds \sum_{i}\mid x_{i}-\text{m}(\mbox{\mathbf{x}})\mid=0, \forall\mbox{\mathbf{x}}\in\text{span}\{\mathbf{1}\}, problem (22) is well-defined if \mbox{\mathbf{x}}\perp\mathbf{1}. Then, the problem in (22) can be recast as:

In it was proved that (22) is an exact relaxation of the Cheeger cut problem and, for any minimizer $\mathbf{x}$, there is a number $\nu$ such that, $\forall i$, the binary solution ${x}_{\nu}(i)=1$ if ${x}(i)>\nu$ and ${x}_{\nu}(i)=0$ for ${x}(i)\leq\nu$, is also a minimizer of the Cheeger cut problem. Then, from the equivalence of problems (22) and (23), this result holds true also for any minimizer of (23). In the sequel, we formulate the problem of finding the Fourier basis minimizing the balanced total variation in both cases of directed and undirected graphs. To this end, let us define the function

where f(\mbox{\mathbf{x}}_{k})=\text{GAV}(\mbox{\mathbf{x}}_{k}) in (4), or f(\mbox{\mathbf{x}}_{k})=\text{GDV}(\mbox{\mathbf{x}}_{k}) in (3), in case of undirected or directed graphs, respectively. According to problem (22), we can find a set of Fourier bases \{\mbox{\mathbf{x}}_{k}\}_{k=1}^{N}, with \mbox{\mathbf{x}}_{1}=b\mathbf{1}, by iteratively solving, for $k=2,\ldots,N$, the following problem

Note that problem $\mathcal{P}^{b}_{k}$ is non-convex in both the constraints set and the objective function. Recently, several algorithms , , , , have been proposed to minimize relaxations of the balanced cut problem that are similar to (22). Typically, these algorithms give excellent numerical performance, although theoretical convergence proofs are not available. For instance, in , the authors proposed an algorithm minimizing (22), along with a proof of convergence to a critical point of the original problem. This method is a new steepest descent algorithm based on the explicit-implicit gradient of the function \text{E}(\mbox{\mathbf{x}}_{k})\triangleq\frac{f(\mbox{\mathbf{x}}_{k})}{\text{B}(\mbox{\mathbf{x}}_{k})} where \text{B}(\mbox{\mathbf{x}}_{k})=\sum_{i}\mid x_{k}(i)-\text{m}(\mbox{\mathbf{x}}_{k})\mid. The explicit-implicit subgradient of the non-smooth function \text{E}(\mbox{\mathbf{x}}_{k}) is given by

Let us now consider the following proximal minimization problem

Any stationary solution of (28) will be also solution of the subgradient equation

Replacing in this last equality the expression of \mbox{\mathbf{x}}_{k}^{n+1} given in (27), we obtain the following set of two equations to be iteratively updated:

The formal description of the iterative optimization method is given in Algorithm 4, where we denote by \text{sign}(\mbox{\mathbf{a}}) and \text{mean}(\mbox{\mathbf{a}}), respectively, the element-wise sign and the mean value of a vector $\mathbf{a}$. The convergence analysis of the algorithm to a critical point of E was derived in , for undirected graphs. However, since for directed graphs f(\mbox{\mathbf{x}}_{k}) preserves all the required properties (i.e., it is non-smooth and convex), the convergence results in , hold also for the minimization of the balanced directed variation.

VI Numerical results

In Fig. 4, we report the optimal basis, computed using Algorithm 2, for the graph with a directed cycle depicted in Fig. 1(c). Interestingly, in this case, there can only be one vector that yields zero directed variation: the constant vector. In fact, the cyclical structure of the graph now prevents the existence of non-constant vectors able to null the directed variation. The properties described above are a unique and an interesting consequence of the edge directivity. In fact, as can be observed from Fig. 5, the optimal bases for the corresponding undirected graph (obtained by simply removing edge directivity) have only one vector with zero variation, the constant vector. Conversely, in the case shown before, we have had three, two, and one vectors yielding zero variation.

Convergence test. Since the optimization problem $\mathcal{P}$ is non-convex, there is of course the possibility that the proposed methods fall into a local minimum. Furthermore, while PAMAL method guarantees convergence, SOC algorithm might also fail to converge because, theoretically speaking, there is no convergence analysis. To test what happens, we considered several independent initializations of both SOC and PAMAL algorithms in the search for a basis for the graph of Fig. 1(a). In Fig. 6, we report the average behavior ($\pm$ the standard deviation) of the directed variation versus the iteration index $m$, which counts the overall number of (outer and inner) iterations for Algorithm 1 and 2. The curves refer to $200$ independent initializations of algorithms SOC and PAMAL, using the same initialization for both. We can observe that in all cases the algorithms converge but indeed there is a spread in the final variation, meaning that both methods can incur into local minima. Nonetheless, the spread is quite limited, which suggests that bases associated to different local minima behave similarly in terms of total variation. Additionally, since the PAMAL algorithm solves the orthogonality constrained, non-convex problem by iteratively updating the primal variables and the multipliers, the objective function evaluated at each (inner and outer) iteration does not necessarily follow a monotonous decay, as can be noticed by the lower subplot in Fig. 6.

Comparison with alternative GFT bases. We compare now the GFT basis found with our methods with the bases associated to either the Laplacian or the adjacency matrix, as proposed in , and references therein. To compare the results, we applied all algorithms to several independent realizations of random graphs. We chose as family of random graphs the so called scale-free graphs, as they are known to fit many situations of practical interest . In the generation of random scale-free graphs, it is possible to set the minimum degree $d_{min}$ of each node. To compare our method with the GFT definition proposed in , since the eigenvectors of an asymmetric matrix can be complex and the directed total variation GDV, as defined in (3), does not represent a valid metric for complex vectors, we restricted the comparison to undirected scale-free graphs, in which case the adjacency and Laplacian matrices are real and symmetric, so that their eigenvectors are real. In the sequel, we will use the notations \text{GAV}(\mbox{\mathbf{X}}):=\sum_{k=1}^{N}\text{GAV}(\mbox{\mathbf{x}}_{k}) and \text{GQV}(\mbox{\mathbf{X}}):=\sum_{k=1}^{N}\text{GQV}(\mbox{\mathbf{x}}_{k}) to denote, respectively, the total graph absolute and quadratic variation of a matrix $\mathbf{X}$. In Fig. 7, we compare the following metrics: a) \text{GAV}(\mbox{\mathbf{X}}^{*}), derived by solving problem $\mathcal{P}$ through the SOC and PAMAL methods; b) \text{GAV}(\mbox{\mathbf{V}}), where $\mathbf{V}$ are the eigenvectors of the adjacency matrix according to the GFT defined in (7); c) \text{GAV}(\mbox{\mathbf{U}}), where $\mathbf{U}$ are the eigenvectors of the Laplacian matrix by assuming the GFT as in (5), that for undirected graphs is equivalent to the GFT defined in (10). More specifically, Fig. 7 shows the previous metrics vs. the minimum degree of the graph averaged over $100$ independent realizations of scale-free graphs of $N=20$ nodes. As we can notice from Fig. 7, the bases built using SOC and PAMAL algorithms yield a significantly lower total variation than the conventional bases built with either adjacency or Laplacian eigenvectors. This is primarily due to the fact that our optimization methods tend to assign constant values within each cluster. Finally, in Fig. 8 we compare the alternative basis vectors using as performance metric the GQV. So, in Fig. 8 we report the \text{GQV}(\mbox{\mathbf{X}}^{*}) metric derived from the SOC and PAMAL methods with \text{GQV}(\mbox{\mathbf{V}}) and \text{GQV}(\mbox{\mathbf{U}}) obtained, respectively, from the eigenvectors of the adjacency and the Laplacian matrix. Again, the results are averaged over $100$ independent realizations of scale-free graphs, vs. the average minimum degree, under the same settings of Fig. 7. Interestingly, even if our basis vectors \mbox{\mathbf{X}}^{*} do not coincide with $\mathbf{V}$ or $\mathbf{U}$, they provide the same GQV, within negligible numerical inaccuracies. Indeed, the invariance of the metric \text{GQV}(\mbox{\mathbf{X}}), for any square, orthogonal matrix $\mathbf{X}$, can be easily proved from the equality \text{GQV}(\mbox{\mathbf{X}})=\sum_{k=1}^{N}\mbox{\mathbf{x}}^{T}_{k}\mbox{\mathbf{L}}\mbox{\mathbf{x}}_{k}=\text{trace}(\mbox{\mathbf{X}}^{T}\mbox{\mathbf{L}}\mbox{\mathbf{X}}), by observing that \text{trace}(\mbox{\mathbf{X}}^{T}\mbox{\mathbf{L}}\mbox{\mathbf{X}})=\text{trace}(\mbox{\mathbf{L}}) for any orthogonal matrix $\mathbf{X}$. Interestingly, this implies that, for undirected graphs, our orthogonal matrix \mbox{\mathbf{X}}^{*} can be obtained by applying an orthogonal transform to the Laplacian eigenvectors basis.

Complexity issues. Clearly, looking at both SOC and PAMAL methods, complexity is a non trivial issue which deserves further investigations, especially when the size of the graph increases. To get an idea of computing time, in Fig. 9 we report the execution time of both SOC and PAMAL algorithms, as a function of the number of vertices in the graph. The results have been obtained running a non-compiled Matlab program, with no optimization of the parameters involved, by setting $\rho^{1}=\beta=20$. The program ran on a laptop having a processor Intel Core i7-4500, CPU 1.8, 2.4 GHz. The graphs under test were generated as geometric random graphs with equal percentage of directed links as $N$ increases.

Examples with real networks. As an application to real graphs, in Fig. 10 we considered the directed graph obtained from the street map of Rome, incorporating the true directions of traffic lanes in the area around Mazzini square. The graph is composed of $239$ nodes. Even though, the scope of this paper is to propose a method to build a GFT basis, so that we do not dig further into applications, this an example that has interesting applications of GSP. The problem in this case is to build a map of vehicular traffic in a city, starting from a subset of measurements collected along road side units or sent by cars equipped with ad hoc equipment. The problem can be interpreted as the reconstruction of the entire graph signal from a subset of samples and then it builds on graph sampling theory . In Fig. 11 we report some basis vectors obtained by using Algorithm $2$ with $\rho^{1}=10$. We can observe that the basis vectors highlight clusters, while capturing the edges’ directivity.

Balanced total variation. In some cases, the solution of the total variation problem in (III) can cut the graph in subsets of very different cardinality. As an extreme case, it may be not uncommon to have a subset composed of only one node and the other set containing all the rest of the network. To prevent such a behavior, Algorithm 4 aims at minimizing the balanced total variation. An example of its application to the graph of Fig. 10 is reported in Fig. 12, where we show some basis vectors computed using Algorithm 4. Comparing these vectors with the corresponding ones obtained with PAMAL algorithm, see, e.g. Fig. 11, we can see how clusters of single nodes are now avoided.

VII Conclusion

where g_{k,n-1}(\mbox{\mathbf{P}})\triangleq\langle\mbox{\mathbf{\Lambda}}^{k},\mbox{\mathbf{P}}-\mbox{\mathbf{X}}^{k,n-1}\rangle+\frac{\rho^{k}}{2}\|\mbox{\mathbf{P}}-\mbox{\mathbf{X}}^{k,n-1}\|^{2}_{F}+\frac{c_{2}^{k,n-1}}{2}\parallel\mbox{\mathbf{P}}-\mbox{\mathbf{P}}^{k,n-1}\parallel^{2}_{F}. Our proof consists of two steps: i) first, we find the stationary solutions by solving the KKT necessary conditions; ii) then, we prove that the resulting closed-form solution is a global minimum of the non-convex problem (31). The Lagrangian function $\mathcal{L}_{P}$ associated to (31) can be written as

where we chose \mbox{\mathbf{\Lambda}}_{1}=\mbox{\mathbf{\Lambda}}_{1}^{T}. Hence, defining \mbox{\mathbf{B}}\triangleq\mbox{\mathbf{I}}+2\mbox{\mathbf{\Lambda}}_{1}/(\rho^{k}+c_{2}^{k,n-1}), from equation a) one gets:

with \mbox{\mathbf{F}}\triangleq\displaystyle\frac{c_{2}^{k,n-1}\mbox{\mathbf{P}}^{k,n-1}+\rho^{k}\mbox{\mathbf{X}}^{k,n-1}-\mbox{\mathbf{\Lambda}}^{k}}{\rho^{k}+c_{2}^{k,n-1}}. Let \mbox{\mathbf{Q}}\mbox{\mathbf{\Sigma}}\mbox{\mathbf{T}}^{T} be the SVD decomposition of $\mathbf{F}$. From (34), it turns out

and, using the orthogonality condition b) in (33), it holds

Therefore, replacing $\mathbf{B}$ in (35), we get

It remains to prove that \mbox{\mathbf{P}}^{\star}=\mbox{\mathbf{P}}^{k,n}=\mbox{\mathbf{Q}}\mbox{\mathbf{T}}^{T} is a global minimum for problem (31). To this end, it is sufficient to show that

i.e., using the equalities \parallel\mbox{\mathbf{P}}^{\star}\parallel_{F}^{2}=\parallel\mbox{\mathbf{P}}\parallel_{F}^{2}=N, we have to prove that \forall\,\mbox{\mathbf{P}}\,:\,\mbox{\mathbf{P}}^{T}\mbox{\mathbf{P}}=\mbox{\mathbf{I}}, it results

Using the above definition of $\mathbf{F}$, (39) reduces to

and since \mbox{\mathbf{P}}^{\star}=\mbox{\mathbf{Q}}\mbox{\mathbf{T}}^{T}, the final inequality to hold true is

Define \mbox{\mathbf{Z}}^{T}:=\mbox{\mathbf{T}}^{T}\mbox{\mathbf{P}}^{T}\mbox{\mathbf{Q}} so that \mbox{\mathbf{Z}}^{T}\mbox{\mathbf{Z}}=\mbox{\mathbf{I}}. Then, from (41) we get

This last inequality holds because $\Sigma_{ii}>0$ and $Z_{ii}\leq\mid Z_{ii}\mid\;\;\leq 1$, $\forall i$, where the latter is implied by \mbox{\mathbf{Z}}^{T}\mbox{\mathbf{Z}}=\mbox{\mathbf{I}} . Additionally, $Z_{ii}=1$, $\forall i$, if and only if \mbox{\mathbf{Z}}=\mbox{\mathbf{I}}, so that the equality in (42) holds if and only if \mbox{\mathbf{Z}}=\mbox{\mathbf{I}} or \mbox{\mathbf{P}}^{\star}=\mbox{\mathbf{Q}}\mbox{\mathbf{T}}^{T}.

-B Proof of Theorem 1

For lack of space, we omit here the details of the proof, which proceeds using similar arguments as in the proof of Proposition $2.5$ in . However, to invoke this correspondence, we need to prove that the following properties hold true: i) the function $\mathcal{L}_{k}$ in (IV-B) satisfies the Kurdyka-Łojasiewicz (K-Ł) property; ii) $\mathcal{L}_{k}$ is a coercive function. To prove point i), let us first introduce some definitions .

where $P_{1},\ldots,P_{k}$ and $Q_{1},\ldots,Q_{l}$ are polynomial in $n$ variables.

It is shown [ cf. , Th. $3$] that the semi-algebraic functions satisfy the K-Ł property.

A function \phi(\mbox{\boldmathx}) satisfies the Kurdyka-Łojasiewicz (K-Ł) property at point \bar{\mbox{\boldmathx}}\in\text{dom}(\partial\phi) if there exists $\theta\in[0,1)$ such that

is bounded around \bar{\mbox{\boldmathx}}.

The global convergence of the PAM method established in requires the objective function to satisfy the K-Ł property. Define \mbox{\mathbf{W}}:=(\mbox{\mathbf{X}},\mbox{\mathbf{P}}) and consider the function $\mathcal{L}_{k}$ in (IV-B), i.e.

where f_{1}(\mbox{\mathbf{X}})=\text{GDV}(\mbox{\mathbf{X}}), f_{2}(\mbox{\mathbf{P}})=\delta_{\mathcal{S}_{t}}(\mbox{\mathbf{P}}) and g_{k}(\mbox{\mathbf{X}},\mbox{\mathbf{P}})=\langle\mbox{\mathbf{\Lambda}}^{k},\mbox{\mathbf{P}}-\mbox{\mathbf{X}}\rangle+\frac{\rho^{k}}{2}\|\mbox{\mathbf{P}}-\mbox{\mathbf{X}}\|^{2}_{F}. Observe that f_{1}(\mbox{\mathbf{X}})=\displaystyle\sum_{i,j=1}^{N}a_{ji}\text{max}(x_{i}-x_{j},0) is the weighted sum of the functions $f_{ij}(x_{i},x_{j})=\text{max}(x_{i}-x_{j},0)$. Being a finite sum of semi-algebraic functions also a semi-algebraic function, it is sufficient to show that $f_{ij}$ is semi-algebraic. Assume, w.l.o.g. $y_{ij}=x_{i}-x_{j}$ so that $z=f_{ij}(y_{ij})=\text{max}(y_{ij},0)$. The graph of $f_{ij}$ becomes

and according to Definition 3 it is a semi-algebraic set. Then f_{1}(\mbox{\mathbf{X}}) as sum of semi-algebraic functions is also semi-algebraic. Since f_{2}(\mbox{\mathbf{P}}) and g_{k}(\mbox{\mathbf{X}},\mbox{\mathbf{P}}) are semi-algebraic functions it follows that \mathcal{L}_{k}(\mbox{\mathbf{W}}) is also semi-algebraic. It remains to prove point ii) to assess that $\mathcal{L}_{k}$ is a coercive function, i.e. \mathcal{L}_{k}(\mbox{\mathbf{W}})\rightarrow\infty when \|\mbox{\mathbf{W}}\|_{\infty}\rightarrow\infty. Clearly, the term f_{2}(\mbox{\mathbf{P}}) is coercive. The remaining terms in (45) can be written as

Since \mbox{\mathbf{P}}\in\mathcal{S}_{t} it holds \parallel\mbox{\mathbf{P}}\parallel_{F}^{2}=N. Thus, from the inequalities \langle\mbox{\mathbf{A}},\mbox{\mathbf{B}}\rangle\geq-\parallel\mbox{\mathbf{A}}\parallel_{F}\parallel\mbox{\mathbf{B}}\parallel_{F} and \parallel\mbox{\mathbf{B}}\parallel_{F}\leq\parallel\mbox{\mathbf{B}}\parallel_{1}, it holds \langle\mbox{\mathbf{\Lambda}}^{k},\mbox{\mathbf{P}}\rangle\geq-\sqrt{N}\parallel\mbox{\mathbf{\Lambda}}^{k}\parallel_{1}, so that one gets

References