New results on the local linear convergence of ADMM: a joint approach

TL;DR

This paper provides new theoretical insights into the local linear convergence of ADMM by analyzing the eigenstructure of key matrices, improving the understanding of its convergence speed and aiding parameter tuning.

Contribution

It introduces a joint approach to estimate the spectral radius of matrix products, enhancing the precision of convergence speed analysis compared to existing methods.

Findings

Develops a technique to locate eigenvalues of symmetric matrix products.

Strengthens the accuracy of ADMM convergence speed estimates.

Provides theoretical tools for better parameter tuning of ADMM.

Abstract

Thanks to its versatility, its simplicity, and its fast convergence, ADMM is among the most widely used approaches for solving a convex problem in distributed form. However, making it running efficiently is an art that requires a fine tuning of system parameters according to the specific application scenario, and which ultimately calls for a thorough understanding of the hidden mechanisms that control the convergence behaviour. In this framework we aim at providing new theoretical insights on the convergence process and specifically on some constituent matrices of ADMM whose eigenstructure provides a close link with the algorithm's convergence speed. One of the key technique that we develop allows to effectively locate the eigenvalues of a (symmetric) matrix product, thus being able to estimate the contraction properties of ADMM. In the comparison with the results available from the literature, we are able to strengthen the precision of our speed estimate thanks to the fact that we are solving a joint problem (i.e., we are identifying the spectral radius of the product of two matrices) in place of two separate problems (the product of two matrix norms).

Figures

Equations

$$min{f}_1(x_1)+f_2(x_2)$$

$$L(x,\lambda )$$

$$x_1^{+}$$

$${\mathcal{P}}_i(u)$$

$$x_1^{+}$$

$$z^{+}=E^{\frac{1}{2}}(y^{+}-b+\tilde{\lambda })$$

$$x_2$$

$$z^{+}=\mathcal{R}(z)=(1-q)z+q{\mathcal{D}}_1{\mathcal{D}}_2(z),$$

$$x_2$$

$$F_i=E^{\frac{1}{2}}{A}_i,M_i=F_i^T{F}_i=A_i^{T}E{A}_i,$$

$$q_M(u)$$

$${\mathcal{P}}_i(E^{-\frac{1}{2}}u)$$

$$g_i^{\ast }(u)=\underset{x_i}{\sup }⟨x_i,u⟩-g_i(x_i).$$

$$\partial g_i^{\ast }(u)=\left\{{\mathrm{argmax}}_{x_i}⟨x_i,u⟩-g_i(x_i)\right\},$$

$${\mathcal{P}}_i(E^{-\frac{1}{2}}u)\in \partial g_i^{\ast }(F_i^{T}u),$$

$${\mathcal{D}}_i(u)\in 2F_i\partial g_i^{\ast }(F_i^{T}u)-u\mp E^{\frac{1}{2}}b.$$

$$0⪯C_i+M_i⪯{\partial }^2{g}_i⪯S_i+M_i.$$

$$0⪯(S_i+M_i{)}^{-1}⪯{\partial }^2{g}_i^{\ast }⪯(C_i+M_i{)}^{-1},$$

$$-I⪯L_i⪯\partial D_i⪯U_i,$$

$$L_i$$

$${\tilde{F}}_i=F_i{M}_i^{-\frac{1}{2}},$$

$$L_i$$

$${\tilde{S}}_i$$

$${\ell }_{i,k}$$

$$h(x)=\frac{1-x}{1+x},$$

$$U_i⪯I$$

$$\Vert {\mathcal{D}}_i(u)-{\mathcal{D}}_i(v)\Vert \le {\mu }_i\Vert u-v\Vert ,$$

$${\mu }_i$$

$$v_{k,j}=\sqrt{ϵ}a{x}_k+\frac{m_j}{\sqrt{ϵ}a}$$

Full text

I Introduction

In recent years the alternating direction method of multipliers (ADMM) has received a considerable attention as an effective method to reach consensus among agents, to implement an optimisation algorithm in a distributed way, or, more generally, to solve networked problems by iteratively applying simple optimisation steps (e.g., see [1]). ADMM is in close relation with, and in many cases it is equivalent to, a wide number of alternative approaches, e.g., Douglas–Rachford splitting, proximal point algorithms, the split Bregman algorithm, that share the same potential, as well as the same techniques for analysis.

The reference problem in ADMM is

[TABLE]

where $f_{i}$ are proper, lower semicontinuous, and convex functions, $x_{i}$ are real valued vectors, and $A_{i}$ are linear operators (i.e., matrices). Solution to ((1)) is found by investigating the associated augmented Lagrangian function

[TABLE]

with $\lambda$ the Lagrange multipliers, $\langle\cdot,\cdot\rangle$ an inner product, and $\|u\|_{E}^{2}=\langle u,Eu\rangle$ a weighted norm that uses a positive definite augmentation matrix $E\succ 0$. For the sake of ease and practical applicability, the augmentation matrix $E$ is typically chosen a diagonal positive definite matrix, e.g., $E=\epsilon I$ for some augmentation factor $\epsilon>0$.

In the above context ADMM is an iterative algorithm that solves ((1)) by looking for saddle points of the augmented Lagrangian ((2)). This is achieved by performing an alternating direction search where each variable is separately updated by keeping the others fixed, that is,

[TABLE]

where the + sign denotes the updated variables. Depending on the choice of functions and linear constraints in ((1)), the above setting enables the possibility of solving a wide range of optimisation problems in distributed form.

Today ADMM is recognised to exhibit a linear (hence fast) convergence provided that parameters are wisely chosen. Unfortunately, parameters selection is a very difficult task. As a consequence, much of the recent effort in ADMM is devoted to deriving a deeper understanding of its convergence properties, as well as their dependence on system parameters. The literature on these aspects is vaste. The $O(1/n)$ convergence of ADMM, with $n$ denoting the iteration number, is addressed in [2, 3, 4, 5] under weak convexity assumptions. Linear convergence under strong convexity and/or Lipschitz continuity (smoothness) is studied in [6, 7, 8, 9, 10, 11, 12, 13], where the derivation is mainly carried out by building inequalities in the dual domain, and by exploiting the properties of proximity functions. A more general statement of linear convergence (under weaker conditions) is available from [12], and a statement on convergence in finite steps under a local polyhedral assumption is available in [13]. Convergence under specific scenarios is also investigated. Multi-block ADMM formalisations, where the target function is split in more than two contributions is available from [14, 15, 16, 12]. The linear convergence in a decentralised consensus optimisation problem, i.e., a specific problem that can be casted in the general ADMM setting to build a distributed algorithm, is proved in [17, 18, 19]. Convergence in the presence of inexact solvers is studied in [8], while a perspective on non-convex problems is available from [20, 21, 22].

Taken from the perspective of the usability of the result, a notable approach among the considered ones is given by the works of Giselsson and Boyd [10, 11], and of Liang, Fadili, and Peyré of [13], which we aim at generalising in some key aspects, and, up to a certain extent, also at merging and at simplifying. More specifically, in our work we provide a number of insights that are useful both on a practical as well as on a theoretical level, namely: we show that the Douglas–Ratchford structure is available in ADMM directly in the primal domain, and not only in the dual domain; we simplify derivations by avoiding side concepts (e.g., co-coercivity), and only relying on simple properties of the Legendre Fenchel transform, which, incidentally, also ensures a straightforward understanding of the convergence process in the presence of non-smooth functions; we relax [10, Assumption 2.ii], which limits the applicability of the results; we jointly consider both the functions of the alternate direction process and accordingly develop an efficient method to locate the eigenvalues of a (symmetric) matrix product; we use and generalise the findings of [13], thus being able to effectively estimate the convergence speed in the vicinity of the target point.

The paper is organised as follows. Section II provides the basis of our approach by showing that a link to a (generalised) Douglas–Rachford splitting is available directly in the primal domain, and that its properties are closely linked to a Legendre–Fenchel transform. These provide the essential structure of ADMM which is used in Section III to establish an equivalent to [10, Corollary 2] under a more general framework where a full-row-rank matrix assumption is dropped. Insights under the presence of discontinuities in the derivatives of functions $f_{i}$ are also provided. In Section IV we study the interdependencies between the two alternating steps in ((3)), and identify a specific matrix $R$ as the key feature of ADMM. The in-depth study of the eigenvalue structure of $R$, which controls the convergence process, is given in Section V, while its (long) technical proof is available in the Appendix. A brief example of application is given Section VI. Section VII concludes the paper.

II Essential structure of ADMM

II-A Douglas–Rachford splitting structure in the primal domain

The starting point is to reveal an alternative representation of ADMM in the form of a (generalised) Douglas–Rachford splitting, which is particularly well suited for studying the convergence behaviour. The Douglas–Rachford representation is a standard result of the ADMM literature, commonly adopted in the dual domain. Here, instead, we reveal it directly in the primal domain. To this aim, in the following we use a compact notation where

[TABLE]

are, respectively, a proximity operator and a reflected proximity operator for function $f_{i}$, $i\in\{1,2\}$. With the above notation, the ADMM updates ((3)) can be written as

[TABLE]

where $\tilde{\lambda}=E^{-1}\lambda$. This reveals that the update of $x_{i}$ is constrained by the proximity operators $\mathcal{P}_{i}$, which ensures a form of coordination among the alternate minimisation steps.

Mimicking the most efficient solutions available in the literature, we further generalise ((5)) in a so called ADMM with scaled variables where the updates are modified in the form

[TABLE]

for some real constant $q\neq 0$. Note that in ((6)) the variable $y$ plays the role of an estimate of $A_{1}x_{1}$ (or, equivalently, of $A_{2}x_{2}+b$), exploiting both $x_{1}$ and $x_{2}$ in a (possibly) more efficient way that linearly combines both the directions of the alternate search. The above corresponds to the standard ADMM algorithm when $q=\frac{1}{2}$.

The Douglas–Rachford structure becomes evident by introducing variable

[TABLE]

in place of variable $\tilde{\lambda}$. Note that, although the presence of $E^{\frac{1}{2}}$ could be dropped, this will guarantee a useful form of symmetry in later derivations, and therefore we keep it despite the redundant notation. Hence, by replacing $\tilde{\lambda}$ with $z^{+}$, the iterative algorithm ((6)) can be expressed in the equivalent form111Note that: the first of ((8)) is obtained by simply replacing ((7)) in the third of ((6)); the third of ((8)) is obtained from the fourth of ((6)) by replacing $\tilde{\lambda}$ according to ((7)), and then $y^{+}$ according to the second of ((6)); the second of ((8)) is obtained from the first of ((6)) by replacing $\tilde{\lambda}$ according to ((7)), then $z^{+}$ according to the third of ((8)), and finally $y^{+}$ according to the second of ((6)).

[TABLE]

showing that $z$ plays the role of a state variable in the ADMM update. Furthermore, by replacing the updates of $x_{1}$ and $x_{2}$ in the third of ((8)), then the update of $z$ can be more compactly expressed through the (generalized) Douglas–Rachford splitting relation

[TABLE]

where $\mathcal{D}_{i}$ are the reflected proximity operators defined in ((4)). We also have

[TABLE]

which express the link between $z$ and the target variables that we are estimating. Note that ((9)) evidences the essential structure of ADMM directly in the primal domain (e.g., note the closeness with its dual counterpart [10, eq. (24)]), and shows that convergence properties are closely linked to the contraction properties of operator $\mathcal{D}_{1}\mathcal{D}_{2}$.

Formally, the three approaches ((5)), ((6)), and ((8)) share the same stationary points, in the sense that there is a one-to-one correspondence between their stationary points. The relation is given by the mappings $y^{*}=A_{1}x_{1}^{*}=A_{2}x_{2}^{*}+b$ and $z^{*}=E^{\frac{1}{2}}(y^{*}-b+E^{-1}\lambda^{*})$, with ∗ denoting the stationary values. In addition, these stationary points coincide with those of the Lagrangian ((2)), i.e., they solve the minimisation problem. Evidently, we are implicitly assuming that such stationary points (i.e., a solution) exist, that is throughout this paper we are taking the following

Assumption 1

A solution to ((1)) exists, i.e., equation ((9)) has at least one stationary point.  $\Box$

II-B Legendre–Fenchel structure of $\mathcal{D}_{i}$

According to ((9)), the convergence behaviour of the ADMM algorithm can be understood (and optimised) by investigating the properties of operators $\mathcal{D}_{i}$. In this respect, it is appropriate to give $\mathcal{D}_{i}$ a neater form, given by a Legendre–Fenchel transform, which better reveals the dependence on (the properties of) functions $f_{i}$. The techniques used in the following are taken from the state-of-the-art literature on ADMM, and wisely exploited.

The starting point for inspecting operator $\mathcal{D}_{i}$ is to first interpret the structure of the proximal operator $\mathcal{P}_{i}$. To do so we introduce the shorthand notations

[TABLE]

and

[TABLE]

which allow writing

[TABLE]

In oder to build upon ((13)), an important requirement is given by

Assumption 2

Functions $g_{i}$ are proper, lower semicontinuous, and convex. $\Box$

This is naturally satisfied if, as discussed in the Introduction, functions $f_{i}$ are proper, lower semicontinuous, and convex. However, Assumption 2 is more general, and further considers the possibility that $f_{i}$ are non-convex, provided that the structure of matrices $M_{i}$ is such to make the sum $f_{i}+q_{M_{i}}$ convex.

We then observe that the second line of ((13)) is related to the Legendre–Fenchel transform $g_{i}^{*}$ of $g_{i}$, defined as (e.g., see [23, §11])

[TABLE]

The explicit link between ((14)) and ((13)) is the standard result [23, Proposition 11.3]

[TABLE]

where $\partial g_{i}^{*}$ expresses the sub-gradient of $g_{i}^{*}$, and the curly brackets identify the sets of arguments that maximise their target function. The validity of ((15)) is guaranteed by Assumption 2. Therefore, the explicit link between $g_{i}^{*}$ and the proximity operator is, by inspection,

[TABLE]

which further provides

[TABLE]

Equation ((17)) shows that the contraction properties of operators $\mathcal{D}_{i}$ are (linearly) linked to those of the Legendre–Fenchel transform $g_{i}^{*}$. Hence, an understanding of the convergence behaviour of ADMM needs investigating $\partial g_{i}^{*}$.

III Contraction properties of the reflected proximity operators $\mathcal{D}_{i}$

In the following we discuss the properties of ((17)) by first concentrating on a simplified context where functions $g_{i}$ are convex and smooth, which essentially leads to the outcomes of [10], although under a more general framework. We then generalise the result to a more realistic context where functions are non-smooth, e.g., in the presence of piecewise functions or boundaries, which clarifies the added value of the proposed approach.

III-A Smooth functions scenario

We first take the simple (but effective) scenario used in [17, 10], namely

Assumption 3

Function $f_{i}$ is $S_{i}$-smooth and $C_{i}$-strongly convex, meaning that both functions $f_{i}-q_{C_{i}}$ and $q_{S_{i}}-f_{i}$ are convex, where $S_{i}$ and $C_{i}$ are symmetric matrices satisfying $S_{i}\succeq C_{i}$, and where $\succeq$ is the positive semidefinite ordering operator. $\Box$

According to Assumption 2, this corresponds to a Lipschitz property that limits the Hessian of $g_{i}$ to the range (see also the definition of $g_{i}$ in ((12)))

[TABLE]

As a consequence of [24, Theorem 18.15], the above reveals a Lipschitz property in the Legendre–Fenchel domain, namely the property

[TABLE]

where the negative exponent stands for a pseudo–inverse operator when matrices are non–invertible. This, in connection with ((17)), evidences that the sub-gradient of operator $\mathcal{D}_{i}$ satisfies the Lipschitz property

[TABLE]

where $I$ is the identity matrix, and where

[TABLE]

are lower and upper bounds. Incidentally, the above (as well as what follows) also holds in the limit for matrices $C_{i}+M_{i}$ and $S_{i}+M_{i}$ having eigenvalues which are arbitrarily small (i.e., $\rightarrow 0$) or arbitrarily large (i.e., $\rightarrow\infty$).

A neater expression for bounds ((21)) can be obtained by introducing matrices

[TABLE]

providing222This is a straightforward consequence of standard properties of the Moore–Penrose pseudo–inverse $F_{i}^{\dagger}=M_{i}^{-1}F_{i}^{T}$.

[TABLE]

with symmetric matrices

[TABLE]

satisfying $\tilde{S}_{i}\succeq\tilde{C}_{i}\succeq-I$ (see also ((18))). Note that $L_{i}$, $U_{i}$, $\tilde{S}_{i}$, and $\tilde{C}_{i}$ are symmetric real matrices that are $\succeq-I$, and their eigenvalues are, therefore, real valued and greater than or equal to $-1$.

The structure revealed in ((23)) is common to a number of findings available in the literature [14, 8, 10]. However, a careful reader can appreciate that we leverage the same essential properties governing ADMM, but the present derivation avoids the difficulties involved with a study performed in the dual domain, and simplifies the proofs. In addition, the present formalisation contains some generalisations that can be fully appreciated by investigating the eigen-structure of matrices. Specifically, from standard linear algebra considerations it is easy to see that the eigenvalues of lower and upper bound matrices ((23)) are given, respectively, by

[TABLE]

where (see a graphical representation in Fig. 1

)

[TABLE]

and where $\tilde{s}_{i,k}$ and $\tilde{c}_{i,k}$ are those (real valued) eigenvalues of matrices $\tilde{S}_{i}$ and $\tilde{C}_{i}$, respectively, that are associated with the span of matrix $M_{i}$, and are therefore in number of ${\rm rank}(M_{i})={\rm rank}(A_{i})$. Incidentally, the eigenvalues $\tilde{s}_{i,k}$ and $\tilde{c}_{i,k}$ can be extracted, respectively, from matrices $S_{i}M_{i}^{-1}$ and $C_{i}M_{i}^{-1}$ which avoid the square–root operation and are diagonalisable by construction. All the eigenvalues in ((25)) are real valued.

It is important to observe that the eigenvalues $-1$ in ((25)) are a consequence of the presence in ((23)) of the contribution $I-\tilde{F}_{i}\tilde{F}_{i}^{T}$ (if active), which is the projector associated to the kernel of $F_{i}$. These values are the most significant novelty we are providing, and are also a straightforward consequence of the fact that we are avoiding [10, Assumption 2.ii], i.e., the request of $A_{i}$ being a full-row-rank matrix. Incidentally, the presence of eigenvalues $-1$ is common in the very many optimisation problems where some of the variables are duplicated more than once.

III-B Smoothness + convexity scenario

If we focus on a scenario where functions $f_{i}$ are convex, which corresponds to setting $C_{i}\succeq 0$, then the natural consequence is that eigenvalues $\tilde{s}_{i,k}$ and $\tilde{c}_{i,k}$ are non negative. Hence, by inspection of the structural function $h(x)$ for $x\geq 0$ (see Fig. 1), all eigenvalues $\ell_{i,k}$ and $\nu_{i,k}$ as defined in ((25)) are constrained in the range $[-1,1]$. The straightforward outcome is that

[TABLE]

holds, which further guarantees that $-I\preceq\partial D_{i}\preceq I$, and, in turn, ensures the continuity of functionals $\mathcal{D}_{i}$ as well as their non expansiveness.

If we further denote with $\ell_{i,\min}\geq-1$ the minimum eigenvalue of $L_{i}$, and with $\nu_{i,\max}\leq 1$ the maximum eigenvalue of $U_{i}$, we also have $I\cdot\ell_{i,\min}\preceq\partial D_{i}\preceq I\cdot\nu_{i,\max}$, whose specific consequence is the contractive property

[TABLE]

with contraction parameter

[TABLE]

satisfying $0<\mu_{i}\leq 1$. Under the assumption of a full-row-ranked matrix $A_{i}$, and although we are using a different notation, this result perfectly corresponds to [10, Proposition 5] (the more general result of the paper).

III-C Non-smooth functions scenario

In very many applications the functions $f_{i}$ involved in ADMM do not satisfy the smoothness property, either because they are non smooth by construction, or because they have a limited domain. In these situations, the formalisation identified by ((17)) is able to fully (and simply) capture the essence of the problem by a simple calculation of a Legendre-Fenchel transform. To make our point, we investigate a context where variables are simply duplicated, and functions $f_{i}$ are expressed as a linear combination of one–dimensional convex contributions. This is equivalent to considering that $E=\epsilon$ and $A_{i}=a$ are scalar values, that $b=0$ (for simplicity, since the value of $b$ is irrelevant for the contracting properties of $\mathcal{D}_{i}$), and that functions $f_{i}(x)$ are one-dimensional. This simplified scenario is able to capture the essence of the most general case.

To begin with we assume a piecewise linear form for $f_{i}$, thus obtaining a tilted-staircase reflected proximity operator $\mathcal{D}_{i}$ as illustrated in Fig. 2

. The expression/structure of $\mathcal{D}_{i}$ was derived by application of ((17)), and by exploiting the identity $\partial g_{i}^{*}=(\partial g_{i})^{-1}$, which is valid thanks to the convexity of $g_{i}$ [23, Proposition 11.3]. The relevant facts that can be observed in Fig. 2 are:

1. a)

the operator $\mathcal{D}_{i}$ has a (rotated) staircase structure where a negative $-1$ slope alternates with a positive $+1$ slope, negative slopes being associated to breaking points $x_{k}$ and positive slopes to the constant derivatives $m_{j}$, the general expression of lines being available in Fig. 2; observe that in a multidimensional (generalised) context, this would map into the fact that the domain is partitioned in regions where slopes can either be $-1$ or $+1$, where $1$ is intended as the all-ones vector; 2. b)

the larger the augmentation variable $\epsilon$ the wider the positive–sloped segments, while the smaller $\epsilon$ the wider the negative–sloped segments, a result which is a consequence of the value

[TABLE]

of the coordinate where a negative–sloped ($x_{k}$) and a positive–sloped ($m_{j}$) curve intersect (see the exemplifying $v_{4,4}$ and $v_{5,4}$ in Fig. 2); 3. c)

the active range $[\ell_{i,\min},\nu_{i,\max}]\subset[-1,1]$ (see Section III.B) may significantly change over different intervals, e.g., in the depicted case where in the region around $v$ we have $\nu_{i,\max}=1$ but $\ell_{i,\min}>-1$; as a natural consequence, the convergence behaviour strongly depends on the limit point.

The above provides an operative method to study convergence bounds $[\ell_{i,\min},\nu_{i,\max}]$ in the vicinity of a target point, e.g., in the vicinity of the target $z^{*}$, when this is known or it can be calculated with some precision.

The findings of Fig. 2 can be generalised to build a simple procedure that identifies the proximity operator $\mathcal{D}_{i}$ when functions $f_{i}$ are non-smooth, but not necessarily piecewise-linear. To do so, observe that, according to Fig. 2, the value $v_{k,j}$ in ((30)) is mapped by the reflected proximity operator into

[TABLE]

If a piecewise–linear curve like the one in Fig. 2 is taken to approximate, in the limit, any given function $f_{i}$, then, by virtue of continuity, we can interpret value $m_{j}$ in ((30)) and ((31)) as $m_{j}=f^{\prime}_{i}(x_{j})$ and we are therefore able to capture the shape of function $\mathcal{D}_{i}(u)$ through the couple of functions

[TABLE]

as well as its derivative through

[TABLE]

This identifies a simple and effective method to practically infer the structure of $\mathcal{D}_{i}$ when $f_{i}$ is a piecewise function, where:

((32)) and ((33)) provide the shape of $\mathcal{D}_{i}$ associated with the smooth pieces of a piecewise function $f_{i}$; 2. 2.

the link between such pieces is captured, for each breakpoint $x_{k}$, by a negative–sloped contribution of the form $2\sqrt{\epsilon}ax_{k}-u$ (see Fig. 2).

Although derived for a convex function, the above method is perfectly applicable to a non–convex setting, provided that Assumption 2 holds, i.e., that $u(x)$ is an increasing function.

IV Joint role of the reflected proximity operators $\mathcal{D}_{1}\mathcal{D}_{2}$

IV-A Merging the reflected proximity operators

The fact that, in the reference equation ((9)), the reflected proximity operators appear in the form $\mathcal{D}_{1}\mathcal{D}_{2}$ suggests an analysis of their joint role. To do so we assume that matrices $L_{i}$ and $U_{i}$ commute, which allows writing the bounds in ((20)) in the form

[TABLE]

for some orthogonal (real–valued and unitary) matrix $V_{i}$, where $\ell_{i}$ and $\nu_{i}$ collect the eigenvalues ((25)), and where, thanks to assumption $L_{i}\preceq U_{i}$, the eigenvalues satisfy the ordering

[TABLE]

We are therefore taking the following:

Assumption 4

Matrices $\tilde{S}_{i}$ and $\tilde{C}_{i}$ in ((24)) commute, or, equivalently, $S_{i}M_{i}^{-1}C_{i}=C_{i}M_{i}^{-1}S_{i}$.  $\Box$

Remark 1

The assumption is satisfied in very many practical cases. Furthermore, it can be always ensured by simply relaxing the true (stricter) bounds declared by $C_{i}$ and $S_{i}$. In practical setups this is straightforwardly obtained, e.g., by making sure by construction that $M_{i}$, $C_{i}$, and $S_{i}$ are diagonal matrices.  $\Box$

Remark 2

In the most general case where $L_{i}$ and $U_{i}$ do not commute, ((34)) holds in the form

[TABLE]

where $V_{i}$ is real–valued and invertible, but not necessarily unitary, and where $\breve{\ell}_{i}$ and $\breve{\nu_{i}}$ are vectors identifying the common diagonal structure between matrices $L_{i}$ and $U_{i}$. This can be proved, e.g., by exploiting the fact that two positive–semidefinite matrices can be simultaneously diagonalised by an invertible matrix (e.g., see [25]). Thanks to the Inertia Theorem, $\breve{\ell}_{i}$ and $\breve{u}_{i}$ keep the same signs of $\ell_{i}$ and $u_{i}$, and some relations on ranges can also be inferred. However, the fact that $V_{i}$ is non–unitary makes some of the derivations of this paper hard, this being the motivation for taking Assumption 4.  $\Box$

The above allows to characterise the Lipschitz properties of $\mathcal{D}_{i}(u)$. Specifically, it is easy to verify that the operator

[TABLE]

satisfies the Lipschitz property333This is a consequence of the fact that $\partial\tilde{\mathcal{D}}_{i}=V_{i}\partial\mathcal{D}_{i}V_{i}^{T}$. $\mathrm{diag}({\ell}_{i})\preceq\partial\tilde{\mathcal{D}}_{i}\preceq\mathrm{diag}({\nu}_{i})$. This corresponds to verifying that the relation

[TABLE]

holds for some vector $\alpha_{i}\in[{\ell}_{i},{\nu}_{i}]$, which is a straightforward consequence of observing (along each direction) an upper and a lower bound on the derivatives. By mapping back onto $\mathcal{D}_{i}$ we further have

[TABLE]

for the same vector $\alpha_{i}$. The use of ((39)) in ((9)) finally provides our target result

[TABLE]

for some $\alpha=[\alpha_{1},\alpha_{2}]$ and $\alpha_{i}\in[{\ell}_{i},{\nu}_{i}]$. Note that, although expressed with a different notation, the structure of $R(\alpha)$ in ((40)) corresponds to [13, Eq. (16)], with the added value that $\alpha$ is not necessarily capturing the behaviour only in the strict vicinity of the target point $z^{*}$, but rather on a larger span. This corresponds to the fact that $\alpha$ belongs to a range as opposed to taking a specific value.

IV-B Global contraction parameter

The above analysis suggests that the contraction parameter of global interest, $\mu$, can be identified as the solution to the non-convex problem

[TABLE]

where the matrix norm $\|\cdot\|$ identifies the largest singular value of $R(\alpha)$. By applying this result to ((40)) we obtain

[TABLE]

so that, by respectively choosing $(u,v)=(z,z^{-})$ and $(u,v)=(z,z^{*})$, with $z^{*}$ the target value, this establishes that

[TABLE]

which corresponds to a contractive property of the ADMM update when $\mu<1$.

Note that this generalises the analysis of [10] where contraction properties ((28))-((29)) are independently exploited to reveal that

[TABLE]

and it is easy to verify that $\tilde{\mu}\geq\mu$ thanks to the fact that we are jointly considering the role of the reflected proximity operators. We are also generalising (to ADMM) the results available for a Douglas–Rachford splitting approach in [11], and in fact we are able to account for the presence of matrices $A_{i}$, as well as to deal with a non-elementary matrix structure for $S_{i}$ and $C_{i}$. The following example shows the quality of the proposed methodology through an improvement over available results, even in a very simple case.

Example

We consider the problem of [11, Sect. 5], that is a convex scenario with $C_{1}=0$, $S_{1}=\beta I\succeq 0$, $C_{2}=\sigma I\succeq 0$, $S_{2}=\infty I$, $A_{1}=A_{2}=I$, and $E=\gamma I$. Use in ((23)) provides

[TABLE]

so that $V_{1}=V_{2}=I$, and in ((40)) we have $N(\alpha)=\alpha_{1}\alpha_{2}I$ with $h({\textstyle\frac{\beta}{\gamma}})\leq\alpha_{1}\leq 1$ and $-1\leq\alpha_{2}\leq h({\textstyle\frac{\sigma}{\gamma}})$. As a consequence, the contraction parameter ((41)) becomes

[TABLE]

One can easily see that the effective range of the product $\alpha_{1}\alpha_{2}$ is captured by the expression

[TABLE]

which can be used to identify the optimum $\gamma$, i.e., the one providing the smallest range. With a little effort444We skip the derivation since it is lengthy but trivial. As a matter of fact it simply implies testing all the possible signs of $h({\textstyle\frac{\beta}{\gamma}})$ and $h({\textstyle\frac{\sigma}{\gamma}})$, along with the constraints $\sigma<\beta$ or $\sigma>\beta$. we find that this optimal parameter is $\gamma^{*}=\sqrt{\sigma\beta}$ providing

[TABLE]

Use of this result in ((46)) further ensures that the minimum contraction parameter is $\mu^{*}=2q^{*}-1$ for $q^{*}=2/(3-H_{1}^{*})$. Hence, by substitution, the optimal contraction parameter in the considered setting is

[TABLE]

We are therefore able to strengthen the result of [11, Theorem 5.6], which identifies the same optimal parameter $\gamma^{*}$, but a weaker bound $(1+\sqrt{\sigma/\beta})^{-1}>\mu^{*}$. The difference is depicted in Fig. 3.

$\Box$

V Characterisation of the local convergence through eigenvalues

Following [13], an alternative to ((41)) is to investigate the measure

[TABLE]

where $\rho(\cdot)$ identifies the the spectral radius of $R(\alpha)$. The underlying idea is that, in the vicinity of the convergence point $R(\alpha)$ becomes stationary, and therefore the most relevant parameter that determines the convergence speed is the spectral radius,555The role of the spectral radius is a consequence of the fact that $\rho(A)=\lim_{k\rightarrow\infty}\|A^{k}\|^{1/k}$ (e.g., see [26]). which satisfies $\rho(\cdot)\leq\|\cdot\|$ and hence is a more precise indicator than the matrix norm.

Note that ((41)) is providing a global information on the convergence speed, while ((48)) is providing a local information in the vicinity of the target point. Hence, both are valuable for understanding how to possibly control the convergence process. However, the spectral radius perspective is particularly welcome thanks to the fact that we are able to effectively locate the eigenvalues of $R(\alpha)$ starting from a loose indication given by the range constraint $\alpha_{i}\in[{\ell}_{i},{\nu}_{i}]$. We show this by first concentrating on a simpler scenario where $\alpha_{i}$ can take two values only, which captures the essence of the problem, and by then generalising the result to a more conventional setting.

V-A A simplified scenario

We begin by investigating the span of the eigenvalues and the singular values of $R(\alpha)$ in a reference scenario where each $\alpha_{i}$ can take either a positive or a negative value (one of which can be zero). This scenario is, e.g., active at convergence for a linear problem, in which case the limit point $z^{*}$ lies somewhere on a staircase function of the form of Fig. 2.a and assumes values $1$ or $-1$. The more general framework that we propose, however, is one where levels can be different from $\pm 1$, which generalises the results available in [13] (see also [27, 28, 29]).

The reference tool for obtaining our results is the CS decomposition of [30, Theorem 2.5.3] (see details in Appendix A), and provides the following

Theorem 1

Assume that the entries of $\alpha$ in ((40)) are ordered and take either a positive or a negative value, that is

[TABLE]

for some positive constants $\overline{p}_{i}\geq 0$ and $\overline{n}_{i}\geq 0$, and where $p_{i}$ and $n_{i}$ denote, respectively, the number of positive and negative entries. Accordingly to the above, consider expressing the matrix product $V_{1}V_{2}^{T}$ in the form

[TABLE]

where $G_{1}$ is $p_{2}\times p_{1}$, $G_{2}$ is $p_{2}\times n_{1}$, $G_{3}$ is $n_{2}\times p_{1}$, and $G_{4}$ is $n_{2}\times n_{1}$. Denote with $c_{1},\ldots,c_{p}$ the singular values of either $G_{1}$ or $G_{4}$, where $G_{1}$ is used if $\min(p_{1},p_{2})\leq\min(n_{1},n_{2})$ and $G_{4}$ is used otherwise. The singular values are in number of $p=\min(p_{1},p_{2},n_{1},n_{2})$, and belong to the range $c_{i}\in[0,1]$. Define constants

[TABLE]

Then the eigenvalues of matrix $N(\alpha)$ are known in the closed form, and we have that:

each singular value $c_{i}$ generates a couple of distinct eigenvalues by the rule

[TABLE]

for a total of $2p$ eigenvalues; 2. 2.

the remaining eigenvalues are real–valued, and can be positive with value $\overline{p}$ or $\breve{p}$ (either of the two, but not both), or negative with value $-\overline{n}$ or $-\breve{n}$ (either of the two, but not both).

$\Box$**

Proof:

See Appendix A. ∎

We observe that, depending on the value of $c_{i}$, the eigenvalues ((52)) can be real and positive, real and negative, or complex valued. The transition is controlled by the threshold values

[TABLE]

which satisfy $\underline{c}\leq\overline{c}$, and which guarantee that the eigenvalues ((52)) satisfy the following:

1. a)

for $c_{i}<\underline{c}$ they are real valued and negative in the range $[-\overline{n},-\breve{n}]$; 2. b)

for $c_{i}>\overline{c}$ they are real valued and positive in the range $[\breve{p},\overline{p}]$; 3. c)

for $c_{i}\in[\underline{c},\overline{c}]$ they form a complex–conjugate couple with absolute value $\overline{r}$.

According to this result, a pictorial representation of the locus of the eigenvalues of $N(\alpha)$, namely the set

[TABLE]

is shown in Fig. 4

.

Incidentally, the above characterisation is complete for an understanding of the spectral behaviour of matrix $R(\alpha)$, i.e., for identifying ((48)). As a matter of fact, the locus ${\mathbb N}$ is mapped in the eigenvalue set of $R(\alpha)$ by a simple linear map, namely

[TABLE]

whose action is exemplified in Fig. 5 for:

a convex problem (Fig. 5, above) where ((27)) ensures that ${\mathbb N}$ is contained in the unit circle, and showing how the map ((55)) can effectively reduce the spectral radius; we used $0<q<1$, which corresponds to a contraction in the complex plane towards the fixed point $1$ (the point is named fixed since the map ${\mathbb N}\rightarrow{\mathbb R}$ does not change its position, this being true for any choice of $q$); 2. 2.

a non–convex problem (Fig. 5, below) where the action of $q$ is unable to map inside the unit circle those points that, broadly speaking, lie on the right of the fixed point $1$.

Note that the above clearly explains the role of $q$, that is, it controls the possibility of shifting the eigenvalues, possibly inside the unit circle and towards [math] (in such a way to increase the convergence speed). It also suggests that different eigenvalues would need in principle a different value of $q$ since, e.g., $\lambda=0.95$ calls for a value $q>1$ to be shifted towards zero, while for a similar effect $\lambda=-0.95$ calls for $0<q<1$, and $\lambda=1.05$ calls for $q<-1$.

From a graphical inspection of Fig. 5 we can also observe that, irrespective of the value of $q$, the spectral radius is determined by the values of ${\mathbb R}$ lying on the real axis, the reasons being that ((55)) maps the circle of radius $\overline{r}$ in a circle of radius $|q|\overline{r}$ centred in $1-q$. This provides the following results whose (simple) proof is left to the reader:

Corollary 2

Under the assumptions of Theorem 1, the spectral radius ((48)) is upper bounded by

[TABLE]

The resulting optimum choice $q_{\rho}^{*}$ minimising ((56)) is therefore

[TABLE]

and the related bound is

[TABLE]

$\Box$**

Proof:

Left to the reader. ∎

Interestingly, for $0<\overline{p}<1$ the bound ((58)) satisfies $\rho_{\max}^{*}<1$, in which case we are guaranteeing a locally linear convergence to the target point. When, instead, $\overline{p}>1$, then the above formalisation is suggesting the inappropriate choice $q=0$, that is, we are not able to build a convergent algorithm.

V-B A more general statement

Although the four–level assumption of ((49)) is linked to a very specific setting observed at convergence, its outcomes are more general than one can expect. As a matter of fact, the structures identified so far (i.e., the shape of Fig. 4) are observed also when the values of $\alpha$ belong to a range. This leads to our general statement, providing a description of the convergence process with a larger view.

Theorem 3

Assume that the entries of $\alpha_{i}$ in ((40)) are limited to the range

[TABLE]

with $0\leq\underline{n}_{i}\leq\overline{n}_{i}$ and $0\leq\underline{p}_{i}\leq\overline{p}_{i}$, and define the positive constants

[TABLE]

with $0\leq\underline{n}\leq\underline{r}\leq\overline{r}\leq\overline{n}$ and $0\leq\underline{p}\leq\underline{r}\leq\overline{r}\leq\overline{p}$. Then, the real–valued eigenvalues of $N(\alpha)$ belong to the range

[TABLE]

while the complex–valued eigenvalues satisfy

[TABLE]

$\Box$**

Proof:

See Appendix B. ∎

The outcomes of Theorem 3 are illustrated in Fig. 6

, where

[TABLE]

The implications on $R(\alpha)$ are, evidently, identical to the simplified case, that is, both ((55)) and Corollary 2 hold also in the general case. Notably, the structure of Fig. 6 is recurring in many fields of study (e.g., in network science [31] or in the characterisation of large random matrices), the reason being that with Theorem 3 we are providing a characterisation of the eigenvalue span of any real matrix since any real matrix can be expressed as the product of two symmetric matrices [26].

The interesting outcome of our result is that by identifying the parameters $\overline{n}_{i},\underline{n}_{i},\underline{p}_{i},\overline{p}_{i}$ in both the alternate directions (whose identification usually corresponds to a very simple task), we are able to characterise the location of the eigenvalues of $R(\alpha)$ without the need to calculate them explicitly, the added value being that calculating the eigenvalues of a large matrix might be a costly operation.

VI Application example

For illustration purposes we consider the simple weighted Lasso minimisation problem used in [10]

[TABLE]

where $\|\cdot\|_{p}$ is the Euclidean norm of order $p$, $\circ$ denotes the entry–wise Hadamard product, $x$, $o$, and $w$ are real–valued vectors of length, respectively, $200$, $300$, and $200$, where $\Omega$ is a sparse $300\times 200$ data matrix with an average of $10$ active entries per row, and where the (active) entries of $\Omega$ and $\omega$ are drawn from a Gaussian distribution with zero mean and unit variance, while the entries of $w$ are taken from a uniform distribution in $[0,1]$.

The problem is solved through ADMM by setting $x_{1}=x_{2}=x$, hence $A_{1}=A_{2}=I$ and $b=0$. We also assume a diagonal augmentation matrix $E=\epsilon I$. This provides a linear reflected proximity operator $\mathcal{D}_{1}$ of the form

[TABLE]

where $\Omega^{T}\Omega\,\epsilon^{-1}=\tilde{S}_{1}=\tilde{C}_{1}$, so that in ((25)) we have $\ell_{1,k}=\nu_{1,k}=h(\lambda_{k}(\Omega^{T}\Omega)/\epsilon)$. This corresponds to a range of the form $\alpha_{1,j}\in[h(\lambda_{\max}/\epsilon),h(\lambda_{\min}/\epsilon)]$ where

[TABLE]

are the maximum and minimum eigenvalue, respectively, of matrix $\Omega^{T}\Omega$. As a straightforward consequence, we can set

[TABLE]

where $[u]^{+}$ is the plus operator providing $u$ for $u>0$ and [math] otherwise. For the reflected proximity operator $\mathcal{D}_{2}$ we can instead exploit the technique developed in Fig. 2, to derive its $j$th entry $\mathcal{D}_{2,j}(u_{j})$ as illustrated in Fig. 7

. The staircase structure suggests that

[TABLE]

on the entire iteration span, but $\underline{p}_{2}$ and $\underline{n}_{2}$ are ensured to be $\simeq 1$ in the vicinity of the target point. By using the weaker (but simpler) constraint ((68)), we obtain

[TABLE]

and $\underline{p}=\underline{n}=\underline{r}=0$, whose symmetry, according to ((57)), identifies the optimum parameter $q=1$.

The study of the algorithm convergence when using $q=1$ and $\epsilon=1$ is shown in Fig. 8

and Fig. 9

for one specific realisation of the random matrices and vectors $\Omega$, $o$, and $w$. In the specific example we have $\lambda_{\max}=60.75$ and $\lambda_{\min}=0.3465$, providing $\overline{p}_{1}=0.4853$ and $\overline{n}_{1}=0.9676$. Note in Fig. 8 the correspondence between the locus of eigenvalues ${\mathbb R}$ (in gray) and the eigenvalues at the limit point $z^{*}$ (dots). The locus ${\mathbb R}$ was identified according to Theorem 3 and ((69)), and perfectly captures the position of the eigenvalues, especially in the real axis where the eigenvalue with the largest absolute value $\lambda^{*}_{\max}=0.9543$ is very close to the maximum allowed value $\overline{p}=0.9676$. Note that the bound given by ${\mathbb R}$ is in any case much stricter in the real axis than in the complex plane (circle $\overline{r}$), this being a common behaviour for randomly wired networks such as the one we are considering. The convergence plot of Fig. 8 further evidences that, although $\lambda^{*}_{\max}$ provides a precise indication of the convergence speed at the limit point, the value of $\rho_{\max}$ identified according to Theorem 3 and ((69)), namely $\rho_{\max}=0.9676$ in the considered case, is an upper value that is often able to capture with greater accuracy the true speed of the convergence process. In addition, $\rho_{\max}$ is easier to calculate since it does not rely on the (precise) knowledge of the limit point $z^{*}$, which in a practical application setup might be hardly known or derived.

VII Conclusions

In this paper we provided new theoretical insights on the convergence process of ADMM, which were derived by interpreting ADMM as a (generalised) Douglas-Rachford splitting in the time domain (as opposed to the classical frequency domain interpretation). The interesting outcome is that convergence speed is linked to the eigenvalues of matrices $R(\alpha)$ whose general structure is evidenced in Fig. 6, where a constraint on the absolute value of complex eigenvalues cohabits with a (looser) constraint on the range of real-valued eigenvalues. The latter constraint is what will ultimately set the convergence speed, and calls for novel speeding methods. For example, under the assumption that $\partial\mathcal{D}_{i}$ is constant, one might consider a cyclostationary $q$ where different values $q_{k}$ (cyclically repeating) are used at each iteration in order to selectively trim the effect of the larger (real-valued) eigenvalues. The solution in a more general and realistic setting is a tough problem which calls for future investigations and novel ideas.

Appendix A Proof of Theorem 1

In the following we will denote with $A\sim B$ a similarity between matrices $A$ and $B$ in the sense that they share the same eigenvalues with the same multiplicities. We will also use the property $AB\sim BA$. The starting point is to identify the similarity (see ((40)))

[TABLE]

where we used the notation ((49)) and ((50)). Note that matrices $H_{i}$ are real valued and diagonal, and matrix $G$ is orthogonal (real valued and unitary).

We then want to identify an operative structure for matrix $G$. To do so we initially assume that

[TABLE]

The general case, which is a simple extension of the outcomes of ((71)), will be discussed at the end of the proof. Hence, we can readily exploit the CS decomposition of [30, Theorem 2.5.3], which we write in the following form:

Proposition 4

Under condition ((71)), matrix $G$ can be expressed as the product

[TABLE]

where $A_{1}$, $A_{2}$, $B_{1}$, and $B_{2}$ are orthogonal matrices of dimension, respectively, $p_{2}\times p_{2}$, $n_{2}\times n_{2}$, $p_{1}\times p_{1}$, $n_{1}\times n_{1}$, and where $M$ is the orthogonal matrix

[TABLE]

where $q_{1}=p_{1}-p_{2}$ and $q_{2}=n_{2}-p_{1}=n_{1}-p_{2}$, and where $C={\rm diag}(c_{1},\ldots,c_{p})$ and $S=\sqrt{I-C^{2}}$ are $p_{1}\times p_{1}$ diagonal, positive definite, and real valued matrices with entries in the range $[0,1]$.  $\Box$

Use of ((72)) in ((70)) provides

[TABLE]

where we exploited the fact that $H_{1}$ commutes with $B^{T}$, and $H_{2}$ commutes with $A^{T}$, as well as the fact that both $A$ and $B$ are orthogonal matrices. By substituting ((73)), we can then explicitly express matrix $H$ in ((74)) in the block–diagonal form

[TABLE]

where

[TABLE]

The eigenvalues in ((52)) are exactly the roots of $H_{0}$, which can be identified by standard techniques since ((76)) has a simple $2\times 2$ block structure.666We skip the derivation since it is trivial, but tedious, and only report the final result. By inspection of ((75)), the remaining eigenvalues are $\overline{n}_{1}\overline{n}_{2}$ and $-\overline{p}_{1}\overline{n}_{2}$, which is consistent with statement 2) given by the theorem.

The extension to the cases where ((71)) is not valid can be performed with a simple substitution of variables. To start with, consider the case where $\min(\overline{p}_{2},\overline{n}_{2})\leq\min(\overline{p}_{1},\overline{n}_{1})$, which entails four different scenarios, namely:

1. a)

$p=p_{2}\leq p_{1}\leq n_{1}\leq n_{2}$, which we already discussed; 2. b)

$p=p_{2}\leq n_{1}\leq p_{1}\leq n_{2}$, 3. c)

$p=n_{2}\leq p_{1}\leq n_{1}\leq p_{2}$, and 4. d)

$p=n_{2}\leq n_{1}\leq p_{1}\leq p_{2}$.

Now, scenario b) can be solved equivalently to scenario a) provided that the matrix $G$ is rearranged in the form

[TABLE]

This corresponds to making the swaps $p_{1}\leftrightarrow n_{1}$ and $\overline{p}_{1}\leftrightarrow\overline{n}_{1}$, and to multiplying by $-1$ the resulting eigenvalues, since the role of positive and negative have been swapped in $H_{1}$. Incidentally, ((77)) implies that the singular values of interest are now those of $G_{2}$. Thanks to the orthogonality of $G$ we have $G_{2}G_{2}^{*}=I-G_{1}G_{1}^{*}$, and therefore the squared singular values of $G_{2}$ are $1-c_{i}^{2}$. This implies a further substitution $c_{i}^{2}\rightarrow 1-c_{i}^{2}$, i.e., the swap $c_{i}\leftrightarrow s_{i}$. Notably, under the swap $\overline{p}_{1}\leftrightarrow\overline{n}_{1}$, the swap $c_{i}\leftrightarrow-s_{i}$, and the multiplication by $-1$, the expression ((52)) remains valid. As a matter of fact, we are simply transforming $H_{0}$ into $H_{0}^{T}$. The other eigenvalues of ((75)), accordingly, become $\overline{n}_{1}\overline{n}_{2}$ and $-\overline{p}_{1}\overline{n}_{2}$, which is consistent with statement 2).

Scenarios c) and d) can be solved equivalently. The remaining cases to consider are those where $\min(\overline{p}_{2},\overline{n}_{2})\geq\min(\overline{p}_{1},\overline{n}_{1})$, in which case we simply need to swap the roles of $i=1$ and $i=2$, and to work with $G^{T}$, but this has no consequences on ((52)), nor on the statement on the remaining eigenvalues.

Appendix B Proof of Theorem 3

In the following, the result ((62)) on the range of complex eigenvalues, and the result ((61)) on the range of real eigenvalues, are separately proved since they require different techniques. Following Appendix A, we denote with $A\sim B$ a similarity between matrices $A$ and $B$ in the sense that they share the same eigenvalues with the same multiplicities, and extensively use the property $AB\sim BA$. The starting point is to identify the similarity (see ((40)))

[TABLE]

where $G=V_{2}V_{1}^{T}$ is as orthogonal matrix, and where $\Lambda_{i}={\rm diag}(\alpha_{i})$ with $\alpha_{i}$ in the ranges ((59)). We also prove the theorem by considering $\underline{n}_{i}>0$ and $\underline{p}_{i}>0$, i.e., that matrices $\Lambda_{i}$ are invertible. The applicability to the case where $\underline{n}_{i}=0$ and/or $\underline{p}_{i}=0$ is guaranteed by taking a limit thanks to the results of perturbation theory (e.g., see [26]).

B-A Complex eigenvalues – Proof of ((62))

Let $x$ be a complex–valued eigenvalue of matrix $G\Lambda_{1}G^{T}\Lambda_{2}$, with unit norm $\|x\|=1$, and with associated eigenvector $\lambda$. This allows writing the general relation

[TABLE]

which is true thanks to the invertibility of $\Lambda_{1}$. In turn, from ((79)) we further obtain

[TABLE]

the first being derived by multiplying ((79)) by $x^{*}$, the second by taking the squared norm. From the second of ((80)) we have

[TABLE]

which is ensured by the invertibility of $\Lambda_{1}$ and $\Lambda_{2}$ (it implies that both $\Lambda_{2}^{2}$ and $\Lambda_{1}^{-2}$ are positive definite, so that both the numerator and the denominator are strictly positive thanks to the assumption $\|x\|=1$). A perfectly equivalent outcome cannot be stated for the first of ((80)), since $x^{*}G\Lambda_{1}^{-1}G^{T}x$ can be equal to zero. In this case we instead have two distinct possibilities, namely: 1)

[TABLE]

which clearly identifies real–valued eigenvalues, and 2)

[TABLE]

which in general identifies complex–valued eigenvalues.

A range for the complex valued eigenvalues can therefore be identified by considering both ((81)) and ((83)), namely the optimisation problem

[TABLE]

which clearly provides and upper bound to $|\lambda|^{2}$. As a simpler alternative we consider relaxing ((84)) by removing the constraint $y=G^{T}x$, which provides two simple optimisation problems of the form

[TABLE]

and

[TABLE]

whose composition still ensures an upper bound to $|\lambda|^{2}$.

We preliminarily concentrate on the solution to ((85)), which we do by denoting $|x_{i}|^{2}=a_{i}$, so that the problem can be equivalently written in the form

[TABLE]

The corresponding Lagrangian assumes the form

[TABLE]

with $\eta_{\geq 0}(x)$ an indicating function for $x\geq 0$. Observe that the Karush Kuhn Tucker (KKT) condition $a^{*}=\max_{a}L(a,\mu^{*})$, where $(a^{*},\mu^{*})$ denotes the primal–dual optimal points, ensures that

[TABLE]

In fact, $g^{*}_{i}<0$ implies $a_{i}^{*}=0$, $g^{*}_{i}=0$ implies any value for $a_{i}^{*}$, while $g^{*}_{i}>0$ implies $a_{i}^{*}=+\infty$ which leads to an infeasible solution. Actually, the KKT condition $\sum_{i}a_{i}^{*}=1$ implies that there must be at least one $g^{*}_{i}=0$, while the KKT condition $\sum_{i}a_{i}^{*}\lambda_{2i}=0$ further implies that the active values must be at least two, one for a positive eigenvalue $\lambda_{2i}$, and one for a negative eigenvalue. Having said so, it is straightforward to see that the quadratic form of $g^{*}_{i}$ as a function of $\lambda_{2i}$ is required to have roots in $\overline{p}_{2}$ and $-\overline{n}_{2}$, in such a way that $g_{i}^{*}=0$ for those $i$ where $\lambda_{2i}=\overline{p}_{2}$ or $-\overline{n}_{2}$, and $g_{i}^{*}<0$ otherwise. As a consequence, it is $a_{i}^{*}=0\,\forall i\hbox{ such that }\lambda_{2i}\neq\overline{p}_{2},-\overline{n}_{2}$. By then using the compact notation

[TABLE]

we can reinterpret ((87)) in the form

[TABLE]

which can be easily solved by exploiting the second constraint, that is by replacing $\alpha=(1-\alpha)\overline{n}_{2}/\overline{p}_{2}$ and $(1-\alpha)=\alpha\overline{p}_{2}/\overline{n}_{2}$ in the target function. The upper bound in ((84)) is therefore simply $\overline{p}_{2}\overline{n}_{2}$.

The solution to ((85)) is perfectly equivalent, with the observation that the counterpart to ((88)) is

[TABLE]

and leads to a choice of $\mu^{*}$ that sets $g^{*}_{i}=0$ for those $i$ associated with eigenvalues $\overline{p}_{1}^{-1}$ and $-\overline{n}_{1}^{-1}$, namely the least positive and greatest negative values, and $g_{i}^{*}>0$ for the rest. This identifies a lower bound $1/(\overline{p}_{1}\overline{n}_{1})$ in ((85)).

By combining the two bounds we therefore have $|\lambda|^{2}\leq\overline{p}_{1}\overline{n}_{1}\overline{p}_{2}\overline{n}_{2}$, which is the upper bound in ((61)). The lower bound can be stated equivalently, by simply investigating the upper bound of the inverse matrix

[TABLE]

where the roles of $\overline{p}_{i}$ and $\overline{n}_{i}$ swap with $\underline{p}_{i}^{-1}$ and $\underline{n}_{i}^{-1}$.

B-B Real eigenvalues – Proof of ((61))

For the real eigenvalues case we exploit the results readily available from Appendix A, and assume that matrices $\Lambda_{i}$ are organised in such a way that the positive values are in the first positions, and are followed by the negative values. This allows writing $\Lambda_{i}$ in the form

[TABLE]

where $H_{i}$ was defined in ((49)) and carries the largest positive and negative values of $\Lambda_{i}$, and where $\breve{\Lambda}_{i}$ is a diagonal positive definite matrix satisfying $0\prec\breve{\Lambda}_{i}\preceq I$. By replacing the equivalence in ((78)), and by exploiting the commutative property between some of the considered matrices, we find that

[TABLE]

where $P_{i}$ are by construction real–valued positive–definite symmetric matrices with $0\prec P_{i}\preceq I$, and where $H$ is the matrix defined in ((75))-((76)).

Now, from Appendix A we know that $H$ is real and diagonalizable, with eigenvalues belonging to the range depicted in Fig. 4 (see details in Section V-A). This ensures that it can be written in the form777$H$ is diagonalisable, hence it can be expressed in the form $UDU^{-1}$ with diagonal $D$ collecting the eigenvalues. The columns of $U$ are the eigenvectors, that can be chosen real–valued for real–valued eigenvalues. For complex valued eigenvalue couples $\{\lambda,\lambda^{*}\}$, instead, the eigenvectors can be chosen to be a complex conjugate couple, i.e., $\{u,u^{*}\}$. To obtain a real–valued form it suffices multiplying the blocks of interest by a normalising unitary matrix $J$, namely

\underbrace{[u,u^{*}]J}_{V}\cdot\underbrace{J^{*}\left[\matrix{\lambda\cr&\lambda^{*}}\right]J}_{B}\cdot J^{*}\;,\quad J=\frac{1}{\sqrt{2}}\left[\matrix{1&i\cr 1&-i}\right]

where matrix $V$ is now real–valued, and where matrix $B$ assumes the real–valued form ((96)).

[TABLE]

where $V_{H}$ is real–valued and invertible, while $\Lambda_{H}$ is a real–valued and normal matrix with a block diagonal structure whose blocks are either:

$1\times 1$ blocks carrying the real–valued eigenvalues of $H$; 2. 2.

$2\times 2$ blocks of the form

[TABLE]

for all the complex–conjugate couples in ((52)).

Hence, we can write

[TABLE]

which reveals the fundamental role of matrix $T$.

Notably, matrix $T$ in ((97)) is a weakly positive matrix, that is, a diagonalisable matrix with non-negative eigenvalues, a result which is a straightforward consequence of the fact that the product of two positive definite matrices is weakly positive [32, 33]. Being weakly positive, $T$ satisfies

[TABLE]

Incidentally, the bounds $P_{1}\preceq I$ and $P_{2}\preceq I$ guarantee that $\|P_{2}P_{1}\|\leq 1$. Since $T\sim P_{2}P_{1}$, the above further implies that the largest eigenvalue of $T$ is bounded by $\lambda_{\max}(T)\leq 1$. Similarly, the invertibility of both $P_{1}$ and $P_{2}$ further guarantee the existence of matrix $T^{-1}$, as well as the fact that its smaller eigenvalue satisfies $\lambda_{\min}(T^{-1})\geq 1$. This implies the property

[TABLE]

which can be simply derived by observing that we can express $T^{-1}$ as $T^{-1}=I+T^{\prime}$ with $T^{\prime}$ some weakly positive matrix.

Given the above, a range for real eigenvalues can be derived by investigating the relation

[TABLE]

where $x$ is a real–valued eigenvector of unit norm, and $\lambda$ is its corresponding real eigenvalue. By multiplying by $x^{T}$, ((100)) allows writing

[TABLE]

where the denominator is ensured by ((99)) to satisfy $x^{T}T^{-1}x\geq 1$, while for the numerator inspection of Fig. 4 ensures $x^{T}\Lambda_{H}x\in[-\overline{n},\overline{p}]$. Overall, ((101)) justifies the extrema of ((61)), namely the relation $\lambda\in[-\overline{n},\overline{p}]$.

The relation $\lambda^{-1}\in[-\underline{n}^{-1},\underline{p}^{-1}]$, which completes the proof, can be derived by an equivalent method that investigates the inverse matrix ((92)).