Giovanni Accongiagioco
email: giovanni.accongiagioco@imtlucca.it

Eitan Altman
email: eitan.altman@inria.fr

Enrico Gregori
email: enrico.gregori@iit.cnr.it

Luciano Lenzini
email: l.lenzini@iet.unipi.it

A Game Theoretical study of Peering vs Transit in the Internet

Keywords: Internet Modeling, Complex Networks, Game Theory, AS-level Internet Topology, Supermodular Games

. Simulation shows that the system, in the vast majority of cases, converges to an equilibrium. Very interestingly, even in asymmetric scenarios the equilibrium reached suggests that players tend to be symmetric with respect to the peering exchange points and send their asymmetric traffic quota via the transit service providers.

I. INTRODUCTION

The Internet consists of thousands of Autonomous Systems (AS), independently administered networks that dynamically connect together to provide end-to-end reachability. Depending on their importance and the offered services, ASes can be sorted in different tiers and categories: content providers, access providers, transit providers and so on [START_REF] Norton | The 2013 Internet Peering Playbook[END_REF]. During the past century, Internet's structure was predominantly hierarchical, with customer ASes paying their providers to carry their traffic, and the latters providing default gateways to reach any requested destination (e.g. an access provider offers network connectivity to the end users, and pays a transit provider to use it as a gateway). This kind of relationship is known as "customer to provider" or transit, and the pricing strategy of providers w.r.t. their customers is typically volume-based [START_REF] Norton | Internet transit prices -historical and projected[END_REF].

Nowadays Internet, as reported in [START_REF] Dhamdhere | The internet is flat: modeling the transition from a transit hierarchy to a peering mesh[END_REF], is evolving from its previous structure and becoming flatter through the introduction of peering meshes, with ASes establishing bilateral agreements to exchange traffic between them for free. This This work has been partially supported by the European Commission within the framework of the CONGAS project FP7-ICT-2011-8-317672.

The first author would like to thank Alexandre Reiffers for the interesting discussions on game theoretic modeling aspects. kind of relationship is known as "peering", and in this case the cost of each AS consists of laying out and maintaining the physical interconnection. In order to avoid exploits, the peering ASes must agree on some policy, which typically forbids the use of each other as default gateway, balances the traffic ratio and so on [START_REF] Norton | A study of 28 peering policies[END_REF]. The massive presence of peering meshes would be impossible if ASes had to lay out physical interconnections one by one, and was achieved thanks to the increasing deployment of Internet Exchange Points (IXPs) [START_REF] Gregori | The impact of IXPs on the AS-level topology structure of the Internet[END_REF]. These public interconnection points can be roughly compared to large peering hubs: when an AS joins an exchange point, it can peer with a subset (or even all) the other ASes which participate to the same IXP. In order to connect to an IXP, customers typically pay a flat, monthly based, fee, which depends on the cost for maintaining the IXP equipment, in proportion to how much they use it: this means that maintainance costs are shared among all the participants of an exchange point [START_REF] Norton | The art of peering -the ix playbook[END_REF]. Please note that the addition of new customers can lower the costs of a single participant, therefore this cost function is quite different from standard modelization, like the one done in [START_REF] Orda | Competitive routing in multiuser communication networks[END_REF].

ASes joining the Internet face the complex question of what is the best strategy for offering their services (e.g. traffic delivery) at the lowest possible cost. While the answer was easy in the past century, due to the existence of transit as the unique interconnection policy, today's answer is much more complex. In fact, peering policies and IXPs brought new variables to the problem, such as the fact that the outcome of an AS decision also depends on what other ASes, dealing with the same problem, do: due to this aspect, we believe that a game-theoretical analysis of the problem would be highly insightful. While, for the sake of mathematical tractability, it is impossible to realistically model the whole decision space of an AS, the model proposed in [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF] keeps into account many realistic factors like port capacity, AS reachability and realistic interconnection costs. The model focuses on access provider ASes who need to set up their connections and route traffic, and highlights the fundamental difference between transit and peering interconnection policies. In fact, a better understanding of this problem is crucial if we really want to grasp the behavior of ASes living in the Internet ecosystem.

The main contributions of this work are twofold. First of all, we extend results presented in [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF], showing the existence of equilibria and an algorithm converging to them even for generalized version of the game therein proposed. Moreover, we use simulation to show that the main findings of [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF], like the emerging competition (first observed in [START_REF] Ager | Anatomy of a large european ixp[END_REF]) between the two facilities enabling either transit or peering connectivity, still hold for more general cases. Even more interestingly, we notice that also for asymmetric scenarios, the equilibrium reached suggests that players tend to be symmetric with respect to the internet exchange points and send their asymmetric traffic quota via the transit service providers.

The remainder of this paper is organized as follows: in Section II we describe the related work, Section III recalls the general model, in Section IV we give convergence results for several scenarios, in Section V we simulate the system's behavior for more general cases, then conclude in Section VI.

II. RELATED WORK

The majority of Internet AS-level topology models are based on graph theory and try to reproduce topological properties observed in the Internet's measured graph, such as its powerlaw degree distribution [START_REF] Faloutsos | On Power-Law Relationships of the Internet Topology[END_REF], the small-world property [START_REF] Watts | Collective dynamics of 'small-world' networks[END_REF] and other structural properties (communities, cliques, etc..) [START_REF] Gregori | k-clique Communities in the Internet AS-level Topology Graph[END_REF]. In order to generate graphs [START_REF] Haddadi | Network Topologies: Inference, Modeling and Generation[END_REF] which exhibit such properties, these kind of models either define attachment criteria, as in [START_REF] Watts | Collective dynamics of 'small-world' networks[END_REF], [START_REF] Yook | Modeling the Internet's large-scale topology[END_REF], [START_REF] Wang | Understanding and Modeling the Internet Topology : Economics and Evolution Perspective[END_REF], [START_REF] Accongiagioco | A Structure-Based Topology Generator for the Internet' s Core[END_REF] or solve constrained optimization problems (see [START_REF] Fabrikant | Heuristically optimized trade-offs: A new paradigm for power laws in the internet[END_REF], [START_REF] Chang | To Peer or not to Peer: Modeling the Evolution of the Internet's AS-level Topology[END_REF]).

The approach of [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF] is quite different, since it tries to understand how connections are formed as the result of a game between ASes. In this context, the realm of network formation games (see [START_REF] Jackson | Social and Economic Networks[END_REF], [START_REF] Goyal | Connections: an introduction to the economics of networks[END_REF]) investigates existence and properties of equilibria in a network created by rational players, each one with their costs and utility functions. Unfortunately, at the present stage, these models cannot be used to study Internet or any other real life network, due to the simplifications needed for mathematical tractability. Taking a different approach, agent-based computational models such as GENESIS [START_REF] Lodhi | GENESIS: An agent-based model of interdomain network formation, traffic flow and economics[END_REF], include more realistic considerations and simulate the dynamic of network formation. These models however, are analytically intractable and it is impossible to give convergence guarantees.

Given all the difficulties arising when trying to model the whole network formation process, [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF] prefers to analyze the decision problem faced by ASes connecting to a pre-existing network, in order to deliver traffic according to some demands. Authors of [START_REF] Jiménez | Competition in access to content[END_REF] present a game-theoretic modelization of a network where access providers optimize their costs when connecting to a subset of content provider. The assumptions made for the cost functions, however, do not take into account the differences between transit and peering interconnection policies, as the work concentrates more on the analysis of network-neutrality issues. An interesting analysis on peering economics and private internet exchanges is made in [START_REF] Shakkottai | Economics of network pricing with multiple isps[END_REF]. Unfortunately the "private peering" therein analyzed is quite different, both in rules and costs, from the "public peering" brought by IXPs, and is unable to explain nowadays phenomenons like the peering meshes (see the Introduction). In order to investigate the present-state of peering the model in [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF] explicitely introduces the concept of exchange point.

ISP 1 NSP 1 NSP l 1 CP 1 ISP I ISP i IXP L IXP l CP n CP N φ n i φ 1 i φ N i c' n t l (x i,l ,x l ) x n i,l c i,l b n l t 1 (x i,1 ,x 1 ) Figure 1. General Model III. GENERAL SCENARIO
In the following we briefly describe the general model proposed in [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF]. The reader can refer to section III of that report for a comprehensive dissertation.

A. Description

First of all, ASes have been categorized according to what follows (as also done by PeeringDB [24]):

Internet Service Provider (ISP)1 : access provider giving to the end users a way to reach the Internet and its contents.

Content Provider (CP) 1 : which has physical access to the contents users are looking after.

Internet Exchange Point (IXP) 1 : a facility that provides peering connection to all its participants, in order to reach every CP connected to the IXP.

Network Service Provider (NSP) 1 : a facility providing transit connection to its customers, in order to reach every possible CP.

The general network has i 2 {1,...,I} ISPs, n 2 {1,...,N} CPs and l 2 {1,...,L} transmit facilities (TF), which can be either NSPs or IXPs. Without loss of generality, it is supposed that TFs j 1 2 {1,...,l 1 } are NSPs, while TFs j 2 2 {l 1 + 1,...,L} are IXPs. Each ISP i has a demand for a CP n, indicated as φ n i . The players of the game are the ISPs, which need to decide how to split their demand among all possible transmit facilities, as shown by figure 1.

Variable x n i,l indicates the flow from ISP i to CP n through TF l, therefore the strategy of ISP i is given by the vector

x i =(x 1 i,1 ,...,x n i,l ,...,x N i,L ) 2 R L⇥N
, while the strategy of all the other players is expressed as x \i 2 R (I-1)⇥L⇥N . The goal of each player is to serve, at the minimum possible cost, his demand

φ i =( φ 1 i ,...,φ N i ) 2 R N by splitting it into several flows x n i,l .
The cost of player i is given by the summation of three costs:

TF usage cost : which depends on the transmit facility (NSP or IXP) used. As shown in [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF], the NSP usage cost is linear in the amount of flow that each player sends to it:

a l P n x n i,l l  l 1 (1)
where a l, ll1 represents the transit price of NSP l per unit of flow. The IXP participation cost, as previously said and furtherly clarified in [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF], depends on the cost for maintaining the IXP equipment, in proportion to how much each player uses it, therefore it can be written as:

P n x n i,l P i P n x n i,l h l (x l ) l>l 1 (2)
where h l represents the maintenance cost of IXP l as a function of the total flow through the IXP. As detailedly shown in [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF], the mantainance cost can be approximated with that of ports, which is shown to be fitted by a square root:

h l (x l )=a l q P i P n x n i,l
, where a l, l>l1 is a constant relating the total flow through IXP l with its mainteinance cost.

TF capacity cost : which arises from the fact that flow on each link between an ISP and a TF cannot exceed a fixed capacity c i,l . With the proper assumptions, illustrated in [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF], the constraint is modeled as a cost growing to infinity as the flow over the link approaches the capacity limit:

1 c i,l - P n x n i,l P n x n i,l (3)
Typically the capacity of an NSP can be assumed to be much larger than that of an IXP: c NSP c IXP (see [START_REF] Norton | Internet transit prices -historical and projected[END_REF], [START_REF] Norton | The art of peering -the ix playbook[END_REF]), moreover we assume that capacities are symmetric w.r.t. the players, therefore c i,l = c l 8i.

CP reachability cost : as previously said, while NSPs can be used to reach every possible CP, through IXPs a player can only reach the subset of CPs connected to the same IXP they participate into. If we indicate with b n l the cost of transporting one unit of flow from TF l to CP n, we can express this constraint by putting:

b n l = ( 0 if l  l 1 _ IXP l ! CP n 1 otherwise (4)
After some algebraic manipulations, explicitely shown in [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF], the cost function of each player can be written as:

C i (x i , x -i )= P l P n x n i,l h f l ⇣ P i P n x n i,l ⌘ + +g l ⇣ P n x n i,l ⌘ + b n l i (5) where: 8 > > < > > : f l ⇣ P i P n x n i,l ⌘ = ( a 1 l  l 1 a l pP i P n x n i,l l>l 1 l ⇣ P n x n i,l ⌘ = 1 c l - P n x n i,l (6)
Equation ( 5) gives the general expression of the cost function for each player and, as already observed in [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF], it is a nonconvex function of x i,l , therefore we cannot directly establish existence of pure equilibria using general frameworks like the one described in [START_REF] Orda | Competitive routing in multiuser communication networks[END_REF].

The best response of player i, BR i (x \i ), is obtained by both minimizing this cost function and satisfying the flow constraints, that is, the total flow served for every CP has to be equal to the demand φ n i . Therefore we have:

( BR i (x \i ) = arg min x i C i (x i , x \i ) s.t. P l x n i,l = φ n i 8n (7)
The vector x ⇤ =( x ⇤ 1 ,...,x ⇤ I ) 2 R I⇥L⇥N is an equilibrium of the game if and only if x ⇤ i 2 BR i (x ⇤ \i ) 8i, that is, if the strategy of each player is a best response to the strategies of other players.

IV. ANALYTIC RESULTS

A. Base Theoretical Results

Theoretical results for the model, regarding the existence of equilibria and algorithms that converge to them, are based on the theory of supermodular games. Here we briefly review the defition of supermodular games and the main results of [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF]. Definition 1. Supermodular games [START_REF] Yao | S-modular games, with queueing applications[END_REF] Consider a generic game G, where players' payoffs are given by an utility function u : R k ! R. When u is twice continuously differentiable, the game is said supermodular if the mixed second derivatives of the utility function are always positive. In the proposed model players are minimizing costs rather than maximizing utilities, therefore this kind of game is supermodular iff:

∂ 2 C(x) ∂x i ∂x j  0 8i 6 = j (8)
Definition 2. Symmetric supermodular games [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF] While supermodular games require property (8) to hold for every possible strategy x i , symmetric supermodular games, as defined in [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF], are those for which [START_REF] Orda | Competitive routing in multiuser communication networks[END_REF] only holds along the symmetric axis, that is when x i = x j 8i 6 = j.

Theorem 1 of [START_REF] Altman | S-modular games and power control in wireless networks[END_REF] proves the existence of equilibria for supermodular games. It also provides a way of computing them, as it shows that due to the monotonicity of best response sequences, they converge to a limit which is shown to be an equilibrium. Along the same line of reasoning, Theorem 4 of [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF] relaxes this result, showing that it holds also for symmetric supermodular games, where equilibria are given as the limit of symmetric best response sequences (defined in Section IV.A of [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF]).

Starting from the general model presented in section III, [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF] defines the Minimal Complexity Model (MCM), which has two ISPs, two CPs, two TFs (an NSP and an IXP), and a fully connected topology. It is shown that the MCM is a symmetric supermodular game (the proof is based on showing that (8) holds along the symmetric axis), and therefore it has at least one pure equilibrium for symmetric demands, given as the limit of a symmetric best response sequence.

While interesting on its own, as it provides a way to highlight basic differences between transit and peering policies (and shed light on the competition between NSPs and IXPs, first observed in [START_REF] Ager | Anatomy of a large european ixp[END_REF]), one might argue that the MCM topology is quite small and simplified to represent the Internet. Here we explicitely tackle this problem by both extending the convergence results to more general scenarios, and using simulation to cover other cases for which we couldn't obtain analytic results, thus analyzing a broad variety of general cases which weren't covered in [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF].

B. Extended Results

Some of the scenarios for which the theorems shown in the previous section did not apply where tested through simulation in [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF], and the observed results brought to the conjecture that convergence was indeed a more general property. As a matter of fact, we are able to generalize Theorem 5 of [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF], showing that symmetric supermodularity holds also for the general scenario of section III.

We consider a system with I ISPs, N CPs, and L TFs, with l 1 =1 , in a possibly disconnected topology. While having multiple IXPs is fundamental for understanding how players aggregate around exchange points, especially in presence of reachability constraints (IXPs may be connected to only a subset of CPs), this is not the case for NSPs, due to the fact that their cost is independent from other players' choice, and they are necessarily connected to all possible CPs. Therefore, without loss of generality for our problem, it is safe to consider only a single NSP for l =1 , and L -1, for l =2 ..L, IXPs as we do here. In the following we will always refer to this system, which is a single NSP version of figure 1.

The cost function of section III can thus be rewritten separating the NSP component from the IXPs, therefore:

C i (x i , x -i )= P n x n i,1 ⇣ a 1 + 1 c1- P n x n i,1 ⌘ + (9)
+ P l6 =1 P n x n i,l ✓ a l pP i P n x n i,l + 1 c l - P n x n i,l + b n l ◆
where the CP reachability cost for the NSP has been removed since we know from (4) that b n 1 =08n. Theorem 3 and Corollary 4, whose complete proofs can be found in the full version [START_REF] Accongiagioco | A game theoretical study of peering vs transit in the internet[END_REF], demonstrate existence of equilibria and convergence of the symmetric best response algorithm for the general case just formalized. Theorem 3. The game defined in ( 9) is symmetric supermodular.

Sketch of the Proof: The proof is based on showing that (8) holds along the symmetric axis for any possible combination of indexes:

∂ 2 C i (x i , x -i ) ∂x n j, l∂x n i,l x n i,l =x n j, l  0 8i 6 = j, 8l, l, n, n (10)
Corollary 4. The game defined in ( 9) has at least one pure equilibrium for symmetric demands, given as the limit of a symmetric best response sequence.

C. Subcases Analysis

After proving existence of equilibria for the general case, here we analyze some subcases in order to understand what kind of equilibria we should expect for specific scenarios. We have two main categories: a) Fully Connected Topologies: Suppose that we have a fully connected topology, meaning that b n l =08l, n. In such a case, we can take the summation over n and consider cumulative flows and demands:

( x i,l = P n x n i,l cum. flow ISP i ! TF l φ i = P n φ n i cum. demand ISP i
We can now substitute these two variables inside cost function [START_REF] Norton | A study of 28 peering policies[END_REF], thus obtaining an equivalent problem where the strategy of each player is a vector x i =( x i,1 ,...,x i,l ,...,x i,L ) 2

a 1 =2, a 2 =2, φ =2, c 1 =10, c 2 =3
R L . This means that, in fully connected topologies, our system is equivalent to another one where we only have a single CP, and each player has to serve a cumulative demand φ i for that CP. This happens because there are no reachability constraints, therefore from a player's perspective the specific CP from which he has to fetch data is not relevant.

b) Symmetric IXPs : Suppose that all the IXPs have the same costs, capacities and reachability matrix:

a l = a IXP , c l = c IXP , b n l = b n IXP 8l 6 =1.
Due to their symmetry, there is an equilibrium where traffic is split equally among them [START_REF] Cheng | Notes on equilibria in symmetric games[END_REF], therefore we might think of transforming this problem in an equivalent one having a single IXP with the same reachability matrix and transformed costs and capacities. Unfortunately, we were unable to perform this conversion due to the form of the cost function for the IXPs. In fact, as we see from ( 5) and ( 6), the non linear port cost h l makes it quite different for players to share small traffic quantities rather than large ones.

The analysis of the two categories highlights that scenarios with multiple CPs can be highly simplified with fully connected topologies, while in the case of multiple IXPs, even if symmetric, the analysis can be quite difficult. To conclude this section, we analyze in more detail the simplified case where we have fully connected topology and just one IXP. Due to the fact that we can handle the multiple CPs as if there was just a single one, this scenario is quite similar to the MCM, except that we have a generic number I of players. However, thanks to Corollary 4, we now know that for symmetric demands φ i = φ 8i we can compute equilibria as the limit of symmetric best response sequences. Just like it was done for the MCM, we can use algorithm 1 of [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF] to compute the equilibria.

With the purpose of analyzing number and position of equilibria, for the case of I =3(it would be difficult to represent more dimensions), we can draw the Best Response Intersection (BRI) graph, shown in Figure 2. Just like in the MCM, the picture shows three equilibrium points, obtained by the intersection of the three surfaces representing the players' best responses. As we see from the straight line crossing all such points, the three equilibria are symmetric, with the leftmost (traffic split between NSP and IXP) and the rightmost (all flows through the NSP) being the stable ones.

V. S IMULATIVE RESULTS

In this section we show results of simulations, performed both for symmetric scenarios, for which convergence has been proven in the general case, and asymmetric scenarios, for which we have no proofs. In fact, as we'll see later on, in this last case it is possible for players to never reach an equilibrium.

Simulations have been performed using the same symmetric best response algorithm shown in [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF]: iteratively, each player performs its best response to the set of other players' strategies. If the simulation converges, the output is the equilibrium for the given input parameters, which are:

• the number of ISPs, TFs and CPs, respectively I, L, N ;

• the cost function parameters a l , c l , b n l and demands φ n i . Moreover, it is important to set a startingpoint, because, as shown in [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF], on startup IXPs need a critical mass, represented by a share of the total traffic in the system, in order to be able to attract players.

A. Symmetric Case

We simulate a scenario with I = 10 ISPs, L =4TFs (L-1 symmetric IXPs) and N =4CPs. The connectivity matrix is:

b n l = ( 1 (l, n)=(2, 2) _ (l, n)=(4, 1) _ (l, n)=(4, 4) 0 otherwise
The cost parameter a l has been chosen in order to be similar to present reality. We use [START_REF] Norton | Internet transit prices -historical and projected[END_REF] for inferring a general transit cost of a 1 = a NSP =3, and [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF] for the IXP cost 8l 6 =1a l = a IXP = 30, using the fit therein obtained for the port costs of the MIX, a specific italian IXP. All users have symmetric flows φ n i = 12.5 8i, n and their capacities to the facilities are c 1 = c NSP = 100 and 8l 6 =1c l = c IXP = 20, so that c NSP c IXP . As already happened in the MCM, depending on the startingpoint we notice the existence of multiple equilibria. In fact, if the initial condition is such that one or more IXPs are underutilized, than at equilibrium those IXPs will not be used. This phenomenon corroborates the outcomes of [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF], showing that indeed even in general scenarios the competition between NSPs and IXPs, and even between IXP themselves, strongly emerges. Differently from the MCM, in this case we observe more than two stable equilibria, since any combination with one or many unused IXPs can be an equilibrium.

Suppose now that the startingpoint is such that flows are split equally among the facilities, so that all IXPs have the critical mass to attract players. Figure 3 shows the scatterplot at equilibrium. In this plot, each dot represents the flow quantity that each user sends on a given path (that is, to a fixed CP through a given IXP). Due to symmetry, we observe that all users will behave symmetrically on the same path, and this is exactly what happens in Figure 3. There is generally a low utilization of the NSP, which rises a little bit for those CPs with a worse reachability matrix (CP 1 ,C P 4 ).

B. Asymmetric Case

We now show the impact of asymmetric players' demands. In this case, convergence of the best response sequence is not guaranteed by Corollary 4, however, we know that if the simulation converges we certainly reach an equilibrium [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF].

We simulate a scenario with exactly the same parameters as in the symmetric case, except that now the demands grow linearly from φ n 1 = 10 8n to φ n I = 15 8n. The average demand is still 12.5, but now the demand of the last player is 1.5 times that of the first one. The scatterplot at equilibrium is shown in Figure 4. Very interestingly, even if demands are asymmetric, paths of flow tend to be almost symmetric for the IXPs, while they spread apart for the NSP. This happens because the benefits of sharing costs at the exchange points is bigger when the traffic ratio is approximately the same between participants, therefore players tend to "symmetrize" around the IXPs. Due to the fact that flows around the IXPs are more or less symmetric, players will send the traffic residual through the NSP, which will see highly asymmetric patterns.

In previous case, the asymmetry in players' demands was not very pronounced. Let's now see what happens when the demands go from φ n 1 =6 .5 8n to φ n I = 18.5 8n, meaning that last player demand is nearly three times that of player one. Again, Figure 5 shows the scatterplot at equilibrium. Due to the heavily unbalanced demands, the symmetric patterns around the IXPs are still present, but much less pronounced. While in the previous case equilibrium was driven by the simple rule of "symmetric behavior", in this case the outcome is more difficult to predict. In general, due to the asymmetry, cost benefits of players' for using exchange points decrease, therefore we observe, on average, an higher quantity of flow going through the NSP.

The phenomenons emerged through this analysis can provide some preliminary insight on how to devise optimal policies to handle peering traffic at IXPs. More specifically, the "symmetric behavior" rule highlights that it is beneficial to balance traffic as much as possible, therefore IXP owners should create few classes of traffic (namely, few different port sizes), and participants should try to aggregate traffic on these ports, since unbalanced flows must be handled at NSPs and bring to suboptimality. While simulations have been carried out with a limited network size due to computational constraints, the conclusions are fairly general, therefore we expect similar results to hold for larger-scale scenarios.

C. Non Convergence

Even quite simple scenarios for which we cannot apply Theorem 3, might lead to a situation where players' behavior oscillates, never reaching an equilibrium. One of such cases can be found in the full version of the paper [START_REF] Accongiagioco | A game theoretical study of peering vs transit in the internet[END_REF]. Consider a system with two symmetric IXPs and an asymmetric starting point s.t. a group of players send more traffic to one of the IXPs and less to the other, while the other group of players do the opposite. Due to the asymmetric assumption we cannot apply Corollary 4, and simulation shows that this scenario might never reach an equilibrium. This happens when players enter a never-ending oscillation between the first and the second group, as detailedly shown in [START_REF] Accongiagioco | A game theoretical study of peering vs transit in the internet[END_REF].

VI. CONCLUSIONS AND FUTURE WORK

This work, starting from [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF], extends the analysis therein proposed for the MCM model to more general scenarios which better fit the complex nature of the Internet. From a game theoretic perspective, we were able to prove that even the general model falls in the category of Symmetric Supermodular games, thus obtaining both existence of equilibria and convergence of symmetric best response sequences for our game. Moreover, we have shown that for asymmetric sequences convergence may fail, by explicitely simulating a specific scenario for which the system enters in a never-ending oscillation. From an engineering perspective, the contributions are twofold. First of all, we corroborate the main result of [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF] by highlighting once more the competition between IXPs, providing customers the ability to lay out peering connections, and NSPs, high level providers selling transit connections. Second, we have shown that also for asymmetric cases the system often reaches an equilibrium. Even more interestingly, the equilibrium reached in asymmetric cases suggests that players tend to "symmetrize" as much as possible with respect to the peering exchange points and send their asymmetric traffic quota via the transit service providers. This observation can provide insights on how to devise optimal policies to handle peering traffic at IXPs.

The work calls for extensions in at least two directions, actively pursued at the time of writing. The first consists of an in-depth study of the properties of asymmetric cases, e.g. the analysis of their optimality using metrics like the price of anarchy, stability and fairness, already considered for the MCM. The second, taken directly from [START_REF] Accongiagioco | Peering vs transit: a game theoretical model for autonomous systems connectivity[END_REF], follows from the dynamic nature of the Internet. In reality, both transit and peering costs can be renegotiated over time, therefore by adding as players the TFs we can analyze a dynamic game where their strategies are intertwined with that of ISPs.

Figure 2 .

 2 Figure 2. BRI for 3 Players:a 1 =2, a 2 =2, φ =2, c 1 =10, c 2 =3

ISPs, NSPs and CPs are typically ASes. IXPs are not ASes, even if their infrastructure is under a single administrative control.


F Lauer

H Ohlsson

This paper deals with sparse phase retrieval, i.e., the problem of estimating a vector from quadratic measurements under the assumption that few components are nonzero. In particular, we consider the problem of finding the sparsest vector consistent with the measurements and reformulate it as a group-sparse optimization problem with linear constraints. Then, we analyze the convex relaxation of the latter based on the minimization of a block ℓ1-norm and show various exact recovery and stability results in the real and complex cases. Invariance to circular shifts and reflections are also discussed for real vectors measured via complex matrices.

Sparse phase retrieval via group-sparse optimization 1 Introduction

The problem of recovering a signal from quadratic measurements is known as phase retrieval. A typical case with many applications, for instance in optics [START_REF] Walther | The question of phase retrieval in optics[END_REF] or crystallography [START_REF] Harrison | Phase problem in crystallography[END_REF], is when the measurements correspond to the (squared) magnitude of the Fourier transform of the signal. Here, we consider the more general setting of arbitrary quadratic measurements, y i = x H Q i x, i = 1, . . . , N , while focusing on the case where the signal x ∈ C n is assumed to be sparse, i.e., with few nonzero entries in x. As in compressive sensing [START_REF] Donoho | Compressed sensing[END_REF][START_REF] Candès | Compressive sampling[END_REF], which deals with the recovery of sparse signals from linear measurements, the sparsity prior reduces the number of measurements required to recover the signal.

Related work. Seminal works on phase retrieval [START_REF] Kohler | Source reconstruction from the modulus of the correlation function: a practical approach to the phase problem of optical coherence theory[END_REF][START_REF] Gerchberg | A practical algorithm for the determination of phase from image and diffraction plane pictures[END_REF][START_REF] Gonsalves | Phase retrieval from modulus data[END_REF][START_REF] Fienup | Phase retrieval algorithms: a comparison[END_REF] did not consider the sparsity prior. Though these methods were able to incorporate prior information on the support, they were not designed to estimate the support with limited information on its size. More recently, matrix lifting techniques were developed in [START_REF] Candès | Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming[END_REF][START_REF] Candès | Phase retrieval via matrix completion[END_REF][START_REF] Waldspurger | Phase recovery, maxcut and complex semidefinite programming[END_REF] for phase retrieval and more particularly for the case of sparse signals in [START_REF] Shechtman | Sparsity based sub-wavelength imaging with partially incoherent light via quadratic compressed sensing[END_REF][START_REF] Jaganathan | Recovery of sparse 1-d signals from the magnitudes of their Fourier transform[END_REF][START_REF] Ohlsson | Quadratic basis pursuit[END_REF]. The basic idea is to apply a change of variable resulting in linearized measurement equations with a rank-1 constraint on the new matrix variable X = xx H . Then, the rank-1 constraint is relaxed to the problem of minimizing the rank of the matrix, which is further relaxed to the minimization of the nuclear norm. In these methods, the sparsity prior is typically incorporated as the minimization of an ℓ 1 -norm, which induces a trade-off between the satisfaction of the rank-1 constraint and the sparsity of the solution. Other methods focusing on sparsity are typically iterative, like Fienup-type methods [START_REF] Mukherjee | An iterative algorithm for phase retrieval with sparsity constraints: application to frequency domain optical coherence tomography[END_REF] using alternate projections or the GESPAR method [START_REF] Shechtman | GESPAR: Efficient phase retrieval of sparse signals[END_REF] implementing a local search strategy with a bi-directional greedy algorithm. These iterative methods usually come without recovery guarantees.

Finally, note that sparse phase retrieval also enters the more general framework of nonlinear compressed sensing, as investigated in [START_REF] Ohlsson | Nonlinear basis pursuit[END_REF] for analytic functions computing the measurements, in [START_REF] Ehler | Quasi-linear compressed sensing[END_REF] for quasi-linear functions and in [START_REF] Lauer | Finding sparse solutions of systems of polynomial equations via group-sparsity optimization[END_REF] for polynomials.

Contribution.

We propose a convex approach to the sparse phase retrieval problem. This approach relies on two main steps: the linearization of the constraints inducing a group-sparse structure on the variables and a convex relaxation of the group-sparse optimization problem enforcing this structure. More precisely, the linearization is based on the Veronese map lifting the signal to a higher dimensional space. This map is invariant to a global sign change and preserves sparsity in the sense that the lifted signal has a group-sparsity similar to the sparsity of the original signal. Thus, the proposed method amounts to estimating a group-sparse signal satisfying linear constraints, from which the original signal can be recovered. This estimation relies on a convex relaxation of the group-sparse problem based on a sum of norms, or block ℓ 1 -norm. Thus, while methods based on matrix lifting lead to semi-definite programming problems, the proposed approach yields a more amenable second-order cone programming formulation. This formulation is also easily extended to deal with noisy measurements. In addition to these algorithmic benefits, we derive exact recovery conditions in the noiseless case and stable recovery guarantees in the presence of noise.

Note that the approach taken here is similar in spirit to the one derived in [START_REF] Lauer | Finding sparse solutions of systems of polynomial equations via group-sparsity optimization[END_REF] for the more general problem of finding sparse solutions of polynomial systems of equations. However, the analysis in [START_REF] Lauer | Finding sparse solutions of systems of polynomial equations via group-sparsity optimization[END_REF] is limited to the real case and does not apply to polynomials without linear terms as the ones found in phase retrieval.

Paper organization. For the sake of clarity, we first detail in Sect. 2 the proposed method in the real case before extending it in Sect. 3 to the complex case. The effect of noise and stability results are discussed in Sect. 4. The case of real signals measured via complex vectors is considered in Sect. 5 which also deals with the invariance of the measurements to circular shifts and reflections. Finally, Section 6 tests the proposed methods in numerical experiments.

Notations. Matrices are written with bold uppercase letters and vectors in bold lowercase letters, except for the ith column A i of a matrix A. The notation (A) i,j denotes the element at the ith row and jth column of a matrix A. ℜ(•) and ℑ(•) denote the real and imaginary parts of a complex number, vector or matrix, and i the imaginary unit. The superscripts T and H denote the transpose and conjugate transpose, respectively, i.e., z H = z T . • p denotes the ℓ p -norm in R n , while • denotes the norm in C n induced by the inner product as z = √ z H z. The ℓ 0 -pseudo-norm of a real or complex vector x of dimension n is defined as x 0 = |{j ∈ {1, . . . , n} : x j = 0}| and denotes the number of nonzero components x j . We also define the ℓ 0 -pseudo-norm of a vector-valued sequence

{u i } N i=1 as {u i } N i=1 0 = |{i ∈ {1, . . . , N } : u i = 0}|.

The real case

We write the sparse phase retrieval problem as

min x∈R n x 0 (1)
s.t. y i = (q T i x) 2 = x T Q i x, i = 1, . . . , N,
where y i ∈ R are the measurements and Q i = q i q T i ∈ R n×n . Due to symmetry, solutions to (1) are defined up to their sign and the goal is to obtain an estimate x = ±x 0 , for x 0 in the solution set of (1). In particular, we are interested in the case where (1) has a unique pair of solutions {x 0 , -x 0 }, while conditions ensuring such a uniqueness are discussed in [START_REF] Ohlsson | On conditions for uniqueness in sparse phase retrieval[END_REF][START_REF] Ranieri | Phase retrieval for sparse signals: Uniqueness conditions[END_REF].

The proposed method relies on two subsequent relaxations. While the first one linearizes the constraints, the second one convexifies the objective function.

First level of relaxation

Let the Veronese map of degree 2, ν : R n → R M , be defined by

ν(x) = [x 2 1 , x 1 x 2 , . . . , x 2 2 , x 2 x 3 , . . . , x 2 n-1 , x n-1 x n , x 2 n ] T ,
and the subscript ij denote the index of its component equal to

x i x j , i.e., ij = min {i,j}-1 k=1 (n -k + 1) + |j -i| + 1 = min {i,j}-1 k=1 (n -k) + min {i, j} + |j -i|. ( 2
)
This notation is symmetric, i.e., ij and ji denote the same index, and will be used throughtout the paper to index the components of vectors of R M or C M . The constraints of the phase retrieval problem (1) can be rewritten as

Aν(x) = y with A ∈ R N ×M and M = n + 1 2 = n(n + 1)/2.
Let W j be an n × M -binary matrix such that1 W j ν(x) = x j x, i.e., W j ν(x) is the vector of n entries corresponding to the monomials in ν(x) including x j . Then, we have

x j = 0 ⇔ W j ν(x) = 0
and the objective function in (1) can be written as

x 0 = {W j ν(x)} n j=1 0 . (3)
Thus, (1) can be reformulated as the nonlinear group-sparse optimization problem

min x∈R n {W j ν(x)} n j=1 0 (4)
s.t. Aν(x) = y.
Note that ν(x) = ν(-x), but that for x = ±x 0 , ν(x) = ν(x 0 ). Thus, the problem can be posed as the one of recovering the value of ν(x 0 ), from which x 0 can be inferred up to its sign.

To estimate ν(x 0 ), we relax (4) to

min v∈R M {W j v} n j=1 0 (5) s.t. Av = y v jj ≥ 0, j = 1, . . . , n,
where the variables in v estimating the components of ν(x) are not constrained to be interdepent monomials of n base variables, but the last constraints in (5) nonetheless ensure that the v jj 's estimating the x 2 j 's are positive.

Convex relaxation

Problem ( 5) is a (linear) group-sparse optimization problem with highly overlapping groups. While groupwise-greedy algorithms, such as the one proposed in [START_REF] Lauer | Finding sparse solutions of systems of polynomial equations via group-sparsity optimization[END_REF], can be applied, their analysis is not available for the case of overlapping groups. Therefore, here, we consider the convex relaxation approach which aims at solving (5) via the following surrogate formulation:

v = arg min v∈R M n j=1 W j W v 2 (6) s.t. Av = y v jj ≥ 0, j = 1, . . . , n,
where we introduced the diagonal matrix W of precompensating weights (W ) i,i = w i = A i 2 , and which can be solved efficiently by off-the-shelf Second-Order Cone Programming (SOCP) solvers.

Then, we easily obtain an estimate of x 0 from the estimate v of ν(x) as x = ν -1 (v), where the inverse mapping ν -1 is defined as

ν -1 (v) =      1 √ v ii [v 1i , v 2i , . . . , v ni ] T , if i > 0 and v 2 ji v ii = v jj , ∀j ∈ {1, . . . , n} 0, otherwise, where i =    min j∈{1,...,n} j, s.t. v jj > 0, if ∃j such that v jj > 0 0, otherwise.
In the above, the first nonzero entry (of index i) is assumed to be positive to fix the signs and make ν -1 injective. This definition also ensures that

ν -1 (ν(x)) = ±x and that ν -1 (v) 0 ≤ {W j v} n j=1 0 (7) since, for ν -1 (v) = 0, W j v = 0 ⇒ v ji = 0 ⇒ ν -1 (v) j = 0.

Analysis

We now turn to theoretical guarantees offered by the proposed approach. First, the following theorem provides the rationale for tackling the sparse phase retrieval problem via the group-sparse optimization formulation [START_REF] Donoho | Compressed sensing[END_REF].

Theorem 1. If the solution v * to (5) is unique and yields x * = ν -1 (v * ) = 0 such that y i = (x * ) T Q i x * , i = 1, . . . , N , then {x * , -x * } is the unique pair of solutions of (1).

Proof. Assume there is an x 0 = ±x * satisfying the constraints of (1) and at least as sparse as x * . Then, Aν(x 0 ) = y and, by using (3) and [START_REF] Ehler | Quasi-linear compressed sensing[END_REF],

{W j ν(x 0 )} n j=1 0 = x 0 0 ≤ x * 0 ≤ {W j v * } n j=1 0 ,
which contradicts the fact that v * is the unique solution to (4) unless ν(x 0 ) = v * . But since x 0 = ±x * , we have x 2 0j = (x * j ) 2 for some j ∈ {1, . . . , n}, which implies (ν(x 0 )) jj = (ν(x * )) jj = ν(ν -1 (v * )) jj . Therefore, by using Lemma 1 in Appendix A with the assumption x * = ν -1 (v * ) = 0, there cannot be such an x 0 .

The following results regarding the convex formulation [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF] are based on the notion of mutual coherence.

Definition 1. The mutual coherence of a matrix

A = [A 1 , . . . , A M ] ∈ R N ×M is µ(A) = max 1≤i<j≤M |A T i A j | A i 2 A j 2 .
With this definition, we can state an exact recovery result (proof given in Appendix B.1).

Theorem 2. Let x 0 be such that y i = x T 0 Q i x 0 , i = 1, . . . , N , and v 0 = ν(x 0 ). If the condition

x 0 0 < 1 2 √ n 1 + 1 µ 2 (A)
holds, then v 0 is the unique solution to [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF].

Corollary 1. Let x 0 be a feasible point of [START_REF] Candès | Compressive sampling[END_REF]. If the condition

x 0 0 < 1 2 √ n 1 + 1 µ 2 (A)
holds, then {x 0 , -x 0 } is the unique pair of solutions to the minimization problem (1) and they can be computed as x 0 = ±ν -1 (v) with v the solution to [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF].

Proof. Assume there exists another solution x 1 = ±x 0 to (1), and thus with x 1 0 ≤ x 0 0 . Then, Theorem 2 implies that both ν(x 1 ) and ν(x 0 ) are unique solutions to [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF] and thus that ν(x 1 ) = ν(x 0 ) = v. But this contradicts the definition of the mapping ν implying ν(x 1 ) = ν(x 0 ) whenever x 1 = ±x 0 . Therefore the assumption x 1 = ±x 0 cannot hold and {x 0 , -x 0 } is the unique pair of solutions to (1), while ν -1 (v) = ν -1 (ν(x 0 )) = ±x 0 .

The complex case

Consider now the problem in complex domain:

min x∈C n x 0 (8) s.t. y i = |q H i x| 2 , i = 1, . . . , N, where y i ∈ R and q i ∈ C n .
The equations in the problem above are invariant to multiplication by a unit complex scalar z with |z| = 1. Thus, there are sets of solutions of the form T (x 0 ) = {x ∈ C n : x = zx 0 , z ∈ C, |z| = 1}, and the goal is to obtain an estimate x ∈ T (x 0 ), from which T (x 0 ) = T (x) can be inferred due to the property of the invariance set:

∀x ∈ T (x 0 ), T (x) = T (x 0 ), (9)
which can be proved as follows. Let x = zx 0 with |z| = 1, then

T (x) = {a ∈ C n : a = bx, b ∈ C, |b| = 1} = {a ∈ C n : a = bzx 0 , b ∈ C, |b| = 1} = {a ∈ C n : a = cx 0 , c ∈ C, |c| = 1} = T (x 0 ).

First level of relaxation

As for the real case, the linearization of the equations will use the Veronese map, which we redefine for complex vectors as follows.

Definition 2 (Complex Veronese map). The complex Veronese map ν :

C n → C M is defined by ν(x) = [x 1 x 1 , x 1 x 2 , . . . , x 2 x 2 , x 2 x 3 , . . . , x n-1 x n-1 , x n-1 x n , x n x n ] T ,
for which the subscript ij , defined as in [START_REF] Candès | Phase retrieval via matrix completion[END_REF], denotes the component index such that (ν(x)) ij equals either x i x j or x j x i (note that, for all pairs (i, j), there is exactly one such component).

Note that ν(x) = ν(x ′ ) for all x ′ ∈ T (x), since for 1 ≤ i ≤ j ≤ n, x ′ i x ′ j = zx j zx j = |z| 2 x i x j = x i x j , but that x ′ / ∈ T (x) implies ν(x) = ν(x ′ ), since x ′ / ∈ T (x) ⇒ |x ′ j | = |x j | ⇒ x ′ j x ′ j = x j x j .
Then, we define the inverse mapping as follows.

Definition 3 (Inverse complex Veronese map). The inverse complex Veronese map, ν -1 :

C M → C n , is defined by ν -1 (v) =    1 √ v ii [v 1i , v 2i , . . . , v ni ] T , if i > 0 and |v ji | 2 v ii = v jj , ∀j ∈ {1, . . . , n} 0, otherwise, where i =    min j∈{1,...,n} j, s.t. ℜ(v jj ) > 0, ℑ(v jj ) = 0, if ∃j such that ℜ(v jj ) > 0, ℑ(v jj ) = 0 0, otherwise.
In particular, we have ν -1 (ν(x)) ∈ T (x).

Note that the square root acts on a real and positive number v ii and is thus also a real positive number. This implies x i = √ v ii ∈ R + and that we arbitrarily set ℑ(x i ) = 0 to fix the value of z in the equation ν -1 (ν(x)) = zx, for a complex number z with |z| = 1, and thus make ν -1 injective.

With these definitions at hand, the equations in ( 8) are reformulated via Lemma 3 (all Lemmas are given in Appendix A) as follows:

y i = x H q i q H i x = 2ℜ ν(q i ) H ν(x) - n j=1 (ν(q i )) jj (ν(x)) jj , i = 1, . . . , N.
Define the vectors a i ∈ C M , i = 1, . . . , N , with components given by (a i ) jk = 2 (ν(q i )) jk for 1 ≤ j < k ≤ n and (a i ) jj = (ν(q i )) jj , j = 1 . . . , n. Then, the equations above, linear wrt. to ν(x), can be rewritten as y i = ℜ a H i ν(x) , i = 1, . . . , N . Additionally define the binary matrices W j such that the vector W j ν(x) contains all the monomials of ν(x) including either x j or x j :

W j ν(x) = [x 1 x j , x 2 x j , . . . , x j-1 x j , x j x j , x j x j+1 , . . . , x j x n ] T ∈ C n .
Then, we have

x 0 = {W j ν(x)} n j=1 0 (10)
and problem ( 8) can be rewritten as the nonlinear group-sparse optimization program

min x∈C n {W j ν(x)} n j=1 0 (11) s.t. y = ℜ (Aν(x)) , where A = [a 1 , . . . , a N ] H .
Next, we relax this formulation by substituting v ∈ C M for ν(x), which yields

min v∈C M {W j v} n j=1 0 (12)
s.t. y = ℜ (Av)

v jj ∈ R + , j = 1, . . . , n,
where the last constraints ensure that the v jj 's estimating the modulus of the base variables are positive real numbers.

The following theorem shows that this relaxation can be used as a proxy to solve the original problem.

Theorem 3. If the solution v * to (12) is unique and yields x * = ν -1 (v * ) = 0 such that y i = |q H i x * | 2 , i = 1, . . . , N , then T (x * )
is the unique set of solutions of (8). Proof. Assume there is an x 0 = T (x * ) satisfying the constraints of (8) and at least as sparse as x * . Then y = ℜ (Aν(x 0 )), and, by using [START_REF] Gonsalves | Phase retrieval from modulus data[END_REF] and Lemma 4,

{W j ν(x 0 )} n j=1 0 = x 0 0 ≤ x * 0 ≤ {W j v * } n j=1 0 , which contradicts the fact that v * is the unique solution to (12) unless ν(x 0 ) = v * . But since x 0 / ∈ T (x * ), we have |x 0j | = |x * j | and thus x 0j x 0j = x * j x * j for some j ∈ {1, . . . , n}, which implies (ν(x 0 )) jj = (ν(x * )) jj = ν(ν -1 (v * )) jj .
Therefore, by using Lemma 5 with the assumption x * = ν -1 (v * ) = 0, there cannot be such an x 0 .

Convex relaxation

We now introduce a second level of relaxation by replacing the ℓ 0 -pseudo norm by a block-ℓ 1 norm. This leads to a convex relaxation in the form of a SOCP:

min v∈C M n j=1 W R j ℜ(v) + iW I j ℑ(v) (13)
s.t. y = ℜ (Av)

v jj ∈ R + , j = 1, . . . , n, where W R j = W j W R and W I j = W j W I with W = W R 0 0 W I a diagonal matrix of precom- pensating weights given by (W R ) i,i = ℜ(A i ) 2 and (W I ) i,i = ℑ(A i ) 2 . Theorem 4. Let x 0 be such that y i = x H 0 q i q H i x 0 , i = 1, . . . , N , v 0 = ν(x 0 ) and à = [ℜ(A), -ℑ(A)]. If the condition x 0 0 < 1 2 √ 2n 1 + 1 µ 2 ( Ã)
holds, then v 0 is the unique solution to [START_REF] Kohler | Source reconstruction from the modulus of the correlation function: a practical approach to the phase problem of optical coherence theory[END_REF].

Proof. The vector v 0 is the unique solution to [START_REF] Kohler | Source reconstruction from the modulus of the correlation function: a practical approach to the phase problem of optical coherence theory[END_REF] if the inequality

n j=1 W R j ℜ(v 0 + δ) + iW I j ℑ(v 0 + δ) > n j=1 W R j ℜ(v 0 ) + iW I j ℑ(v 0 ) holds for all δ ∈ C M such that ℜ(A(v 0 + δ)) = y, which implies the constraint ℜ (Aδ) = 0 on δ.
The inequality above can be rewritten as

j∈I0 W R j ℜ(δ) + iW I j ℑ(δ) + j / ∈I0 W R j ℜ(v 0 + δ) + iW I j ℑ(v 0 + δ) -W R j ℜ(v 0 ) + iW I j ℑ(v 0 ) > 0,
where

I 0 = {j ∈ {1, . . . , n} : W R j ℜ(v 0 ) = 0 ∧ W I j ℑ(v 0 ) = 0}. By the triangle inequality, a + b -a ≥ -b with a = W R j ℜ(v 0 ) + iW I j ℑ(v 0 ), this condition is met if j∈I0 W R j ℜ(δ) + iW I j ℑ(δ) - j / ∈I0 W R j ℜ(δ) + iW I j ℑ(δ) > 0 or n j=1 W R j ℜ(δ) + iW I j ℑ(δ) -2 j / ∈I0 W R j ℜ(δ) + iW I j ℑ(δ) > 0. ( 14
)
By defining G j as the set of indexes corresponding to nonzero columns of W j , Lemma 6 yields

W R j ℜ(δ) + iW I j ℑ(δ) 2 = i∈Gj (w R i ) 2 ℜ(δ i ) 2 + (w I i ) 2 ℑ(δ i ) 2 ≤ 2n µ 2 ( Ã) 1 + µ 2 ( Ã) W R ℜ(δ) + iW I ℑ(δ) 2 .
Due to the fact that k∈{1,...,n} G k = {1, . . . , M }, we also have

W R ℜ(δ) + iW I ℑ(δ) 2 = M i=1 (w R i ) 2 ℜ(δ i ) 2 + (w I i ) 2 ℑ(δ i ) 2 ≤ n k=1 i∈G k (w R i ) 2 ℜ(δ i ) 2 + (w I i ) 2 ℑ(δ i ) 2 = n k=1 W R k ℜ(δ) + iW I k ℑ(δ) 2 ≤ n k=1 W R k ℜ(δ) + iW I k ℑ(δ) 2 ,
which then leads to

W R j ℜ(δ) + iW I j ℑ(δ) 2 ≤ 2n µ 2 ( Ã) 1 + µ 2 ( Ã) n k=1 W R k ℜ(δ) + iW I k ℑ(δ) 2 .
Introducing this result in [START_REF] Lauer | Finding sparse solutions of systems of polynomial equations via group-sparsity optimization[END_REF] gives the condition

n j=1 W R j ℜ(δ) + iW I j ℑ(δ) -2(n -|I 0 |) µ( Ã) √ 2n 1 + µ 2 ( Ã) n k=1 W R k ℜ(δ) + iW I k ℑ(δ) > 0.
Finally, given that

|I 0 | = n -x 0 0 , this yields n j=1 W R j ℜ(δ) + iW I j ℑ(δ) -2 x 0 0 µ( Ã) √ 2n 1 + µ 2 ( Ã) n k=1 W R k ℜ(δ) + iW I k ℑ(δ) > 0. ( 15
)
or, for δ = 0,

x 0 0 < 1 + µ 2 ( Ã) 2µ( Ã) √ 2n ,
which can be rewritten as in the statement of the Theorem.

Corollary 2. Let x 0 be a feasible point of [START_REF] Fienup | Phase retrieval algorithms: a comparison[END_REF]. If the condition

x 0 0 < 1 2 √ 2n 1 + 1 µ 2 ( Ã)
holds, then T (x 0 ) is the unique set of solutions to the minimization problem [START_REF] Fienup | Phase retrieval algorithms: a comparison[END_REF] and it can be computed as T (x 0 ) = T (ν -1 (v)) with v the solution to [START_REF] Kohler | Source reconstruction from the modulus of the correlation function: a practical approach to the phase problem of optical coherence theory[END_REF].

Proof. Assume there exists another solution x 1 / ∈ T (x 0 ) to ( 8), and thus with x 1 0 ≤ x 0 0 . Then, Theorem 4 implies that both ν(x 1 ) and ν(x 0 ) are unique solutions to [START_REF] Kohler | Source reconstruction from the modulus of the correlation function: a practical approach to the phase problem of optical coherence theory[END_REF] and thus that ν(x 1 ) = ν(x 0 ) = v. But this contradicts Definition 2 implying ν(x 1 ) = ν(x 0 ) whenever x 1 / ∈ T (x 0 ). Therefore the assumption x 1 / ∈ T (x 0 ) cannot hold and T (x 0 ) is the unique set of solutions to (8), while ν -1 (v) = ν -1 (ν(x 0 )) ∈ T (x 0 ), which, by using ( 9), implies T (x 0 ) = T (ν -1 (v)).

Stable recovery in the presence of noise

Consider now the case where the measurements y are perturbed by an additive noise e ∈ R N of bounded ℓ 2 -norm, e 2 ≤ ε. Then, the equations in (1) and ( 8) are of the form y i = |q H i x| 2 + e i with the noise terms e i to be estimated together with the sparse signal. Note that in this context, multiple solutions with different noise vectors can be valid. Thus, we aim at stability results bounding the error on the estimates by a function of ε rather exact recovery ones. Details on the proposed method to achieve these goals are given below, first for real signals and then for complex ones.

Stability in the real case

In the noisy case with real data, the problem that we need to solve becomes min

x∈R n ,e∈R N x 0 (16) s.t. y i = |q T i x| 2 + e i , i = 1, . . . , N, e 2 ≤ ε,
where y i ∈ R and q i ∈ R n . Following the approach of Sect. 2 leads to a convex relaxation in the form of the SOCP min

v∈R M n j=1 W j W v 2 (17) s.t. y -Av 2 ≤ ε v jj ≥ 0, j = 1, . . . , n,
for which we have the following stability result.

Theorem 5. Let (x 0 , e 0 ) denote a solution to [START_REF] Ohlsson | On conditions for uniqueness in sparse phase retrieval[END_REF]. If the inequality

x 0 0 < 1 2n 2 (n + 1) 1 + 1 µ(A) (18)
holds, then the solution v to (17) must obey

W (v -ν(x 0 )) 2 2 ≤ 4nε 2 1 -µ(A)[2n 2 (n + 1) x 0 0 -1] . ( 19
)
If, in addition, ε = 0, then x = ±x 0 .

We omit the proof which is similar to the one of Theorem 6 in [START_REF] Lauer | Finding sparse solutions of systems of polynomial equations via group-sparsity optimization[END_REF], except for the last statement concluding on the stability of x, and which closely follows the one for the complex case of Theorem 6 to be detailed below.

Stability in the complex case

Consider now the complex variant of the problem perturbed by noise:

min x∈C n ,e∈R N x 0 (20) s.t. y i = |q H i x| 2 + e i , i = 1, . . . , N, e 2 ≤ ε,
where y i ∈ R and q i ∈ C n . The solution to this problem can be approached via the convex relaxation

min v∈C M n j=1 W R j ℜ(v) + iW I j ℑ(v) (21)
s.t. y -ℜ (Av) 2 ≤ ε v jj ∈ R + , j = 1, . . . , n.
As for the real case, we have a stability result for the estimation of ν(x).

Theorem 6. Let (x 0 , e 0 ) denote a solution to [START_REF] Shechtman | GESPAR: Efficient phase retrieval of sparse signals[END_REF]. If the inequality

x 0 0 < 1 2n 2 (n + 1) 1 + 1 µ( Ã) (22)
holds, then the solution v to (21) must obey

W R ℜ(v -ν(x 0 )) + iW I ℑ(v -ν(x 0 )) 2 ≤ 4nε 2 1 -µ( Ã)[2n 2 (n + 1) x 0 0 -1] . ( 23
)
If, in addition, ε = 0, then x = ν -1 (v) ∈ T (x 0 ).
Proof. Assume (20) has a solution (x 0 , e 0 ). Let define v 0 = ν(x 0 ) and δ = vv 0 . The proof follows a path similar to that of Theorem 3.1 in [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF], which was adapted in [START_REF] Lauer | Finding sparse solutions of systems of polynomial equations via group-sparsity optimization[END_REF] to the group-sparse setting and which is here further extended to the complex case. Due to the definition of v as a minimizer of (21), δ must satisfy either

n j=1 W R j ℜ(v 0 + δ) + iW I j ℑ(v 0 + δ) < n j=1 W R j ℜ(v 0 ) + iW I j ℑ(v 0 )
or δ = 0, in which case the statement is obvious. The inequality above can be rewritten as

j∈I0 W R j ℜ(δ) + iW I j ℑ(δ) + j / ∈I0 W R j ℜ(v 0 + δ) + iW I j ℑ(v 0 + δ) -W R j ℜ(v 0 ) + iW I j ℑ(v 0 ) < 0,
where

I 0 = {j ∈ {1, . . . , n} : W R j ℜ(v 0 ) = 0 ∧ W I j ℑ(v 0 ) = 0}. By the triangle inequality, a + b -a ≥ -b , this implies j∈I0 W R j ℜ(δ) + iW I j ℑ(δ) - j / ∈I0 W R j ℜ(δ) + iW I j ℑ(δ) < 0. ( 24
)
In addition, δ must satisfy the constraints in [START_REF] Shechtman | Sparsity based sub-wavelength imaging with partially incoherent light via quadratic compressed sensing[END_REF] as

ℜ(A(v 0 + δ)) -y 2 ≤ ε,
in which y can be replaced by ℜ(Av 0 ) + e 0 , leading to

ℜ(Aδ) -e 0 2 ≤ ε. Using a 2 ≤ a -b 2 + b 2 , this implies ℜ(Aδ) 2 ≤ 2ε, which further gives 4ε 2 ≥ ℜ(Aδ) 2 2 = Ãδ 2 2 = Ã W -1 W δ 2 2 = ( W δ) T W -1 ÃT Ã W -1 ( W δ) = W δ 2 2 + ( W δ) T ( W -1 ÃT Ã W -1 -I)( W δ) ≥ W δ 2 2 -( W δ) T ( W -1 ÃT Ã W -1 -I)( W δ) ≥ W δ 2 2 -| W δ| T | W -1 ÃT Ã W -1 -I|| W δ| ≥ W δ 2 2 -µ( Ã)( W δ 2 1 -W δ 2 2 ) = (1 + µ( Ã)) W δ 2 2 -µ( Ã) W δ 2 1 ( 25
)
where we used W -1 W = I and the fact that the diagonal entries of | W -1 ÃT Ã W -1 -I| are zeros while off-diagonal entries are bounded from above by µ( Ã).

Due to W R j ℜ(δ) + iW I j ℑ(δ) being a vector with a subset of entries from W R ℜ(δ)

+ iW I ℑ(δ), we have W δ 2 2 = W R ℜ(δ) + iW I ℑ(δ) 2 ≥ W R j ℜ(δ) + iW I j ℑ(δ) 2 , j = 1, . . . , n, and thus W δ 2 2 ≥ 1 n n j=1 W R j ℜ(δ) + iW I j ℑ(δ) 2 . ( 26
)
Since the groups defined by the W j 's overlap,

W δ 2 = W R ℜ(δ) + iW I ℑ(δ) ≤ n j=1 W R j ℜ(δ) + iW I j ℑ(δ)
, and the squared ℓ 1 -norm in (25) can be bounded by

W δ 2 1 ≤ M W δ 2 2 ≤ M   n j=1 W R j ℜ(δ) + iW I j ℑ(δ)   2 . ( 27
)
Introducing the bounds (26)-( 27) in (25) yields

1 + µ( Ã) n n j=1 W R j ℜ(δ) + iW I j ℑ(δ) 2 -µ( Ã)M   n j=1 W R j ℜ(δ) + iW I j ℑ(δ)   2 ≤ 4ε 2 . ( 28
)
We will now use this inequality to derive an upper bound on

n j=1 W R j ℜ(δ) + iW I j ℑ(δ) 2 , which will also apply to W δ 2 2 ≤ n j=1 W R j ℜ(δ) + iW I j ℑ(δ) 2 = n j=1 ( W R j ℜ(δ) 2 2 + W I j ℑ(δ) 2
2 ), since the groups overlap and the squared components of W δ are summed multiple times in the right-hand side. To derive the upper bound, we first introduce a few notations:

a = j∈I0 W R j ℜ(δ) + iW I j ℑ(δ) , b = j / ∈I0 W R j ℜ(δ) + iW I j ℑ(δ) , and
c 0 = j∈I0 W R j ℜ(δ) + iW I j ℑ(δ) 2 j∈I0 W R j ℜ(δ) + iW I j ℑ(δ) 2 ∈ 1 |I 0 | , 1 , c 1 = j / ∈I0 W R j ℜ(δ) + iW I j ℑ(δ) 2 j / ∈I0 W R j ℜ(δ) + iW I j ℑ(δ) 2 ∈ 1 n -|I 0 | , 1 ,
where the box bounds are obtained by classical relations between the ℓ 1 and ℓ 2 norms2 . With these notations, the term to bound is rewritten as

n j=1 W R j ℜ(δ) + iW I j ℑ(δ) 2 = c 0 a 2 + c 1 b 2 ,
while the inequality (28) becomes

1 + µ( Ã) n (c 0 a 2 + c 1 b 2 ) -µ( Ã)M (a + b) 2 ≤ 4ε 2 .
We further reformulate this constraint by letting a = ρb:

1 + µ( Ã) n (c 0 ρ 2 + c 1 )b 2 -µ( Ã)M (1 + ρ) 2 b 2 ≤ 4ε 2 . ( 29
) Let γ = (1 + ρ) 2 /(c 0 ρ 2 + c 1 )
. Due to (24), we have a < b and thus ρ ∈ [0, 1), which, together with the bounds on c 0 and c 1 , gives the constraints 1

≤ γ ≤ 4(n -|I 0 |). By setting V = (c 0 ρ 2 + c 1 )b 2 , ( 29
) is rewritten as 1 + µ( Ã) n V -µ( Ã)M γV ≤ 4ε 2 , where 1 + µ( Ã) n -µ( Ã)M γ ≥ 1 + µ( Ã) n -4(n -|I 0 |)µ( Ã)M > 0, since γ ≤ 4(n -|I 0 |
) and the positivity is ensured by the condition [START_REF] Waldspurger | Phase recovery, maxcut and complex semidefinite programming[END_REF] and the fact that

x 0 0 = n -|I 0 |. Thus, W δ 2 2 ≤ n j=1 W R j ℜ(δ) + iW I j ℑ(δ) 2 = V ≤ 4nε 2 1 + µ( Ã) -4µ( Ã)nM x 0 0 , which proves the stability result in (23) since W δ 2 2 = W R ℜ(δ) + iW I ℑ(δ) 2 .
To conclude in the case ε = 0, it remains to see that [START_REF] Walther | The question of phase retrieval in optics[END_REF] implies v = ν(x 0 ) and that Definition 3 ensures ν -1 (ν(x 0 )) ∈ T (x 0 ).

Complex data, but real solutions

Consider now the following problem:

min x∈R n x 0 (30) s.t. y i = |q H i x| 2 , i = 1, . . . , N,
where the signal x ∈ R n and measurements y i ∈ R are assumed to be real while the vectors q i ∈ C n can be complex. In this case, v = ν(x) is also a real vector and solutions can be approximated via a dedicated version of (13):

min v∈R M n j=1 W R j v 2 (31) s.t. y = ℜ (Av) = ℜ (A) v v jj ≥ 0, j = 1, . . . , n,
where A ∈ C N ×M is defined as in Sect. 3.1 and W R j = W j W R with precompensating weights given by (W R ) i,i = ℜ(A i ) 2 .

For this particular case, Theorem 7 below is similar in spirit to Theorem 4, but allows us to gain a √ 2 factor by using Lemma 7 instead of Lemma 6 in order to take into account that δ belongs to R M (detailed proof given in Appendix B.2). This results in a bound on x 0 0 similar to the one in Theorem 2 for the real case and based on the mutual coherence of the real part of A. Since µ(ℜ(A)) ≤ µ( Ã), this also improves (relaxes) the bound compared with the one of Theorem 4. Theorem 7. Let x 0 ∈ R n be such that y i = x T 0 q i q H i x 0 , i = 1, . . . , N , and v 0 = ν(x 0 ). If the condition

x 0 0 < 1 2 √ n 1 + 1 µ 2 (ℜ(A))
holds, then v 0 is the unique solution to (31).

In a typical instance of Problem (30), the measurements y correspond to the power spectrum of a real signal. However, in this case, Theorem 7 does not apply since the solution is known not to be unique due to invariances of the measurements. The following subsections first describe a practical technique to deal with such invariances and then focus on the power spectrum case.

Invariances

Consider the case where the measurements are invariant to some transformation of the signal. A typical example is when x ∈ R n and {q i ∈ C n } forms a Fourier basis. Then, the equations y i = |q H i x| 2 are invariant to circular shifts of the components of x. This is problematic since shifted versions of x lead to shuffled 3 versions of ν(x). Thus, multiple shuffled v's satisfy the linearized constraints, y i = ℜ(A)v, and so does any convex combination of them. These convex combinations need not be sparse as they combine vectors with different sparsity patterns, but lead to lower or equal values of the convex cost function of (31).

To circumvent this issue, we need to linearize the constraints by a shift-invariant transformation φ, such that φ(x 0 ) = φ(x 1 ) for x 0 and x 1 two shifted versions of the same vector. To lead to an effective estimation method, the transformation φ must also retain sparsity in a sense similar to ν.

Consider the transformation φ : R n → R M defined by 4

φ(x) = ν(shift(x, 1 -k)), with k = arg max i∈{1,...,n} |x i |, (32)
where shift(•, k) stands for the k-steps circular shift operator. The transformation φ in (32) first shifts the vector x such that the first component is the one with maximal magnitude. This results in a shift-invariant transformation which retains sparsity and the linearization ability via the mapping ν. The linearized constraints remain the same, i.e., y i = ℜ(A)φ, but another constraint must be added to (31) in order to account for the shift-invariance. More precisely, the definition of φ ensures that (φ(x)) 11 ≥ (φ(x)) jj , j = 2, . . . , n.

Thus, a valid vector φ estimating φ(x 0 ) can be found by solving φ = arg min

φ∈R M n j=2 W R j φ 2 (33) s.t. y = ℜ (A) φ φ 11 ≥ φ jj ≥ 0, j = 2, . . . , n.
3 Circular shifts of x lead to rearrangements of the components of ν(x) which are not exactly circular shifts. For example, with x 0 = [1, 1, 0, 0] T and x 1 = [0, 1, 1, 0] T , we have ν(x 0 ) = [1, 1, 0, 0, 1, 0, 0, 0, 0, 0] T and ν(x 1 ) = [0, 0, 0, 0, 1, 1, 0, 1, 0, 0] T . 4 In the case where the argmax in (32) is not a singleton, k is arbitrarily set to the minimum of the indexes in the argmax and φ cannot be shift-invariant. Assumptions regarding this issue will be made clear in Proposition 1 and Theorem 8 below.

Note that the the cost function does not involve the first group of variables since x 1 is assumed to be nonzero. In comparison with (31), this formulation is still convex but cannot have multiple solutions that are shifted/shuffled versions of one another.

Finally, the set of shifted solutions {x k } n k=1 to (30) is approximated by the set

φ -1 ( φ) = shift(ν -1 ( φ), k) n k=1
.

Following a similar approach, we can additionally take into account reflections by defining

               φ(x) = ν(x 2 ) x 2 = ϕ(x) = x 1 , if n/2 i=2 |x 1i | 2 ≥ n i=2+ n 2 |x 1i | 2
shift(reflection(x 1 ), 1), otherwise

x 1 = shift(x, 1 -k) k = arg max i∈{1,...,n} |x i |, (34)
where we assume n to be even and reflection(•) is the reflection operator defined by reflection(x) = [x n , x n-1 , . . . , x 2 , x 1 ] T . In plain words, the transformation φ in (34) first shifts the vector x such that the first component is the one with maximal magnitude. Then, it applies a centered reflection to the shifted x, named x 1 , only if the sum of squares over the first entries of x 1 (without the first one) is smaller than the one over the last ones. If this is the case, the result is shifted again to recover the first component of x 2 with maximal magnitude. Finally, the Veronese map is applied to x 2 to give φ(x).

Invariance of φ = ν • ϕ to circular shifts and reflections is implied by the invariance of ϕ given in the Proposition below (proof in Appendix B.3).

Proposition 1. For all x ∈ C n such that arg max i∈{1,...,n} |x i | = 1, the following statements hold for the transformation ϕ defined in (34):

1. ϕ is idempotent, i.e., ϕ • ϕ(x) = ϕ(x);
and, for all x additionally satisfying

n/2 i=2 |x 1i | 2 = n i=2+ n 2 |x 1i | 2 with x 1 = shift(x, arg max i∈{1,...,n} |x i |),
2. ϕ is shift-invariant, i.e., ∀s ∈ Z, ϕ(x) = ϕ(shift(x, s));

3. ϕ is reflection-invariant, i.e., ϕ(x) = ϕ(reflection(x)).

Note that statements 2 and 3 imply that ϕ is invariant to any combination of shifts and reflections.

Proposition 1 also shows the idempotence of ϕ, which of course does not transfer to φ directly but which is however very useful. Indeed, this allows us to test if a candidate x 2 , supposed to be the result of ϕ applied to some vector, is consistent with the definition of ϕ as x 2 ? = ϕ(x 2 ), which can be checked via simple inequalities. Since these inequalities only involve the (squared) magnitude of the entries in the vector, they can also be easily embedded as linear constraints in (33) to compute an estimate φ that is consistent with the transformation (34). This yields the convex program φ = arg min

φ∈R M n j=2 W R j φ 2 (35) s.t. y = ℜ (A) φ (data fitting) φ 11 ≥ φ jj ≥ 0, j = 2, . . . , n (shift-invariance) n/2 i=2 φ ii ≥ n i=2+ n 2 φ ii (reflection-invariance),
where φ jj estimates |x j | 2 . The advantage of using ( 35) is that it has a single solution independently of the number of shifted/reflected solutions to (30), in the sense of the next theorem.

Theorem 8. If Problem (30) has a unique set of shifted/reflected solutions S(x 0 ) = {shift(x 0 , k) : k ∈ {1, . . . , n}} ∪ {shift(reflection(x 0 ), k) : k ∈ {1, . . . , n}}, and if arg max i∈{1,...,n} |x 0i | = 1, then there is exactly one vector ν(x) with x ∈ S(x 0 ) in the feasible set of (35).

Proof. Since all x ∈ S(x 0 ) are solutions to (30), they satisfy the first constraint in (35) as y = ℜ (A) ν(x). By Proposition 1, ϕ is constant over the set S(x 0 ) and equal to one of the vectors from this set, say x 0 , for which we have in particular x 0 = ϕ(x 0 ). Then, ν(x 0 ) = ν • ϕ(x 0 ) = φ(x 0 ) and thus ν(x 0 ) satisfies all the constraints of (35) and is a feasible point. Now take an x ∈ S(x 0 ) with x = x 0 . Then, either x = shift(x 0 , k) or x = shift(reflection(x 0 ), k) for some k ∈ {1, . . . , n}. This implies arg max i (ν(x)

) ii = arg max i |x i | = arg max i |x 0i | = arg max i (ν(x 0 )) ii .
Since x 0 satisfies the constraints in (35), we have arg max i (ν(x 0 )) ii = 1 and thus arg max i (ν(x)) ii = 1, which shows that ν(x) violates the second constraint and is not a feasible point of (35).

Illustrative example. Consider the vectors (x 1 and x 2 are not related to the notations of (34))

x 1 = [1, 2, 3, 4, 0, 0] T , x 2 = [4, 0, 0, 1, 2, 3] T , x 3 = [0, 0, 4, 3, 2, 1] T , x 4 = [2, 1, 0, 0, 4, 3] T , x 5 = [4, 3, 2, 1, 0, 0] T ,
which are all obtained by shifts and reflections of the same vector. They all lead to the same φ(x i ) = ν(x 5 ). Indeed, the transformation φ first applies the required shift and reflect operations to map x i to x 5 which has a first component with the largest magnitude and the largest half-sum over its first entries. Then, φ computes the Veronese map of x 5 . Also note that all the ν(x i ) are feasible with respect to y = ℜ (A) ν(x i ), but only ν(x 5 ) satisfies the additional constraints implementing shift and reflection invariance in (35). Therefore, we can obtain φ = ν(x 5 ) as the unique solution to (35).

Support recovery from the power spectrum

When the measurements correspond to the squared magnitude of the Fourier transform of the signal, i.e., its (squared) power spectrum, there is another issue beside shift/reflection-invariances. Even with the correct support S = supp(x 0 ), the linear system y = ℜ(A S )φ S , where A S is the submatrix of A with the corresponding columns, is under-determined. More precisely, for all i and all j, |q ij | = 1, which implies (a i ) jj = (ν(q i )) jj = 1 and A jj = 1. Thus, even when limited to the support S, A S has |S| similar columns and is rank-defficient.

Therefore, in this case the proposed approach cannot exactly recover x 0 and (35) is used to estimate the support S. Indeed, knowing the correct support can be useful for other methods dedicated to the classical phase retrieval problem [START_REF] Fienup | Phase retrieval algorithms: a comparison[END_REF].

Though this case seems unfavorable, if the number of measurements satisfies N ≥ 2n -1, the autocorreleation of x 0 (with zero padding),

r k = n i=1 x i x i+k , k = -n + 1, . . . , n -1,
can be computed via the inverse Fourier transform of y and thus is available. This information can be included as linear constraints on the components φ i(i+k) of φ estimating the cross products x i x i+k to drive the solution towards a satisfactory one.

In addition, following [START_REF] Jaganathan | Recovery of sparse 1-d signals from the magnitudes of their Fourier transform[END_REF][START_REF] Shechtman | GESPAR: Efficient phase retrieval of sparse signals[END_REF], this can be used to restrict the support of the solution as follows. First, note that with N ≥ 2n -1, the measurements are not invariant to circular shifts of x but of x with zero-padding. Thus, we cannot assume |x 1 | ≥ |x j |, j = 2, . . . , n. However, we can fix x 1 = 0. Then, for all base variable index j ∈ {1, . . . , n},

r k = 0, k = j -1, . . . , n -1 ⇒ x j = 0. (36)
Thus, the corresponding variables in φ can be set to zero, i.e., W j φ = 0. Moreover, let j m be the minimum of the j's satisfying the condition in (36), then x jm-1 = 0, which can be favored by removing the corresponding term in the cost function. Combining all the information on the solution in a convex program leads to the final formulation φ = arg min

φ∈R M jm-2 j=2 W R j φ 2 (37)
s.t. y = ℜ (A) φ (data fitting)

φ jj ≥ 0, j = 1, . . . , j m -1 (structural knowledge) n/2 i=2 φ ii ≥ n i=2+ n 2 φ ii (reflection-invariance) n-k i=1 φ i(i+k) = r k , k = 0, . . . , n -1 (autocorrelation)
W j φ = 0, j = j m , . . . , n (restricted support).

A noise-tolerant version of (37) is obtained by replacing the constraint y = ℜ (A) φ by y -ℜ (A) φ 2 ≤ ε. This modification also applies to (31) and (35).

Experiments

For the experiments, we consider two variants of the proposed approach: the convex method described in details in the previous sections and the greedy method for solving the group-sparse optimization problems ( 5) and [START_REF] Jaganathan | Recovery of sparse 1-d signals from the magnitudes of their Fourier transform[END_REF]. Implementation details for these two methods are as described in [START_REF] Lauer | Finding sparse solutions of systems of polynomial equations via group-sparsity optimization[END_REF], with slight modifications to handle complex variables. In particular, the convex method uses the iterative reweighting of [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF] adapted to the group-sparse setting to enhance the sparsity of the solution. The greedy method starts with an empty support and, at each iteration, adds the group of variables in v corresponding to the base variable x j that results in the best approximation of y.

For complex signals, x 0 , we measure the relative error corresponding to the normalized distance between the estimate x and the set T (x 0 ), i.e., min x∈T (x0) x-x / x . Exact recovery is detected when this error is smaller than 10 -6 x 0 .

Exact recovery in the noiseless case. We estimate the probability of exact recovery at various sparsity levels, x 0 0 , in the complex case where n = 20 and N = 50. For each sparsity level, the probability is estimated as the percentage of successful trials over a Monte Carlo experiment with 100 trials. In each trial, N complex vectors q i ∈ C n are drawn from a zero-mean Gaussian distribution of unit variance to measure a random signal x 0 ∈ C n with x 0 0 nonzero entries at random locations whose values are drawn from a zero-mean Gaussian distribution of unit variance. Results shown in Fig. 1 indicate that the proposed approach can exactly recover sufficiently sparse signals with high probability.

Similar experiments are performed to evaluate the influence of the number of measurements N on the probability of exact recovery. Results reported in the right plot of Fig. 1 show that, with a sparsity level of x 0 0 = 4, the convex method requires more measurements for perfect recovery than the greedy strategy. However, N ≈ x 0 3 0 measurements are already sufficient to exactly recover x 0 in all trials.

Stable recovery in the noisy case. We now consider the noisy case where y i = |q H i x 0 | 2 + e i , i = 1, . . . , N , and e 2 ≤ ε. Over 100 trials, Figure 2 reports the mean relative error and the rate of successful support recovery for N = 50 and ε = 3. These results show that sufficiently sparse signals can be accurately estimated from noisy quadratic measurements and that the probability of support recovery in this case follows a curve similar to the one of the probability of exact recovery in the noiseless case. This means that the methods are robust to noise regarding the recovery of the correct support. Support estimation from the power spectrum. We now test if the convex method is robust to invariances, such as the ones discussed in Sect. 5.2, when estimating the support. In particular, we start with a setting in which x ∈ R n , N = n = 20 and q H i , i = 1, . . . , N , are the rows of the n-point Fourier matrix. In this case, the autocorrelation cannot be computed and the convex formulation ( 35) is used to estimate the support. Then, we perform similar experiments but with oversampling (N = 2n), thus allowing for the computation of the autocorrelation and the use of (37). Results shown in Fig. 3 indicate that by using (35) we can recover the correct support for sufficiently sparse signals without oversampling, i.e., with as many measurements as unknowns, while using more information extracted from the autocorrelation of the signal only slightly helps to recover larger supports.

Conclusions

The paper proposed a new approach to phase retrieval of sparse signals. This approach is based on a group-sparse optimization formulation of the problem with linearized constraints. Exact and stable recovery results were shown for a convex relaxation of this formulation both in the real and complex case. Invariances to circular shifts and reflections that are common in phase retrieval problems were also discussed and a practical technique was given to prevent these from breaking the sparsity of the solution.

Future work will focus on deriving theoretical guarantees for the invariance issue. Another direction of research concerns the analysis of greedy algorithms for the group-sparse optimization problem, which proved as valuable as the convex relaxations in experiments for measurements without invariance. How to deal with invariances in these methods will also be investigated.

A Lemmas

Lemma 1. Let ν and ν -1 be defined as in Section 2. For all v ∈ R M such that ν

-1 (v) = 0 and v ′ ∈ R M , if ν(ν -1 (v)) jj = v ′ jj for some j ∈ {1, . . . , n}, then v = v ′ .
Proof. Assume that ν -1 (v) = 0 and let i = min j j, s.t. v jj > 0. Then, ∀j ∈ {1, . . . , n}, we have

ν(ν -1 (v)) jj = ν -1 (v) 2 j = v 2 ji v ii = v jj ,
where the last equality follows from the definition of ν

-1 . Therefore, if ν(ν -1 (v)) jj = v jj = v ′ jj for some j ∈ {1, . . . , n}, we obtain v = v ′ . Lemma 2. Let A = [A 1 , . . . , A M ] be an N × M real matrix with mutual coherence µ(A) as in Definition 1. Let W be the M × M -diagonal matrix of entries w i = A i 2 .
Then, for all δ ∈ Ker(A) and i ∈ {1, . . . , M }, the bound

w 2 i δ 2 i ≤ µ 2 (A) 1 + µ 2 (A) W δ 2 2 (38) holds.
Proof. See Lemma 2 in [START_REF] Lauer | Finding sparse solutions of systems of polynomial equations via group-sparsity optimization[END_REF].

Lemma 3. Let the complex Veronese map ν be as in Definition 2. Then, for all q ∈ C n and x ∈ C n , the following equality holds:

x H qq H x = 2ℜ ν(q) H ν(x) - n j=1 (ν(q)) jj (ν(x)) jj .
Proof.

x

H qq H x = n k=1 x k q k   n j=1 q j x j   = n j=1 q j x j n k=1 q k x k = n j=1 n k=1 q j q k x j x k = n j=1 q j q j x j x j + n j=1 k =j q j q k x j x k = n j=1 q j q j x j x j + n j=1   k<j q j q k x j x k + k>j q j q k x j x k   = n j=1 q j q j x j x j + n j=1 k>j (q k q j x k x j + q j q k x j x k ) = n j=1
q j q j x j x j + n j=1 k>j q j q k x j x k + q j q k x j x k At this point, we use the fact that z + z = 2ℜ(z), which yields

x H qq H x = n j=1 q j q j x j x j + 2 n j=1 k>j ℜ (q j q k x j x k ) = n j=1 (ν(q)) jj (ν(x)) jj + 2 n j=1 k>j ℜ ν(q) jk (ν(x)) jk with W = W R 0 0
W I a diagonal matrix of entries ( W ) i,i = Ãi 2 , and

(w I i ) 2 ℑ(δ i ) 2 ≤ µ 2 ( Ã) 1 + µ 2 ( Ã) W δ 2 2 . These inequalities lead to (w R i ) 2 ℜ(δ i ) 2 + (w I i ) 2 ℑ(δ i ) 2 ≤ 2µ 2 ( Ã) 1 + µ 2 ( Ã) W δ 2 2 , where W δ 2 2 = W R ℜ(δ) 2 2 + W I ℑ(δ) 2 2 = W R ℜ(δ) + iW I ℑ(δ) 2 .
Lemma 7. For all δ ∈ R M such that ℜ(A)δ = 0, the inequality

(w R i ) 2 δ 2 i ≤ µ 2 (ℜ(A)) 1 + µ 2 (ℜ(A)) W R δ 2 2 ,
where W R is a diagonal matrix of entries

(W R ) i,i = w R i = ℜ(A i ) 2 .
Proof. This is a direct consequence of Lemma 2 applied to ℜ(A) and δ ∈ Ker(ℜ(A)).

B Proofs B.1 Proof of Theorem 2

This proof is similar to the one of Theorem 3 in [START_REF] Lauer | Finding sparse solutions of systems of polynomial equations via group-sparsity optimization[END_REF].

Proof. The vector v 0 is the unique solution to (6) if the inequality

n j=1 W j W (v 0 + δ) 2 > n j=1 W j W v 0 2
holds for all δ = 0 satisfying Aδ = 0. The inequality above can be rewritten as j∈I0

W j W δ 2 + j / ∈I0 W j W (v 0 + δ) 2 -W j W v 0 2 > 0, where I 0 = {j ∈ {1, . . . , n} : W j W v 0 = 0}. By the triangle inequality, a + b 2 -a 2 ≥ -b 2 , this condition is met if j∈I0 W j W δ 2 - j / ∈I0 W j W δ 2 > 0 or n j=1 W j W δ 2 -2 j / ∈I0 W j W δ 2 > 0. ( 40
)
By defining G j as the set of indexes corresponding to nonzero columns of W j , Lemma 2 yields

W j W δ 2 2 = i∈Gj w 2 i δ 2 i ≤ n µ 2 (A) 1 + µ 2 (A) W δ 2 2 ,
Due to the fact that k∈{1,...,n} G k = {1, . . . , M }, we also have

W δ 2 2 = M i=1 w 2 i δ 2 i ≤ n k=1 i∈G k w 2 i δ 2 i = n k=1 W k W δ 2 2 ≤ n k=1 W k W δ 2 2 ,
which then leads to

W j W δ 2 2 ≤ n µ 2 (A) 1 + µ 2 (A) n k=1 W k W δ 2 2 .
Introducing this result in (40) gives the condition

n j=1 W j W δ 2 -2(n -|I 0 |) µ(A) √ n 1 + µ 2 (A) n k=1 W k W δ 2 > 0.
Finally, given that

|I 0 | = n -{W j v 0 } n j=1 0 = n -x 0 0 , this yields n j=1 W j W δ 2 -2 x 0 0 µ(A) √ n 1 + µ 2 (A) n k=1 W k W δ 2 > 0.
or, after rearranging the terms,

x 0 0 < 1 + µ 2 (A) 2µ(A) √ n ,
which can be rewritten as in the statement of the Theorem.

B.2 Proof of Theorem 7

This proof is very similar to the ones of Theorems 2 and 4.

Proof. The vector v 0 is the unique solution to (31) if the inequality

n j=1 W R j (v 0 + δ) 2 > n j=1 W R j v 0 2
holds for all δ ∈ R M such that ℜ(A)(v 0 + δ) = y, which implies the constraint ℜ (A) δ = 0 on δ.

The inequality above can be rewritten as

j∈I0 W R j δ 2 + j / ∈I0 W R j (v 0 + δ) 2 -W R j v 0 2 > 0,
where I 0 = {j ∈ {1, . . . , n} : W R j v 0 = 0}. By the triangle inequality, a + ba ≥b with a = W R j v 0 , this condition is met if

n j=1 W R j δ 2 -2 j / ∈I0 W R j δ 2 > 0. ( 41
)
By defining G j as the set of indexes corresponding to nonzero columns of W j , Lemma 7 yields

W R j δ 2 2 = i∈Gj (w R i ) 2 δ 2 i ≤ n µ 2 (ℜ(A)) 1 + µ 2 (ℜ(A)) W R δ 2 2 .
Due to the fact that k∈{1,...,n} G k = {1, . . . , M }, we also have Finally, given that |I 0 | = nx 0 0 , this yields , for δ = 0,

W R δ 2 2 = M i=1 (w R i ) 2 δ 2 i ≤ n k=1 i∈G k (w R i ) 2 δ 2 i = n k=1 W R k δ 2 2 ≤ n k=1 W R k δ 2
x 0 0 < 1 + µ 2 (ℜ(A)) 2µ(ℜ(A)) √ n ,
which can be rewritten as in the statement of the Theorem.

B.3 Proof of Proposition 1

Proof. To prove Statement 1, note that the operation x 2 = shift(reflection(x 1 ), 1) is equivalent to

x 2 =
|x 2i | 2 = n/2-2 i=0 |x 1(n-i) | 2 = n i=2+ n 2 |x 1i | 2 > n/2 i=2 |x 1i | 2 = n/2-2 i=0 |x 2(n-i) | 2 = n i=2+ n 2 |x 2i | 2 ,
where the inequality holds whenever this operation is performed in (34) (i.e., if this is not the case, then x 2 = x 1 ). Therefore, for all x, we obtain a vector x 2 = ϕ(x) such that and ϕ applies a reflection to x if and only if it does not apply one to x ′ (the inequalities above are strict by assumption). Thus, in the case x is reflected by ϕ, we have

Figure 1 :

 1 Figure1: Estimated probability of exact recovery in the noiseless case for the convex relaxation (13) (plain line) and the greedy method applied to (12) (dashed line) versus the number of nonzeros (left) and the number of measurements (right).

Figure 2 :

 2 Figure 2: Mean relative error (left) and estimated probability of support recovery (right) in the noisy case for the convex relaxation (13) (plain line) and the greedy method applied to (12) (dashed line).

Figure 3 :

 3 Figure3: Estimated probability of support recovery from the magnitude of the Fourier transform when using (35) for N = n (plain line) and when using (37) for N = 2n (dashed line) versus the sparsity level.

2 ,

 2 which then leads toW R j δ 2 2 ≤ n µ 2 (ℜ(A)) 1 + µ 2 (ℜ(A))

n/ 2 i=2 |x 2i | 2 ≥ n i=2+ n 2 |x 2i | 2 .

 222 Since arg max i∈{1,...,n} |x 2i | = 1 and shift(x 2 , 1 -arg max i∈{1,...,n} |x 2i |) = x 2 , this implies that ϕ(x 2 ) = x 2 , i.e., ϕ(ϕ(x)) = ϕ(x).Statement 2 is easily seen from the fact that x 1 does not change when x is shifted.To prove Statement 3, let x ′ = reflection(x) = [x n , x n-1 , . . . , x 1 ] T , k = arg max i∈{1,...,n} |x i |, and k ′ = arg max i∈{1,...,n} |x ′ i |. Then, k ′ = n -k + 1 andx ′ 1 = shift(x ′ , 1 -k ′ ) = shift(x ′ , k -n) = shift(x ′ , k) = [x ′ n-k+1 , x ′ n-k+2 , . . . , x ′ n , x ′ 1 , x ′ 2 , . . . , x ′ n-k ] T = [x k , x k-1 , . . . , x 1 , x n , x n-1 , . . . , x k+1 ] T . Define x 1 = shift(x, 1 -k) = [x k , x k+1 , . . . , x n , x 1 , x 2 , . . . , x k-1 ] T . If k < n/2, we have

ϕ

  (x ′ ) = shift(reflection(x), k) = reflection(shift(x, n -k)) and ϕ(x) = shift(reflection(shift(x, 1 -k)), 1) = reflection(shift(shift(x, 1 -k), n -1)) = reflection(shift(x, n -k)) = ϕ(x ′ ),while in the case x is not reflected by ϕ, we haveϕ(x) = shift(x, 1 -k)and ϕ(x ′ ) = shift(reflection(shift(reflection(x), k)), 1) = reflection(shift(shift(reflection(x), k), n -1)) = reflection(shift(reflection(x), n + k -1)) = reflection(reflection(shift(x, -k + 1))) = shift(x, 1 -k) = ϕ(x).

  x 11 reflection(x 1 ) , x1 = [x 12 , x 13 , . . . , x 1n ] T ,where the reflection is centered on x1(n/2) = x 1(1+n/2) . Thus,x 2 = [x 11 , x 1n , x1(n-1) , . . . , x 12 ] T

	and
	n/2
	i=2

More precisely, the entries of W j are given by (W j ) k,l = δ l,jk , l = 1, . . . , M , k = 1, . . . , n, where δ is the Kronecker delta and jk is an index as in[START_REF] Candès | Phase retrieval via matrix completion[END_REF].

Let u ∈ R |I 0 | with u j = W R j ℜ(δ) + iW I j ℑ(δ) . Then, c 0 = ( u 2 / u 1 )2 and the bounds are obtained by the classical relation ∀u∈ R k , u 2 ≤ u 1 ≤ √ k u 2 .

Since (ν(q)) jj (ν(x)) jj = q j q j x j x j = |q j x j | 2 is a real number, it can be introduced in the second sum as

where the last equality is due to

Lemma 4. Let ν -1 be as in Definition 3 and the matrices W j as in Sect. 3.1. Then, for all

Proof.

On the other hand, v jj belongs to a single group of variables generated by the W j , i.e., (W k

However, note that the converse is not true: we can have ν -1 (v) 0 = {W j v} n j=1 0 , e.g., when v jj = 0, j = 1, . . . , n, and v ij = 0 for some i and j. Lemma 5. Let ν and ν -1 be defined as in Definitions 2 and 3. For all v ∈ C M such that

Proof. Assume that ν -1 (v) = 0 and let i = min j j, s.t. ℜ(v jj ) > 0, ℑ(v jj ) = 0. Then, ∀j ∈ {1, . . . , n}, we have

where the last equality follows from Definition 3. Therefore, if ν(ν -1 (v)) jj = v jj = v ′ jj for some j ∈ {1, . . . , n}, we obtain v = v ′ . Lemma 6. For all δ ∈ C M such that ℜ(Aδ) = 0, the inequality

Proof. We rewrite the assumption as

where δ = [ℜ(δ) T , ℑ(δ) T ] T ∈ R 2M . Then, we can use Lemma 2 to bound the entries in δ as